preconditioned conjugate gradient solver

[Cplusc]

I wrote a preconditioned conjugate gradient solver but when calling the function it just gives me zero. any help would be appreciated.

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

[Handy Andy]

Hello resabzr,

but when calling the function it just gives me zero

It take this to mean that return converged; is returning zero?

Without the rest of the code and anything else that goes with the program there is no way for me to test it. This also means any input file the program might use.

The variable "converged" is changes in 2 places each within an if statement. So the first thing to do is check to see if either if statement is ever true.

If the function returns zero that tells me that neither if statement is true and "converged" is never changed.

I would hope that you are compiling to at least the C++11 standards. To make your life a little easier there is the uniform initializer, the {}s.

You could rewrite some of your code as:

double rho_1{}, rho{}; // stores (z _ '* r_) at subsequent iterations
double alpha{}; // solution increment along the search direction
double beta{};
double resid{}, resid2{}, maxincr{}, maxval{};
double normb{};
int iter{ 1 }, ret{ 1 };
int converged{};

An empty set of {} will allow the compiler to choose the correct form of zero based on the variable's type.

And a line like if (rho == 0) can be shortened to if (!rho)
Just some tips that you might find useful.

Andy

[Cplusc]

Hello Handy Andy,
Thanks for response and good points. I meant that the solution is a vector in which all elements are zero.

[TheIdeasMan]

Hi resabzr,

After all your time here, have you yet learnt to use a debugger?

if (rho == 0)

rho is of type double, so this will be problematic. What if rho == 1e-15 ? I always write double literal values with digits on both sides of the decimal point, so if (rho == 0.0) becomes an obvious red flag.

Why is converged not of type bool?

VecDbl_t diag(x.size());

You have 7 of these: why not assign x.size() to a variable and use that, rather than call the same function 7 times?

Line 56: Is rho_1 zero for the first iteration? divide by zero error? Always check that denominators are not zero or close enough to zero to cause a problem.

normb = dotproduct(b, b);

This looks suspicious. Is it supposed to be normalised b ? I have always disliked the term normalised in the context of vectors, IMO it should be unit vector - a vector of magnitude 1.0. That is not what you have.

[TheIdeasMan]

What does the ctor of VecDbl_t do?

if (rho == 0.0 ) could be true IF AND ONLY IF both z and r are zero vectors, which might be true if the ctor of VecDbl_t sets each value to zero. If this is the case, then large sections of the code are a furphy.

[lastchance]

Your code will only ever do one iteration, and it won't change sol() on that iteration.


L14 - VecDbl_t z(x.size()); presumably sets z to 0 (deliberately or otherwise).

The next thing you do to z is on L35 when you divide by the matrix's diagonal coefficients on the first iteration; z remains as zero.

L38 - rho += dotproduct(z, r); then simply leaves rho as 0 ...

... so it enters the first if block:

	if (rho == 0) {
...
			converged = true;

and as you return at the end of the first iteration anyway (Line 101 - return converged;) you never get to change your solution vector.

I'm guessing that sol() on entry probably contained your obligatory initial guess and is supposed to be the answer you require at the end, but you don't say so.



Either give enough code for people to run it, or provide the equations / algorithm that you are trying to invoke. There are many different ways of preconditioning the conjugate gradient algorithm.