public class CG {

    /* Solve Ax = b using a preconditioned conjugate-gradient iteration.
     * 
     * Arguments:
     *  aMat - Matrix A, an nxn Sparse matrix on which matvec is called
     *  mMat - Matrix M, an nxn Sparse matrix on which psolve is called
     *  b     - System right-hand side, an n-element vector
     *  x     - An n-ement vector for storing the result
     *  rtol  - Tolerance on the relative residual (||r||/||b||)
     *  n     - System size
     *  rhist - Residual norm history.  Should be at least maxiter doubles
     *            if it is not null.
     *  maxiter - Maximum number of iterations before returning
     *
     * Returns:
     *  Iteration count on success, -1 on failure
     */

    public static int precondCG(Sparse aMat, Preconditioner pMat, 
				double [1d] b, double [1d] x,
				double rtol, int n, 
				double [1d] rhist, int maxiter)
    {
	double bnorm2;     /* ||b||^2 */
	double rnorm2;     /* Residual norm squared */
	double rz, rzold;  /* r'*z from two successive iterations */
	double alpha, beta;

	double [1d] s;  /* Search direction */
	double [1d] r;  /* Residual         */
	double [1d] z;  /* Temporary vector */

	int i;      /* Current iteration */

	s  = new double [0:n-1];
	r  = new double [0:n-1];
	z  = new double [0:n-1];

	bnorm2 = Vector.ddot(b, b);
	x.set(0);
	r.copy(b);
	pMat.psolve(z, r);
	s.copy(z);

	rz     = Vector.ddot(r, z);
	rnorm2 = Vector.ddot(r, r);

	for (i = 0; i < maxiter && rnorm2 > bnorm2 * rtol*rtol; ++i) {
	    if (rhist != null)
		rhist[i] = Math.sqrt(rnorm2/bnorm2);

	    aMat.matvec(z, s);
	    alpha = rz / Vector.ddot(s, z);
	    Vector.axpy(x,  alpha, s, x);
	    Vector.axpy(r, -alpha, z, r);

	    pMat.psolve(z, r);
	    rzold  = rz;
	    rz     = Vector.ddot(r, z);
	    rnorm2 = Vector.ddot(r, r);

	    beta = -rz/rzold;
	    Vector.axpy(s, -beta, s, z);

	}

	z = null; // Allow temporary space to be freed (not necessary)
	r = null;
	s = null;

	if (i >= maxiter) 
	    return -1;

	if (rhist != null)
	    rhist[i] = Math.sqrt(rnorm2/bnorm2);

	return i;
    }
}