au_qmr_ilu_s(tr,X,i)

function Y = au_qmr_ilu_s(tr,X,i)
%
%  Solves shifted linear systems with the matrix A or its transpose A' by
%  the QMR iteration with ILU preconditioning:
%
%  for tr = 'N':
%
%    Y = inv(A+p(i)*I)*X;
%
%  for tr = 'T':
%
%    Y = inv(A.'+p(i)*I)*X.
%
%  The incomplete LU factors of the ILU preconditioner for A+p(i)*I are 
%  provided as global data if mc has been set to 'C' when calling 
%  'au_qmr_ilu_m_i'. In this case, this data must be generated by calling
%  'au_qmr_ilu_s_i' before calling this routine! If mc has been set to 'M', 
%  the incomplete factorizations are computed within this routine to save
%  memory. However, the computational cost is increased in this case.
%
%  Calling sequence:
%
%    Y = au_qmr_ilu_s(tr,X,i)
%
%  Input:
%
%    tr        (= 'N' or 'T') determines whether shifted systems with 
%              A or A' should be solved;
%    X         a matrix of proper size;
%    i         the index of the shift parameter.
%
%  Output:
%
%    Y         the resulting solution matrix.
%  
%
%   LYAPACK 1.0 (Thilo Penzl, August 1999)

if nargin~=3
  error('Wrong number of input arguments.');
end

global LP_A LP_P LP_MC LP_MAX_IT_QMR LP_TOL_QMR LP_TOL_ILU LP_INFO_QMR

if ~length(LP_A) | ~length(LP_P) | ~length(LP_MC) | ~length(LP_MAX_IT_QMR) ...
  | ~length(LP_TOL_QMR) | ~length(LP_TOL_ILU) | ~length(LP_INFO_QMR)
  error('This routine needs global data which must be generated by calling ''au_qmr_ilu_m_i'' first.');
end 

[n,m] = size(X);

if LP_MC=='C'

  eval(lp_e( 'global LP_L',i,' LP_U',i ));

  is_init1 = eval(lp_e( 'length(LP_L',i,')' ));
  is_init2 = eval(lp_e( 'length(LP_U',i,')' ));
  if ~is_init1 | ~is_init2
    error('This routine needs global data which must be generated by calling ''au_qmr_ilu_s_i'' first.');
  end 

else

  eval(lp_e( '[LP_L',i,',LP_U',i,'] = luinc(LP_A+LP_P(i)*speye(n),LP_TOL_ILU);' ));
  
%
%  This shows the amount of fill-in created by the ILU preconditioner:
%
  if LP_INFO_QMR >= 5
    nnA = nnz(LP_A);
    eval(lp_e( 'nnLU = nnz(LP_L',i,')+nnz(LP_U',i,');'  ));   
    disp(sprintf('ILU-Ratio (nnz(L)+nnz(U))/nnz(A) = %6g',nnLU/nnA)); 
  end

end

Y = zeros(n,m);

if tr=='N'
  for j = 1:m
    eval(lp_e(  '[Y(:,j),fl,rr] = qmr(LP_A+LP_P(i)*speye(n),X(:,j),LP_TOL_QMR,LP_MAX_IT_QMR,LP_L',i,',LP_U',i,');'  ));
    if LP_INFO_QMR & fl
      disp('WARNING in ''au_qmr_ilu_s'': QMR terminated not properly!');
      disp(sprintf('  ( normalized precond. QMR residual = %6g )',rr));  
      rr0 = norm(LP_A*Y(:,j)+LP_P(i)*Y(:,j)-X(:,j))/norm(X(:,j));  
      disp(sprintf('  ( normalized UNprecond. QMR residual = %6g )',rr0));  
    end 
  end
elseif tr=='T'
  for j = 1:m
    eval(lp_e(  '[Y(:,j),fl,rr] = qmr((LP_A.''+LP_P(i)*speye(n)),X(:,j),LP_TOL_QMR,LP_MAX_IT_QMR,LP_U',i,'.'',LP_L',i,'.'');'  ));
    if LP_INFO_QMR & fl
      disp('WARNING in ''au_qmr_ilu_s'': QMR terminated not properly!');
      disp(sprintf('  ( normalized QMR residual = %6g )',rr));  
      rr0 = norm(LP_A.'*Y(:,j)+LP_P(i)*Y(:,j)-X(:,j))/norm(X(:,j));  
      disp(sprintf('  ( normalized UNprecond. QMR residual = %6g )',rr0));  
    end 
  end
else
  error('tp must be either ''N'' or ''T''.');
end