NAME

       minres - minimum residual algorithm (rheolef-7.1)

SYNOPSIS

       template <class Matrix, class Vector, class Preconditioner>
       int minres (const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
         const solver_option& sopt = solver_option())

DESCRIPTION

       This function solves the symmetric positive but possibly singular linear system A*x=b with the minimal
       residual method. The minres function follows the algorithm described in

           C. C. Paige and M. A. Saunders,
           Solution of sparse indefinite systems of linear equations",
           SIAM J. Numer. Anal., 12(4), 1975.

        For more, see http://www.stanford.edu/group/SOL/software.html and also at page 60 of the PhD report:

          S.-C. T. Choi,
          Iterative methods for singular linear equations and least-squares problems,
          Stanford University, 2006,
          http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf

EXAMPLE

           solver_option sopt;
           sopt.max_iter = 100;
           sopt.tol = 1e-7;
           int status = minres (A, x, b, eye(), sopt);

        The fourth argument of minres is a preconditionner: here, the eye(5) one is a do-nothing
       preconditionner, for simplicity. Finally, the solver_option(4) variable sopt transmits the stopping
       criterion with sopt.tol and sopt.max_iter.

       On return, the sopt.residue and sopt.n_iter indicate the reached residue and the number of iterations
       effectively performed. The return status is zero when the prescribed tolerance tol has been obtained, and
       non-zero otherwise. Also, the x variable contains the approximate solution. See also the solver_option(4)
       for more controls upon the stopping criterion.

IMPLEMENTATION

       This documentation has been generated from file linalg/lib/minres.h

       The present template implementation is inspired from the IML++ 1.2 iterative method library,
       http://math.nist.gov/iml++

       template <class Matrix, class Vector, class Preconditioner>
       int minres (const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
         const solver_option& sopt = solver_option())

       {
         // Size &max_iter, Real &tol, odiststream *p_derr = 0
         typedef typename Vector::size_type  Size;
         typedef typename Vector::float_type Real;
         std::string label = (sopt.label != "" ? sopt.label : "minres");
         Vector b = M.solve(Mb);
         Real norm_b = sqrt(fabs(dot(Mb,b)));
         if (sopt.absolute_stopping || 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));
         sopt.residue = norm_r/norm_b;
         if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl
                                       << "[" << label << "] 0 " << sopt.residue << std::endl;
         if (beta2 < 0 || sopt.residue <= sopt.tol) {
           dis_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.ownership(), 0.);
         Vector Mv_old (x.ownership(), 0.);
         Vector w      (x.ownership(), 0.);
         Vector w_old  (x.ownership(), 0.);
         Vector w_old2 (x.ownership(), 0.);
         for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
           // 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) {
               dis_warning_macro ("minres: machine precision problem");
               sopt.residue = norm_r/norm_b;
               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;
           sopt.residue = norm_r/norm_b;
           if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
           if (sopt.residue <= sopt.tol) return 0;
         }
         return 1;
       }

AUTHOR

       Pierre  Saramito  <Pierre.Saramito@imag.fr>

COPYRIGHT

       Copyright   (C)  2000-2018  Pierre  Saramito  <Pierre.Saramito@imag.fr> GPLv3+: GNU GPL version 3 or
       later  <http://gnu.org/licenses/gpl.html>.  This  is  free  software:  you  are free to change and
       redistribute it.  There is NO WARRANTY, to the extent permitted by law.

Version 7.1                                      Sun Oct 31 2021                                minres(5rheolef)