librheolef-dev 6.5-1build1 / usr / include / rheolef / pminres.h

# ifndef _SKIT_PMINRES_H
# define _SKIT_PMINRES_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
#include "rheolef/diststream.h"

namespace rheolef { 
/*D:pminres
NAME: @code{pminres} -- conjugate gradient algorithm.
@findex pminres
@cindex conjugate gradient algorithm
@cindex iterative solver
@cindex preconditioner
SYNOPSIS:
  @example
    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, odiststream *p_derr=0);
  @end example
EXAMPLE:
  @noindent
  The simplest call to 'pminres' has the folling form:
  @example
    size_t max_iter = 100;
    double tol = 1e-7;
    int status = pminres(a, x, b, EYE, max_iter, tol, &derr);
  @end example
DESCRIPTION:       
  @noindent
  @code{pminres} solves the symmetric positive definite linear
  system Ax=b using the Conjugate Gradient method.

  @noindent
  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:
  @table @code
    @item x
	approximate solution to Ax = b

    @item max_iter
	the number of iterations performed before the tolerance was reached

    @item tol
	the residual after the final iteration
  @end table
NOTE: 
  @noindent
  @code{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.
  @noindent
  The present implementation style is inspired from @code{IML++ 1.2} iterative method library,
  @url{http://math.nist.gov/iml++}.
AUTHOR: 
    Pierre Saramito
    | Pierre.Saramito@imag.fr
    LJK-IMAG, 38041 Grenoble cedex 9, France
DATE: 
    22 april 2009
METHODS: @pminres
End:
*/
//<pminres:
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, odiststream *p_derr = 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_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
  if (p_derr) (*p_derr) << "[" << label << "] 0 " << norm_r/norm_b << std::endl;
  if (beta2 < 0 || norm_r <= tol*norm_b) {
    tol = norm_r/norm_b;
    max_iter = 0;
    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 (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) { 
	dis_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_derr) (*p_derr) << "[" << 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;
}
//>pminres:
}// namespace rheolef
# endif // _SKIT_PMINRES_H