/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
/// 
/// =========================================================================

#include "rheolef/vec_expr_v2.h"

namespace rheolef {

/*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, odiststream *p_derr=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, &derr);
  @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
    @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{pcg} is an iterative template routine.

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

  @emph{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}.

  @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:
*/

//<pcg:
template <class Matrix, class Vector, class Vector2, class Preconditioner, class Real, class Size>
int pcg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
	Size &max_iter, Real &tol, odiststream *p_derr = 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_derr) (*p_derr) << "[" << 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_derr) (*p_derr) << "[" << 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 6.7-6 / usr / include / rheolef / pcg.h
