No BSD License  

Highlights from
Templates for the Solution of Linear Systems

image thumbnail

Templates for the Solution of Linear Systems

by

 

19 Aug 2002 (Updated )

Companion Software

qmr( A, x, b, M, max_it, tol )
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

Contact us