42#ifndef STOKHOS_TPETRA_CG_HPP
43#define STOKHOS_TPETRA_CG_HPP
48template <
typename Matrix,
51bool CG_Solve(
const Matrix& A, Vector& x,
const Vector& b,
52 typename Vector::mag_type tol,
Ordinal max_its,
53 std::ostream* out = 0)
55 typedef typename Vector::mag_type mag_type;
56 typedef typename Vector::dot_type dot_type;
62 r.update(1.0, b, -1.0);
64 dot_type rho = r.dot(r);
65 dot_type rho_old = rho;
66 mag_type nrm = std::sqrt(rho);
68 dot_type alpha, beta, pAp;
69 while (k < max_its && nrm > tol) {
71 p.update(1.0, r, 0.0);
75 p.update(1.0, r, beta);
80 x.update( alpha, p, 1.0);
81 r.update(-alpha, w, 1.0);
89 *out << k <<
": " << nrm << std::endl;
98template <
typename Matrix,
102bool PCG_Solve(
const Matrix& A, Vector& x,
const Vector& b,
const Prec& M,
103 typename Vector::mag_type tol,
Ordinal max_its,
104 std::ostream* out = 0)
106 typedef typename Vector::mag_type mag_type;
107 typedef typename Vector::dot_type dot_type;
109 Vector r(b.getMap());
110 Vector p(x.getMap());
111 Vector w(x.getMap());
113 r.update(1.0, b, -1.0);
115 mag_type nrm = r.norm2();
118 dot_type rho_old, alpha, beta, pAp;
119 while (k < max_its && nrm > tol) {
125 p.update(1.0, w, 0.0);
128 beta = rho / rho_old;
129 p.update(1.0, w, beta);
134 x.update( alpha, p, 1.0);
135 r.update(-alpha, w, 1.0);
141 *out << k <<
": " << nrm << std::endl;
Top-level namespace for Stokhos classes and functions.
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
bool PCG_Solve(const Matrix &A, Vector &x, const Vector &b, const Prec &M, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)