SIMPLE SOLUTIONS

PCG(4RHEOLEF) - Linux man page online | Special files

Chapter
rheolef-6.7
pcg(4rheolef) rheolef-6.7 pcg(4rheolef)

NAME

pcg -- conjugate gradient algorithm.

SYNOPSIS

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);

EXAMPLE

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

DESCRIPTION

pcg solves the symmetric positive definite linear system Ax=b using the Conjugate Gradient method. The return value indicates convergence within max_iter (input) iterations (0), or no con‐ vergence within max_iter iterations (1). Upon successful return, output arguments have the following values: x approximate solution to Ax = b max_iter the number of iterations performed before the tolerance was reached tol the residual after the final iteration

NOTE

pcg is an iterative template routine. pcg follows the algorithm described on p. 15 in 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, ftp.netlib.org/templates/templates.ps. The present implementation is inspired from IML++ 1.2 iterative method library, http://math.nist.gov/iml++.

IMPLEMENTATION

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; }
rheolef-6.7 rheolef-6.7 pcg(4rheolef)
This manual Reference Other manuals
pcg(4rheolef) referred by dia(2rheolef)
refer to
Download raw manual
Index rheolef-6.7 (+82) rheolef-6.7 (+82) № 4 (+981)
Go top