function [x, error, iter, flag] = qmr( A, x, b, M, max_it, tol )
% -- Iterative template routine --
% Univ. of Tennessee and Oak Ridge National Laboratory
% October 1, 1993
% Details of this algorithm are described in "Templates for the
% Solution of Linear Systems: Building Blocks for Iterative
% Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
% Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
% 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
% [x, error, iter, flag] = qmr( A, x, b, M, max_it, tol )
%
% qmr.m solves the linear system Ax=b using the
% Quasi Minimal Residual Method with preconditioning.
%
% input A REAL matrix
% x REAL initial guess vector
% b REAL right hand side vector
% M REAL preconditioner
% max_it INTEGER maximum number of iterations
% tol REAL error tolerance
%
% output x REAL solution vector
% error REAL error norm
% iter INTEGER number of iterations performed
% flag INTEGER: 0: solution found to tolerance
% 1: no convergence given max_it
% breakdown:
% -1: rho
% -2: beta
% -3: gamma
% -4: delta
% -5: ep
% -6: xi
iter = 0; % initialization
flag = 0;
bnrm2 = norm( b );
if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end
r = b - A*x;
error = norm( r ) / bnrm2;
if ( error < tol ) return, end
[M1,M2] = lu( M );
v_tld = r;
y = M1 \ v_tld;
rho = norm( y );
w_tld = r;
z = M2' \ w_tld;
xi = norm( z );
gamma = 1.0;
eta = -1.0;
theta = 0.0;
for iter = 1:max_it, % begin iteration
if ( rho == 0.0 | xi == 0.0 ), break, end
v = v_tld / rho;
y = y / rho;
w = w_tld / xi;
z = z / xi;
delta = z'*y;
if ( delta == 0.0 ), break, end
y_tld = M2 \ y;
z_tld = M1'\ z;
if ( iter > 1 ), % direction vector
p = y_tld - ( xi*delta / ep )*p;
q = z_tld - ( rho*delta / ep )*q;
else
p = y_tld;
q = z_tld;
end
p_tld = A*p;
ep = q'*p_tld;
if ( ep == 0.0 ), break, end
beta = ep / delta;
if ( beta == 0.0 ), break, end
v_tld = p_tld - beta*v;
y = M1 \ v_tld;
rho_1 = rho;
rho = norm( y );
w_tld = ( A'*q ) - ( beta*w );
z = M2' \ w_tld;
xi = norm( z );
gamma_1 = gamma;
theta_1 = theta;
theta = rho / ( gamma_1*beta );
gamma = 1.0 / sqrt( 1.0 + (theta^2) );
if ( gamma == 0.0 ), break, end
eta = -eta*rho_1*(gamma^2) / ( beta*(gamma_1^2) );
if ( iter > 1 ), % compute adjustment
d = eta*p + (( theta_1*gamma )^2)*d;
s = eta*p_tld + (( theta_1*gamma )^2)*s;
else
d = eta*p;
s = eta*p_tld;
end
x = x + d; % update approximation
r = r - s; % update residual
error = norm( r ) / bnrm2; % check convergence
if ( error <= tol ), break, end
end
if ( error <= tol ), % converged
flag = 0;
elseif ( rho == 0.0 ), % breakdown
flag = -1;
elseif ( beta == 0.0 ),
flag = -2;
elseif ( gamma == 0.0 ),
flag = -3;
elseif ( delta == 0.0 ),
flag = -4;
elseif ( ep == 0.0 ),
flag = -5;
elseif ( xi == 0.0 ),
flag = -6;
else % no convergence
flag = 1;
end
% END qmr.m
.