function Y = au_qmr_ilu_l(tr,X)
%
% Solves linear systems with the real matrix A or its transpose A' using
% the QMR iteration and the ILU preconditioner.
%
% for tr = 'N':
%
% Y = inv(A)*X,
%
% for tr = 'T':
%
% Y = inv(A')*X.
%
% The LU factors of the ILU preconditioner for A are provided as global
% data. This data must be generated by calling 'au_qmr_ilu_l_i' before
% calling this routine!
%
% Calling sequence:
%
% Y = au_qmr_ilu_l(tr,X)
%
% Input:
%
% tr (= 'N' or 'T') determines whether systems with A or A'
% should be solved;
% X matrix of proper size.
%
% Output:
%
% Y the solution matrix.
%
%
% LYAPACK 1.0 (Thilo Penzl, August 1999)
if nargin~=2
error('Wrong number of input arguments.');
end
global LP_A LP_L LP_U LP_MAX_IT_QMR LP_TOL_QMR LP_TOL_ILU LP_INFO_QMR
if ~length(LP_A) | ~length(LP_L) | ~length(LP_U) | ~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_l_i'' first.');
end
[n,m] = size(X);
Y = zeros(n,m);
if tr=='N'
for i = 1:m
[Y(:,i),fl,rr] = qmr(LP_A,X(:,i),LP_TOL_QMR,LP_MAX_IT_QMR,LP_L,LP_U);
if LP_INFO_QMR & fl
disp('WARNING in ''au_qmr_ilu_l'': QMR terminated not properly!');
disp(sprintf(' ( normalized precond. QMR residual = %6g )',rr));
rr0 = norm(LP_A*Y(:,i)-X(:,i))/norm(X(:,i));
disp(sprintf(' ( normalized UNprecond. QMR residual = %6g )',rr0));
end
end
elseif tr=='T'
for i = 1:m
[Y(:,i),fl,rr] = qmr(LP_A.',X(:,i),LP_TOL_QMR,LP_MAX_IT_QMR,LP_U.',LP_L.');
if LP_INFO_QMR & fl
disp('WARNING in ''au_qmr_ilu_l'': QMR terminated not properly!');
disp(sprintf(' ( normalized precond. QMR residual = %6g )',rr));
rr0 = norm(LP_A.'*Y(:,i)-X(:,i))/norm(X(:,i));
disp(sprintf(' ( normalized UNprecond. QMR residual = %6g )',rr0));
end
end
else
error('tp must be either ''N'' or ''T''.');
end