/usr/include/rheolef/pcg.h

# ifndef _SKIT_PCG_H
# define _SKIT_PCG_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
/// 
/// =========================================================================
/*D:pcg
NAME: @code{pcg} -- conjugate gradient algorithm.
@findex pcg
@cindex conjugate gradient algorithm
@cindex iterative solver
@cindex preconditioner
SYNOPSIS:
  @example
    template <class Matrix, class Vector, class Preconditioner, class Real>
    int pcg (const Matrix &A, Vector &x, const Vector &b,
      const Preconditioner &M, int &max_iter, Real &tol, std::ostream *p_cerr=0);
  @end example

EXAMPLE:
  @noindent
  The simplest call to 'pcg' has the folling form:
  @example
    size_t max_iter = 100;
    double tol = 1e-7;
    int status = pcg(a, x, b, EYE, max_iter, tol, &cerr);
  @end example

DESCRIPTION:       
  @noindent
  @code{pcg} 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
    @itemx x
	approximate solution to Ax = b

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

    @itemx tol
	the residual after the final iteration
  @end table

NOTE: 

  @noindent
  @code{pcg} is an iterative template routine.

  @noindent
  @code{pcg} follows the algorithm described on p. 15 in

  @quotation
	Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 
	2nd Edition, 
        R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout,
	R. Pozo, C. Romine, H. Van der Vorst,
        SIAM, 1994, 
	@url{ftp.netlib.org/templates/templates.ps}.
  @end quotation

  @noindent
  The present implementation 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: 
    20 april 2009
METHODS: @pcg
End:
*/

namespace rheolef { 
//<pcg:
template < class Matrix, class Vector, class Preconditioner, class Real, class Size>
int pcg(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
	Size &max_iter, Real &tol, std::ostream *p_cerr = 0, std::string label = "cg")
{
    Vector b = M.solve(Mb);
    Real norm2_b = dot(Mb,b);
    if (norm2_b == Real(0)) norm2_b = 1;
    Vector Mr = Mb - A*x;
    Real  norm2_r = 0;
    if (p_cerr) (*p_cerr) << "[" << label << "] #iteration residue" << std::endl;
    Vector p;
    for (Size n = 0; n <= max_iter; n++) {
	Vector r = M.solve(Mr);
        Real prev_norm2_r = norm2_r;
        norm2_r = dot(Mr, r);
	if (p_cerr) (*p_cerr) << "[" << label << "] " << n << " " << ::sqrt(norm2_r/norm2_b) << std::endl;
        if (norm2_r <= sqr(tol)*norm2_b) {
          tol = ::sqrt(norm2_r/norm2_b);
          max_iter = n;
          return 0;     
        }
        if (n == 0) {
          p = r;
        } else {
          Real beta = norm2_r/prev_norm2_r;
          p = r + beta*p;
        }
        Vector Mq = A*p;
        Real alpha = norm2_r/dot(Mq, p);
        x  += alpha*p;
        Mr -= alpha*Mq;
    }
    tol = ::sqrt(norm2_r/norm2_b);
    return 1;
}
//>pcg:
}// namespace rheolef
# endif // _SKIT_PCG_H
librheolef-dev 5.93-2 / usr / include / rheolef / pcg.h
