tminres
tminres.hpp
Go to the documentation of this file.
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 
 All Classes Files Functions Friends