No BSD License  

Highlights from
LYAPACK

from LYAPACK by Volker Mehrmann
LYAPACK toolbox provides solutions for certain large scale problems related to Lyapunov equations.

au_qmr_ilu_l(tr,X)
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

Contact us at files@mathworks.com