tminres
minres.m
00001 function [ x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm ] = ...
00002            minres( A, b, M, shift, show, check, itnlim, rtol )
00003 
00004 %        [ x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm ] = ...
00005 %          minres( A, b, M, shift, show, check, itnlim, rtol )
00006 %
00007 % minres solves the n x n system of linear equations Ax = b
00008 % or the n x n least squares problem           min ||Ax - b||_2^2,
00009 % where A is a symmetric matrix (possibly indefinite or singular)
00010 % and b is a given vector.  The dimension n is defined by length(b).
00011 %
00012 % INPUT:
00013 %
00014 % "A" may be a dense or sparse matrix (preferably sparse!)
00015 % or a function handle such that y = A(x) returns the product
00016 % y = A*x for any given n-vector x.
00017 %
00018 % If M = [], preconditioning is not used.  Otherwise,
00019 % "M" defines a positive-definite preconditioner M = C*C'.
00020 % "M" may be a dense or sparse matrix (preferably sparse!)
00021 % or a function handle such that y = M(x) solves the system
00022 % My = x given any n-vector x.
00023 %
00024 % If shift ~= 0, minres really solves (A - shift*I)x = b
00025 % (or the corresponding least-squares problem if shift is an
00026 % eigenvalue of A).
00027 %
00028 % When M = C*C' exists, minres implicitly solves the system
00029 %
00030 %            P(A - shift*I)P'xbar = Pb,
00031 %    i.e.               Abar xbar = bbar,
00032 %    where                      P = inv(C),
00033 %                            Abar = P(A - shift*I)P',
00034 %                            bbar = Pb,
00035 %
00036 % and returns the solution      x = P'xbar.
00037 % The associated residual is rbar = bbar - Abar xbar
00038 %                                 = P(b - (A - shift*I)x)
00039 %                                 = Pr.
00040 %
00041 % OUTPUT:
00042 %
00043 % x      is the final estimate of the required solution
00044 %        after k iterations, where k is return in itn.
00045 % istop  is a value from [-1:9] to indicate the reason for termination.
00046 %        The reason is summarized in msg[istop+2] below.
00047 % itn    gives the final value of k (the iteration number).
00048 % rnorm  estimates norm(r_k)  or norm(rbar_k) if M exists.
00049 % Arnorm estimates norm(Ar_{k-1}) or norm(Abar rbar_{k-1}) if M exists.
00050 %        NOTE THAT Arnorm LAGS AN ITERATION BEHIND rnorm.
00051 
00052 % Code authors:Michael Saunders, SOL, Stanford University
00053 %              Sou Cheng Choi,  SCCM, Stanford University
00054 %
00055 % 02 Sep 2003: Date of Fortran 77 version, based on 
00056 %              C. C. Paige and M. A. Saunders (1975),
00057 %              Solution of sparse indefinite systems of linear equations,
00058 %              SIAM J. Numer. Anal. 12(4), pp. 617-629.
00059 %
00060 % 02 Sep 2003: ||Ar|| now estimated as Arnorm.
00061 % 17 Oct 2003: f77 version converted to MATLAB.
00062 % 03 Apr 2005: A must be a matrix or a function handle.
00063 % 10 May 2009: Parameter list shortened.
00064 %              Documentation updated following suggestions from
00065 %              Jeffery Kline <jeffery.kline@gmail.com>
00066 %              (author of new Python versions of minres, symmlq, lsqr).
00067 % 06 Jul 2009: Michael Chen <mc462@cornell.edu> reports divide by zero
00068 %              when beta = 0 (in this case it was beta_2 = 0).
00069 %              Realized that the istop values were out of sync.
00070 %              They should be right now.
00071 % 02 Sep 2011: David Fong reports error in Acond when alpha1=0.
00072 %              gmax and gmin should be initialized before itn 1.
00073 % 02 Sep 2011: ynorm = norm(x) is now computed directly instead of
00074 %              being updated (incorrectly).
00075 
00076 % Known bugs:
00077 %  1. As Jeff Kline pointed out, Arnorm = ||A r_{k-1}|| lags behind
00078 %     rnorm = ||r_k||.  On singular systems, this means that a good
00079 %     least-squares solution exists before Arnorm is small enough
00080 %     to recognize it.  The solution x_{k-1} gets updated to x_k
00081 %     (possibly a very large solution) before Arnorm shuts things
00082 %     down the next iteration.  It would be better to keep x_{k-1}.
00083 %------------------------------------------------------------------
00084 
00085 %  Initialize                               
00086 
00087 msg = [' beta2 = 0.  If M = I, b and x are eigenvectors '   % -1
00088        ' beta1 = 0.  The exact solution is  x = 0       '   %  0
00089        ' A solution to Ax = b was found, given rtol     '   %  1
00090        ' A least-squares solution was found, given rtol '   %  2
00091        ' Reasonable accuracy achieved, given eps        '   %  3
00092        ' x has converged to an eigenvector              '   %  4
00093        ' acond has exceeded 0.1/eps                     '   %  5
00094        ' The iteration limit was reached                '   %  6
00095        ' A  does not define a symmetric matrix          '   %  7
00096        ' M  does not define a symmetric matrix          '   %  8
00097        ' M  does not define a pos-def preconditioner    ']; %  9
00098  
00099 n      = length(b);
00100 precon = ~isempty(M);
00101 if show
00102    fprintf('\n minres.m   SOL, Stanford University   Version of 06 Jul 2009')
00103    fprintf('\n Solution of symmetric Ax = b or (A-shift*I)x = b')
00104    fprintf('\n\n n      =%8g    shift =%22.14e', n,shift)
00105    fprintf('\n itnlim =%8g    rtol  =%10.2e\n', itnlim,rtol)
00106 end
00107 
00108 istop = 0;   itn   = 0;   Anorm = 0;    Acond = 0;
00109 rnorm = 0;   ynorm = 0;   done  = false;
00110 x     = zeros(n,1);
00111 
00112 %---------------------------------------------------------------------
00113 % Decode A.
00114 %---------------------------------------------------------------------
00115 if isa(A,'double')         % A is an explicit matrix A.
00116   if issparse(A)
00117     nnzA = nnz(A);
00118     fprintf('\n A is an explicit sparse matrix')
00119     fprintf('\n nnz(A) =%8g', nnzA)
00120   else
00121     fprintf('\n A is an explicit dense matrix' )
00122   end
00123 elseif isa(A,'function_handle')
00124   disp(['The matrix A is defined by function_handle ' func2str(A)])
00125 else
00126   error('minres','A must be a matrix or a function handle')
00127 end
00128 
00129 %---------------------------------------------------------------------
00130 % Decode M.
00131 %---------------------------------------------------------------------
00132 if precon
00133   if isa(M,'double')       % M is an explicit matrix M.
00134     if issparse(M)
00135       nnzM = nnz(M);
00136       fprintf('\n The matrix M is an explicit sparse matrix')
00137       fprintf('\n nnz(M) =%8g', nnzM)
00138     else
00139       fprintf('\n The matrix M is an explicit dense matrix' )
00140     end
00141   elseif isa(M,'function_handle')
00142     disp(['The matrix M is defined by function_handle ' func2str(M)])
00143   else
00144     error('minres','M must be a matrix or a function handle')
00145   end
00146 end
00147 
00148 %------------------------------------------------------------------
00149 % Set up y and v for the first Lanczos vector v1.
00150 % y  =  beta1 P' v1,  where  P = C**(-1).
00151 % v is really P' v1.
00152 %------------------------------------------------------------------
00153 y     = b;
00154 r1    = b;
00155 if precon, y = minresxxxM( M,b ); end
00156 beta1 = b'*y;
00157                                                 
00158 %  Test for an indefinite preconditioner.
00159 %  If b = 0 exactly, stop with x = 0.
00160 
00161 if beta1< 0, istop = 9;  show = true;  done = true; end
00162 if beta1==0,             show = true;  done = true; end
00163 
00164 if beta1> 0
00165   beta1  = sqrt(beta1);       % Normalize y to get v1 later.
00166 
00167   % See if M is symmetric.
00168 
00169   if check & precon
00170     r2     = minresxxxM( M,y );
00171     s      = y' *y;
00172     t      = r1'*r2;
00173     z      = abs(s-t);
00174     epsa   = (s+eps)*eps^(1/3);
00175     if z > epsa, istop = 8;  show = true;  done = true; end
00176   end
00177 
00178   % See if A is symmetric.
00179 
00180   if check
00181      w    = minresxxxA( A,y );
00182      r2   = minresxxxA( A,w );
00183      s    = w'*w;
00184      t    = y'*r2;
00185      z    = abs(s-t);
00186      epsa = (s+eps)*eps^(1/3);
00187      if z > epsa, istop = 7;  done  = true;  show = true;  end
00188   end
00189 end
00190 
00191 %------------------------------------------------------------------
00192 % Initialize other quantities.
00193 % ------------------------------------------------------------------
00194 oldb   = 0;       beta   = beta1;   dbar   = 0;       epsln  = 0;
00195 qrnorm = beta1;   phibar = beta1;   rhs1   = beta1;
00196 rhs2   = 0;       tnorm2 = 0;       gmax   = 0;       gmin   = realmax;
00197 cs     = -1;      sn     = 0;
00198 w      = zeros(n,1);
00199 w2     = zeros(n,1);
00200 r2     = r1;
00201 
00202 if show
00203   fprintf('\n\n   Itn     x(1)     Compatible    LS       norm(A)  cond(A)')
00204   fprintf(' gbar/|A|\n')   %%%%%% Check gbar
00205 end
00206 
00207 %---------------------------------------------------------------------
00208 % Main iteration loop.
00209 % --------------------------------------------------------------------
00210 if ~done                              % k = itn = 1 first time through
00211   while itn < itnlim
00212     itn    = itn+1;
00213 
00214     %-----------------------------------------------------------------
00215     % Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
00216     % The general iteration is similar to the case k = 1 with v0 = 0:
00217     %
00218     %   p1      = Operator * v1  -  beta1 * v0,
00219     %   alpha1  = v1'p1,
00220     %   q2      = p2  -  alpha1 * v1,
00221     %   beta2^2 = q2'q2,
00222     %   v2      = (1/beta2) q2.
00223     %
00224     % Again, y = betak P vk,  where  P = C**(-1).
00225     % .... more description needed.
00226     %-----------------------------------------------------------------
00227     s = 1/beta;                 % Normalize previous vector (in y).
00228     v = s*y;                    % v = vk if P = I
00229 
00230     y = minresxxxA( A,v ) - shift*v;
00231     if itn >= 2, y = y - (beta/oldb)*r1; end
00232 
00233     alfa   = v'*y;              % alphak
00234     y      = (- alfa/beta)*r2 + y;
00235     r1     = r2;
00236     r2     = y;
00237     if precon,  y = minresxxxM( M,r2 );  end
00238     oldb   = beta;              % oldb = betak
00239     beta   = r2'*y;             % beta = betak+1^2
00240     if beta < 0, istop = 9;  break;  end
00241     beta   = sqrt(beta);
00242     tnorm2 = tnorm2 + alfa^2 + oldb^2 + beta^2;
00243 
00244     if itn==1                   % Initialize a few things.
00245       if beta/beta1 <= 10*eps   % beta2 = 0 or ~ 0.
00246         istop = -1;             % Terminate later.
00247       end
00248     end
00249 
00250     % Apply previous rotation Qk-1 to get
00251     %   [deltak epslnk+1] = [cs  sn][dbark    0   ]
00252     %   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
00253 
00254     oldeps = epsln;
00255     delta  = cs*dbar + sn*alfa; % delta1 = 0         deltak
00256     gbar   = sn*dbar - cs*alfa; % gbar 1 = alfa1     gbar k
00257     epsln  =           sn*beta; % epsln2 = 0         epslnk+1
00258     dbar   =         - cs*beta; % dbar 2 = beta2     dbar k+1
00259     root   = norm([gbar dbar]);
00260     Arnorm = phibar*root;       % ||Ar{k-1}||
00261 
00262     % Compute the next plane rotation Qk
00263 
00264     gamma  = norm([gbar beta]); % gammak
00265     gamma  = max([gamma eps]);
00266     cs     = gbar/gamma;        % ck
00267     sn     = beta/gamma;        % sk
00268     phi    = cs*phibar ;        % phik
00269     phibar = sn*phibar ;        % phibark+1
00270 
00271     % Update  x.
00272 
00273     denom = 1/gamma;
00274     w1    = w2;
00275     w2    = w;
00276     w     = (v - oldeps*w1 - delta*w2)*denom;
00277     x     = x + phi*w;
00278 
00279     % Go round again.
00280 
00281     gmax   = max([gmax gamma]);
00282     gmin   = min([gmin gamma]);
00283     z      = rhs1/gamma;
00284     rhs1   = rhs2 - delta*z;
00285     rhs2   =      - epsln*z;
00286 
00287     % Estimate various norms.
00288 
00289     Anorm  = sqrt( tnorm2 );
00290     ynorm  = norm(x);
00291     epsa   = Anorm*eps;
00292     epsx   = Anorm*ynorm*eps;
00293     epsr   = Anorm*ynorm*rtol;
00294     diag   = gbar;
00295     if diag==0, diag = epsa; end
00296 
00297     qrnorm = phibar;
00298     rnorm  = qrnorm;
00299     test1  = rnorm/(Anorm*ynorm);    %  ||r|| / (||A|| ||x||)
00300     test2  = root / Anorm;      % ||Ar{k-1}|| / (||A|| ||r_{k-1}||)
00301 
00302     % Estimate  cond(A).
00303     % In this version we look at the diagonals of  R  in the
00304     % factorization of the lower Hessenberg matrix,  Q * H = R,
00305     % where H is the tridiagonal matrix from Lanczos with one
00306     % extra row, beta(k+1) e_k^T.
00307 
00308     Acond  = gmax/gmin;
00309 
00310     % See if any of the stopping criteria are satisfied.
00311     % In rare cases, istop is already -1 from above (Abar = const*I).
00312 
00313     if istop==0
00314       t1 = 1 + test1;      % These tests work if rtol < eps
00315       t2 = 1 + test2;
00316       if t2    <= 1      , istop = 2; end
00317       if t1    <= 1      , istop = 1; end
00318       
00319       if itn   >= itnlim , istop = 6; end
00320       if Acond >= 0.1/eps, istop = 4; end
00321       if epsx  >= beta1  , istop = 3; end
00322      %if rnorm <= epsx   , istop = 2; end
00323      %if rnorm <= epsr   , istop = 1; end
00324       if test2 <= rtol   , istop = 2; end
00325       if test1 <= rtol   , istop = 1; end
00326     end
00327 
00328     % See if it is time to print something.
00329 
00330     prnt   = false;
00331     if n      <= 40       , prnt = true; end
00332     if itn    <= 10       , prnt = true; end
00333     if itn    >= itnlim-10, prnt = true; end
00334     if mod(itn,10)==0     , prnt = true; end
00335     if qrnorm <= 10*epsx  , prnt = true; end
00336     if qrnorm <= 10*epsr  , prnt = true; end
00337     if Acond  <= 1e-2/eps , prnt = true; end
00338     if istop  ~=  0       , prnt = true; end
00339 
00340     if show & prnt
00341       if mod(itn,10)==0, disp(' '); end
00342       str1 = sprintf('%6g %12.5e %10.3e', itn,x(1),test1);
00343       str2 = sprintf(' %10.3e',           test2);
00344       str3 = sprintf(' %8.1e %8.1e',      Anorm,Acond);
00345       str4 = sprintf(' %8.1e',            gbar/Anorm);
00346       str  = [str1 str2 str3 str4];
00347       fprintf('\n %s', str)
00348 
00349       debug = false;  % true;
00350       if debug   % Print true Arnorm.
00351                  % This works only if no preconditioning.
00352         vv = b - minresxxxA(A,x)  + shift*x;    % vv = b - (A - shift*I)*x
00353         ww =     minresxxxA(A,vv) - shift*vv;   % ww = (A - shift*I)*vv = "Ar"
00354         trueArnorm = norm(ww);
00355         fprintf('\n Arnorm = %12.4e   True ||Ar|| = %12.4e', Arnorm,trueArnorm)
00356       end
00357     end % show & prnt
00358 
00359     if istop ~= 0, break; end
00360 
00361   end % main loop
00362 end % if ~done early
00363 
00364 % Display final status.
00365 
00366 if show
00367   fprintf('\n\n istop   =%3g               itn   =%6g', istop,itn  )
00368   fprintf('\n Anorm   =%12.4e      Acond =%12.4e', Anorm,Acond)
00369   fprintf('\n rnorm   =%12.4e      ynorm =%12.4e', rnorm,ynorm)
00370   fprintf('\n Arnorm  =%12.4e\n', Arnorm)
00371   disp(msg(istop+2,:))
00372 end
00373 %-----------------------------------------------------------------------
00374 % End function pdco.m
00375 %-----------------------------------------------------------------------
00376 
00377 
00378 function y = minresxxxA( A,x )
00379 
00380 % sets   y = A*x for a matrix A defined by input parameter A.
00381 %
00382 % 10 May 2009: A is now a matrix or a function handle.
00383 
00384   if isa(A,'double')
00385     y = A*x;
00386   else
00387     y = A(x);
00388   end
00389 %-----------------------------------------------------------------------
00390 % End private function minresxxxA
00391 %-----------------------------------------------------------------------
00392 
00393 
00394 function y = minresxxxM( M,x )
00395 
00396 % solves My = x, where M is the preconditioner.
00397 %
00398 % 10 May 2009: M is now a matrix or a function handle.
00399 
00400   if isa(M,'double')
00401     y = M\x;
00402   else
00403     y = M(x);
00404   end
00405 %-----------------------------------------------------------------------
00406 % End private function minresxxxM
00407 %-----------------------------------------------------------------------
00408 
 All Classes Files Functions Friends