Homec4science
Diffusion lammps (lammps-icms) Browse

lammps/lib/atc/CG.h
b0007d45b732lammps-icms

Last Change
Raw File

CG.h
View Options

//*****************************************************************
// Iterative template routine -- CG
//
// CG solves the symmetric positive definite linear
// system Ax=b using the Conjugate Gradient method.
//
// CG follows the algorithm described on p. 15 in the
// SIAM Templates book.
//
// 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:
//
// 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
//
//*****************************************************************
template < class Matrix, class Vector, class DataVector, class Preconditioner, class Real >
int
CG(const Matrix &A, Vector &x, const DataVector &b, const Preconditioner &M, int &max_iter, Real &tol) {
Real resid;
DenseVector<Real> p, z, q;
Real alpha, beta, rho, rho_1;
DenseVector<Real> tmp;
tmp.reset(b.size());
p.reset(b.size());
z.reset(b.size());
q.reset(b.size());
Real normb = b.norm();
DenseVector<Real> r;
tmp = A*x;
r = b - tmp;
// Implicit assumption that only diagonal matrices are being used for preconditioning
Preconditioner Minv = M.inv();
if (normb == 0.0)
normb = 1;
if ((resid = r.norm() / normb) <= tol) {
tol = resid;
max_iter = 0;
return 0;
}
for (int i = 0; i < max_iter; i++) {
z = Minv*r;
rho = r.dot(z);
if (i == 0)
p = z;
else {
beta = rho / rho_1;
tmp = p*beta;
p = z + tmp;
}
q = A*p;
alpha = rho / p.dot(q);
x += p*alpha;
r -= q*alpha;
if ((resid = r.norm() / normb) <= tol)
{
tol = resid;
max_iter = i+1;
return 0;
}
rho_1 = rho;
}
tol = resid;
return 1;
}
c4science · Help