Provided by: librheolef-dev_5.93-2_i386 bug

NAME

       pminres -- conjugate gradient algorithm.

SYNOPSIS

           template <class Matrix, class Vector, class Preconditioner, class Real>
           int pminres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
             int &max_iter, Real &tol, std::ostream *p_cerr=0);

EXAMPLE

       The simplest call to 'pminres' has the folling form:

           size_t max_iter = 100;
           double tol = 1e-7;
           int status = pminres(a, x, b, EYE, max_iter, tol, &cerr);

DESCRIPTION

       pminres solves the symmetric positive definite linear system Ax=b using
       the Conjugate Gradient method.

       The  return  value  indicates  convergence  within   max_iter   (input)
       iterations (0), or no convergence within max_iter iterations (1).  Upon
       successful return, output arguments have the following values:

       x      approximate solution to Ax = b

       max_iter
              the number of iterations  performed  before  the  tolerance  was
              reached

       tol    the residual after the final iteration

NOTE

       pminres   follows  the  algorithm  described  in  "Solution  of  sparse
       indefinite systems  of  linear  equations",  C.  C.  Paige  and  M.  A.
       Saunders,   SIAM   J.   Numer.  Anal.,  12(4),  1975.   For  more,  see
       http://www.stanford.edu/group/SOL/software.html  and   also   the   PhD
       "Iterative  methods  for  singular  linear  equations and least-squares
       problems",    S.-C.    T.    Choi,    Stanford    University,     2006,
       http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-
       thesis.pdf at page 60.  The present implementation  style  is  inspired
       from IML++ 1.2 iterative method library, http://math.nist.gov/iml++.

IMPLEMENTATION

       template <class Matrix, class Vector, class Preconditioner, class Real, class Size>
       int pminres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
         Size &max_iter, Real &tol, std::ostream *p_cerr = 0, std::string label = "minres")
       {
         Vector b = M.solve(Mb);
         Real norm_b = sqrt(fabs(dot(Mb,b)));
         if (norm_b == Real(0.)) norm_b = 1;

         Vector Mr = Mb - A*x;
         Vector z = M.solve(Mr);
         Real beta2 = dot(Mr, z);
         Real norm_r = sqrt(fabs(beta2));
         if (p_cerr) (*p_cerr) << "[" << label << "] #iteration residue" << std::endl;
         if (p_cerr) (*p_cerr) << "[" << label << "] 0 " << norm_r/norm_b << std::endl;
         if (beta2 < 0 || norm_r <= tol*norm_b) {
           tol = norm_r/norm_b;
           max_iter = 0;
           warning_macro ("beta2 = " << beta2 << " < 0: stop");
           return 0;
         }
         Real beta = sqrt(beta2);
         Real eta = beta;
         Vector Mv = Mr/beta;
         Vector  u =  z/beta;
         Real c_old = 1.;
         Real s_old = 0.;
         Real c = 1.;
         Real s = 0.;
         Vector u_old  (x.size(), 0.);
         Vector Mv_old (x.size(), 0.);
         Vector w      (x.size(), 0.);
         Vector w_old  (x.size(), 0.);
         Vector w_old2 (x.size(), 0.);
         for (Size n = 1; n <= max_iter; n++) {
           // Lanczos
           Mr = A*u;
           z = M.solve(Mr);
           Real alpha = dot(Mr, u);
           Mr = Mr - alpha*Mv - beta*Mv_old;
           z  =  z - alpha*u  - beta*u_old;
           beta2 = dot(Mr, z);
           if (beta2 < 0) {
               warning_macro ("pminres: machine precision problem");
               tol = norm_r/norm_b;
               max_iter = n;
               return 2;
           }
           Real beta_old = beta;
           beta = sqrt(beta2);
           // QR factorisation
           Real c_old2 = c_old;
           Real s_old2 = s_old;
           c_old = c;
           s_old = s;
           Real rho0 = c_old*alpha - c_old2*s_old*beta_old;
           Real rho2 = s_old*alpha + c_old2*c_old*beta_old;
           Real rho1 = sqrt(sqr(rho0) + sqr(beta));
           Real rho3 = s_old2 * beta_old;
           // Givens rotation
           c = rho0 / rho1;
           s = beta / rho1;
           // update
           w_old2 = w_old;
           w_old  = w;
           w = (u - rho2*w_old - rho3*w_old2)/rho1;
           x += c*eta*w;
           eta = -s*eta;
           Mv_old = Mv;
            u_old = u;
           Mv = Mr/beta;
            u =  z/beta;
           // check residue
           norm_r *= s;
           if (p_cerr) (*p_cerr) << "[" << label << "] " << n << " " << norm_r/norm_b << std::endl;
           if (norm_r <= tol*norm_b) {
             tol = norm_r/norm_b;
             max_iter = n;
             return 0;
           }
         }
         tol = norm_r/norm_b;
         return 1;
       }