librheolef-dev 6.5-1+b1 / usr / include / rheolef / pgmres.h

# ifndef _SKIT_GMRES_H
# define _SKIT_GMRES_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
/// 
/// =========================================================================

#ifdef HAVE_CONFIG_H
#include "rheolef/compiler.h"
#endif
#include <cmath> 

/*Class:pgmres
NAME: @code{pgmres} -- generalized minimum residual method 
@findex qmr
@cindex generalized minimum residual method
@cindex iterative solver
SYNOPSIS:
  @example
   template <class Matrix, class Vector, class Preconditioner,
     class SmallMatrix, class SmallVector, class Real, class Size>
   int pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
     SmallMatrix &H, const SmallVector& dummy,
     Size m, Size &max_iter, Real &tol);
  @end example
EXAMPLE:
  @noindent
  The simplest call to @code{pgmres} has the folling form:
  @example
        double tol = 1e-7;
        size_t max_iter = 100;
  	size_t m = 6;
  	boost::numeric::ublas::matrix<double> H(m+1,m+1);
  	vec<double,sequential> dummy;
  	int status = pgmres (a, x, b, ic0(a), H, dummy, m, max_iter, tol);
  @end example
DESCRIPTION:       
  @noindent
  @code{pgmres} solves the unsymmetric linear system Ax = b 
  using the generalized minimum residual 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
  @noindent
  In addition, M specifies a preconditioner, H specifies a matrix
  to hold the coefficients of the upper Hessenberg matrix constructed
  by the @code{pgmres} iterations, @code{m} specifies the number of iterations
  for each restart.

  @noindent
  @code{pgmres} requires two matrices as input, A and H.
  The matrix A, which will typically be a sparse matrix) corresponds
  to the matrix in the linear system Ax=b.
  The matrix H, which will be typically a dense matrix, corresponds
  to the upper Hessenberg matrix H that is constructed during the 
  @code{pgmres} iterations. Within @code{pgmres}, H is used in a different way
  than A, so its class must supply different functionality.
  That is, A is only accessed though its matrix-vector and 
  transpose-matrix-vector multiplication functions.
  On the other hand, @code{pgmres} solves a dense upper triangular linear 
  system of equations on H. Therefore, the class
  to which H belongs must provide H(i,j) operator for element acess.

NOTE: 
  @noindent
  It is important to remember that we use the convention that indices
  are 0-based. That is H(0,0) is the first component of the
  matrix H. Also, the type of the matrix must be compatible with the
  type of single vector entry. That is, operations such as
  H(i,j)*x(j) must be able to be carried out.
  
  @noindent
  @code{pgmres} is an iterative template routine.

  @noindent
  @code{pgmres} follows the algorithm described on p. 20 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
    LMC-IMAG, 38041 Grenoble cedex 9, France

DATE: 
    
    12 march 1997

METHODS: @pgmres
End:
*/
// ========== [ method implementation ] ====================================
namespace rheolef { 

//<pgmres:
template <class SmallMatrix, class Vector, class SmallVector, class Size>
void 
Update(Vector &x, Size k, SmallMatrix &h, SmallVector &s, Vector v[])
{
  SmallVector y = s;
  // Back solve:  
  for (int i = k; i >= 0; i--) {
    y(i) /= h(i,i);
    for (int j = i - 1; j >= 0; j--)
      y(j) -= h(j,i) * y(i);
  }
  for (Size j = 0; j <= k; j++) {
    x += v[j] * y(j);
  }
}
template<class Real> 
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
  if (dy == Real(0)) {
    cs = 1.0;
    sn = 0.0;
  } else if (abs(dy) > abs(dx)) {
    Real temp = dx / dy;
    sn = 1.0 / ::sqrt( 1.0 + temp*temp );
    cs = temp * sn;
  } else {
    Real temp = dy / dx;
    cs = 1.0 / ::sqrt( 1.0 + temp*temp );
    sn = temp * cs;
  }
}
template<class Real> 
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
  Real temp  =  cs * dx + sn * dy;
  dy = -sn * dx + cs * dy;
  dx = temp;
}
template <class Matrix, class Vector, class Preconditioner,
          class SmallMatrix, class SmallVector, class Real, class Size>
int
pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
      SmallMatrix &H, const SmallVector&, const Size &m, Size &max_iter, Real &tol)
{
  odiststream* p_derr = &derr;
  std::string label = "pgmres";
  if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
  Vector w;
  SmallVector s(m+1), cs(m+1), sn(m+1);
  Size i;
  Size j = 1;
  Size k;
  Real resid;
  Real normb = norm(M.solve(b));
  Vector r = M.solve(b - A * x);
  Real beta = norm(r);
  if (normb == Real(0)) {
    normb = 1;
  } 
  resid = norm(r) / normb;
  if (resid  <= tol) {
    tol = resid;
    max_iter = 0;
    return 0;
  }
  Vector *v = new Vector[m+1];
  while (j <= max_iter) {
    v[0] = r * (1.0 / beta);    // ??? r / beta
    s = 0.0;
    s(0) = beta;
    for (i = 0; i < m && j <= max_iter; i++, j++) {
      w = M.solve(A * v[i]);
      for (k = 0; k <= i; k++) {
        H(k, i) = dot(w, v[k]);
        w -= H(k, i) * v[k];
      }
      H(i+1, i) = norm(w);
      v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
      for (k = 0; k < i; k++) {
        ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
      }
      GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
      ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
      ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
      resid = abs(s(i+1)) / normb;
      if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
      if (resid  < tol) {
        Update(x, i, H, s, v);
        tol = resid;
        max_iter = j;
        delete [] v;
        return 0;
      }
    }
    Update(x, m - 1, H, s, v);
    r = M.solve(b - A * x);
    beta = norm(r);
    resid = beta / normb;
    if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
    if (resid < tol) {
      tol = resid;
      max_iter = j;
      delete [] v;
      return 0;
    }
  }
  tol = resid;
  delete [] v;
  return 1;
}
//>pgmres:
}// namespace rheolef
# endif // _SKIT_GMRES_H