1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

int PCG(SparseMatrix& A, VecDbl_t& x, VecDbl_t& sol, size_t& max_iter, double & tol, double & tol2) {
double rho_1(0.0), rho(0.0); // stores (z _ '* r_) at subsequent iterations
double alpha(0.0); // solution increment along the search direction
double beta(0.0);
double resid(0.0), resid2(0.0), maxincr(0.0), maxval(0.0);
double normb(0.0);
int iter(1), ret(1);
int converged = 0;
VecDbl_t diag(x.size());
VecDbl_t b(x.size());
VecDbl_t r(x.size());
VecDbl_t z(x.size());
VecDbl_t p(x.size());
VecDbl_t q(x.size());
VecDbl_t aux(x.size());
for (size_t i = 0; i < diag.size(); i++) {
diag[i] = A.val(i, i);
}
b = x;
normb = dotproduct(b, b);
MxV(A, sol, aux);
for (size_t i = 0; i < diag.size(); i++) {
r[i] = x[i]  aux[i];
}
while (iter <= max_iter) // CG iteration loop
{
iter++;
for (size_t i = 0; i < diag.size(); i++) {
z[i] = z[i] / diag[i];
}
rho += dotproduct(z, r);
if (rho == 0) {
resid += dotproduct(r, r);
resid = std::sqrt(resid / normb);
if (resid <= tol) {
ret = 0;
}
else {
ret = 2;
}
converged = true;
}
else {
if (iter == 1) {
p = z;
MxV(A, p, q);
}
else {
beta = rho / rho_1;
for (size_t i = 0; i < diag.size(); i++) {
p[i] = z[i] + beta * p[i];
}
MxV(A, p, q);
}
alpha += dotproduct(p, q);
alpha = rho / alpha; // alpha = rho / (p _ '* q_)
resid = 0.0, maxincr = 0.0, maxval = 0.0;
double maxincr_i(0.0), maxval_i(0.0), resid_i(0.0);
{
double dx, fabs_dx, fabs_x, sub_rho(0.0);
for (size_t ii = 0; ii < x.size(); ++ii) {
dx = alpha * p[i];
sol[i] += dx;
r[i] = alpha * q[i];
sub_rho += r[i] * r[i];
fabs_dx = std::fabs(dx); fabs_x = std::fabs(sol[i]);
fabs_dx = std::fabs(dx); fabs_x = std::fabs(sol[i]);
if (fabs_dx > maxincr_i) maxincr_i = fabs_dx;
if (fabs_x > maxval_i) maxval_i = fabs_x;
}
resid_i += sub_rho;
}
int num_procs, myrank;
resid = std::sqrt(resid / normb);
resid2 = maxincr / (maxval + tol);
if ((resid <= tol) && (resid2 <= tol2)) {
converged = 1;
tol = resid;
tol2 = resid2;
max_iter = iter;
break;
}
rho_1 = rho;
rho = 0.0; alpha = 0.0;
}
return converged;
}
}
 