tminres
|
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