tminres
|
00001 // tminres is free software; you can redistribute it and/or modify it under the 00002 // terms of the GNU Lesser General Public License (as published by the Free 00003 // Software Foundation) version 2.1 dated February 1999. 00004 // 00005 // Authors: 00006 // - Umberto Villa, Emory University - uvilla@emory.edu 00007 // - Michael Saunders, Stanford University 00008 // - Santiago Akle, Stanford University 00009 00095 #include <limits> 00096 #include <memory> 00097 #include <iostream> 00098 #include <iomanip> 00099 #include <vector> 00100 #include <cmath> 00101 #include <algorithm> 00102 #include <limits> 00103 00104 template< typename Operator, typename Vector, typename Preconditioner> 00105 int 00106 MINRES(const Operator &A, Vector &x, const Vector &b, 00107 const Preconditioner * M, double shift, int &max_iter, double &tol, bool show) 00108 { 00109 00110 typedef Operator operator_Type; 00111 typedef Vector vector_Type; 00112 00113 double eps(std::numeric_limits<double>::epsilon()); 00114 00115 std::vector<std::string> msg(11); 00116 msg[0] = " beta1 = 0. The exact solution is x = 0 "; 00117 msg[1] = " A solution to Ax = b was found, given tol "; 00118 msg[2] = " A least-squares solution was found, given tol "; 00119 msg[3] = " Reasonable accuracy achieved, given eps "; 00120 msg[4] = " x has converged to an eigenvector "; 00121 msg[5] = " acond has exceeded 0.1/eps "; 00122 msg[6] = " The iteration limit was reached "; 00123 msg[7] = " A does not define a symmetric matrix "; 00124 msg[8] = " M does not define a symmetric matrix "; 00125 msg[9] = " M does not define a pos-def preconditioner "; 00126 msg[10] = " beta2 = 0. If M = I, b and x are eigenvectors "; 00127 00128 if(show) 00129 { 00130 std::cout<<std::setfill('-')<<std::setw(80)<<"-"<< "\n"; 00131 std::cout<<"| Adapted from minres.m, Stanford University, 06 Jul 2009 |\n"; 00132 std::cout<<"| Solution of symmetric Ax=b or (A-shift*I)x = b |\n"; 00133 std::cout<<std::setfill('-')<<std::setw(80)<<"-"<<"\n"; 00134 std::cout<<std::setfill(' '); 00135 std::cout<<"shift = "<< shift << "; tol = " << tol << "; max iter = " << max_iter<<"\n"; 00136 } 00137 00138 int istop(0), itn(0); 00139 double Anorm(0.0), Acond(0.0), Arnorm(0.0); 00140 double rnorm(0.0), ynorm(0.0); 00141 bool done(false); 00142 00143 // Step 1 00144 /* 00145 * Set up y and v for the first Lanczos vector v1. 00146 * y = beta1 P' v1, where P = C^(-1). 00147 * v is really P'v1 00148 */ 00149 00150 const std::auto_ptr<Vector> r1(x.Clone()); 00151 const std::auto_ptr<Vector> y(x.Clone()); 00152 00153 A.Apply(x, *r1); 00154 add(*r1, -shift, x, *r1); 00155 subtract(b, *r1, *r1); 00156 00157 if(M) 00158 M->Apply(*r1, *y); 00159 else 00160 *y = *r1; 00161 00162 double beta1(0.0); 00163 beta1 = InnerProduct(*r1, *y); 00164 00165 // Test for an indefined preconditioner 00166 // If b = 0 exactly stop with x = x0. 00167 00168 if(beta1 < 0.0) 00169 { 00170 istop = 9; 00171 show = true; 00172 done = true; 00173 } 00174 else 00175 { 00176 if(beta1 == 0.0) 00177 { 00178 show = true; 00179 done = true; 00180 } 00181 else 00182 beta1 = std::sqrt(beta1); // Normalize y to get v1 later 00183 } 00184 00185 // TODO: port symmetry checks for A and M 00186 00187 // STEP 2 00188 /* Initialize other quantities */ 00189 double oldb(0.0), beta(beta1), dbar(0.0), epsln(0.0), oldeps(0.0); 00190 double qrnorm(beta1), phi(0.0), phibar(beta1), rhs1(beta1); 00191 double rhs2(0.0), tnorm2(0.0), ynorm2(0.0); 00192 double cs(-1.0), sn(0.0); 00193 double gmax(0.0), gmin(std::numeric_limits<double>::max( )); 00194 double alpha(0.0), gamma(0.0); 00195 double delta(0.0), gbar(0.0); 00196 double z(0.0); 00197 00198 const std::auto_ptr<Vector> w(x.Clone()), w1(x.Clone()), w2(x.Clone()), r2(x.Clone()); 00199 ( *w ) = (* w1 ) = ( *w2 ) = 0.0; 00200 *r2 = *r1; 00201 const std::auto_ptr<Vector> v(x.Clone()); 00202 00203 if(show) 00204 std::cout<<std::setw(6)<<"Itn" 00205 << std::setw(14) << "Compatible" 00206 << std::setw(14) << "LS" 00207 << std::setw(14) << "norm(A)" 00208 << std::setw(14) << "cond(A)" 00209 << std::setw(14) << "gbar/|A|"<<"\n"; 00210 00211 /* Main Iteration */ 00212 if(!done) 00213 { 00214 for(itn = 0; itn < max_iter; ++itn) 00215 { 00216 // STEP 3 00217 /* 00218 ----------------------------------------------------------------- 00219 Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,... 00220 The general iteration is similar to the case k = 1 with v0 = 0: 00221 00222 p1 = Operator * v1 - beta1 * v0, 00223 alpha1 = v1'p1, 00224 q2 = p2 - alpha1 * v1, 00225 beta2^2 = q2'q2, 00226 v2 = (1/beta2) q2. 00227 00228 Again, y = betak P vk, where P = C**(-1). 00229 .... more description needed. 00230 ----------------------------------------------------------------- 00231 */ 00232 double s(1./beta); //Normalize previous vector (in y) 00233 (*v) = (*y); 00234 v->Scale(s); // v = vk if P = I 00235 00236 A.Apply(*v,*y); 00237 add(*y, -shift, *v, *y); 00238 if(itn) 00239 add(*y, -beta/oldb, *r1, *y); 00240 00241 alpha = InnerProduct(*v, *y); // alphak 00242 add(*y, -alpha/beta, *r2, *y); 00243 (*r1) = (*r2); 00244 (*r2) = (*y); 00245 if(M) 00246 M->Apply(*r2, *y); 00247 else 00248 *y = *r2; 00249 00250 oldb = beta; //oldb = betak 00251 beta = InnerProduct(*r2, *y); 00252 00253 if(beta < 0) 00254 { 00255 istop = 9; 00256 break; 00257 } 00258 00259 beta = sqrt(beta); 00260 tnorm2 += alpha*alpha + oldb*oldb + beta*beta; 00261 00262 if(itn == 0) //Initialize a few things 00263 { 00264 if(beta/beta1 < 10.0*eps) 00265 istop = 10; 00266 } 00267 00268 // Apply previous rotation Q_{k-1} to get 00269 // [delta_k epsln_{k+1}] = [cs sn] [dbar_k 0] 00270 // [gbar_k dbar_{k+1}] [sn -cs] [alpha_k beta_{k+1}]. 00271 oldeps = epsln; 00272 delta = cs*dbar + sn*alpha; 00273 gbar = sn*dbar - cs*alpha; 00274 epsln = sn*beta; 00275 dbar = - cs*beta; 00276 double root(sqrt(gbar*gbar + dbar*dbar)); 00277 Arnorm = phibar * root; // ||Ar_{k-1}|| 00278 00279 // Compute next plane rotation Q_k 00280 gamma = sqrt(gbar*gbar + beta*beta); // gamma_k 00281 gamma = std::max(gamma, eps); 00282 cs = gbar/gamma; // c_k 00283 sn = beta/gamma; // s_k 00284 phi = cs*phibar; // phi_k 00285 phibar = sn*phibar; // phibar_{k+1} 00286 00287 00288 // Update x 00289 double denom(1./gamma); 00290 (*w1) = (*w2); 00291 (*w2) = (*w); 00292 add(-oldeps, *w1, -delta, *w2, *w); 00293 add(denom, *v, *w, *w); 00294 add(x, phi, *w, x); 00295 00296 // go round again 00297 gmax = std::max(gmax, gamma); 00298 gmin = std::min(gmin, gamma); 00299 z = rhs1/gamma; 00300 rhs1 = rhs2 - delta*z; 00301 rhs2 = - epsln*z; 00302 00303 // Estimate various norms 00304 00305 Anorm = sqrt(tnorm2); 00306 ynorm2 = InnerProduct(x, x); 00307 ynorm = sqrt(ynorm2); 00308 double epsa(Anorm*eps); 00309 double epsx(epsa*ynorm); 00310 double epsr(Anorm*ynorm*tol); 00311 double diag(gbar); 00312 if(0 == diag) 00313 diag = epsa; 00314 00315 qrnorm = phibar; 00316 rnorm = qrnorm; 00317 double test1(0.0), test2(0.0); 00318 test1 = rnorm/(Anorm*ynorm); // ||r||/(||A|| ||x||) 00319 test2 = root/ Anorm; // ||A r_{k-1}|| / (||A|| ||r_{k-1}||) 00320 00321 // Estimate cond(A) 00322 /* 00323 In this version we look at the diagonals of R in the 00324 factorization of the lower Hessenberg matrix, Q * H = R, 00325 where H is the tridiagonal matrix from Lanczos with one 00326 extra row, beta(k+1) e_k^T. 00327 */ 00328 Acond = gmax/gmin; 00329 00330 //See if any of the stopping criteria is satisfied 00331 if(0 == istop) 00332 { 00333 double t1(1.0+test1), t2(1.0+test2); //This test work if tol < eps 00334 if(t2 <= 1.) istop = 2; 00335 if(t1 <= 1.) istop = 1; 00336 00337 if(itn >= max_iter-1) istop = 6; 00338 if(Acond >= .1/eps) istop = 4; 00339 if(epsx >= beta1) istop = 3; 00340 if(test2 <= tol ) istop = 2; 00341 if(test1 <= tol) istop = 1; 00342 } 00343 00344 if(show) 00345 std::cout<< std::setw(6) << itn 00346 << std::setw(14) << test1 00347 << std::setw(14) << test2 00348 << std::setw(14) << Anorm 00349 << std::setw(14) << Acond 00350 << std::setw(14) << gbar/Anorm << std::endl; 00351 00352 if(0 != istop) 00353 break; 00354 00355 } 00356 } 00357 00358 // Display final status 00359 if(show) 00360 { 00361 std::cout << std::setfill('-') << std::setw(80) << "-" << "\n"; 00362 std::cout << msg[istop] << "\n"; 00363 std::cout << " Number of iterations: " << itn << "\n"; 00364 std::cout << " Anorm = " << Anorm << "\t Acond = " << Acond << "\n"; 00365 std::cout << " rnorm = " << rnorm << "\t ynorm = " << ynorm << "\n"; 00366 std::cout << " Arnorm = " << Arnorm << "\n"; 00367 std::cout << std::setfill('-') << std::setw(80) << "-" << std::endl; 00368 std::cout << std::setfill(' '); 00369 } 00370 00371 return istop; 00372 } 00373