//***************************************************************** // Iterative template routine -- CHEBY // // CHEBY solves the symmetric positive definite linear // system Ax = b using the Preconditioned Chebyshev Method // // CHEBY follows the algorithm described on p. 30 of 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 Preconditioner, class Real, class Type > int CHEBY(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol, Type eigmin, Type eigmax) { Real resid; Type alpha, beta, c, d; Vector p, q, z; Real normb = norm(b); Vector r = b - A * x; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } c = (eigmax - eigmin) / 2.0; d = (eigmax + eigmin) / 2.0; for (int i = 1; i <= max_iter; i++) { z = M.solve(r); // apply preconditioner if (i == 1) { p = z; alpha = 2.0 / d; } else { beta = c * alpha / 2.0; // calculate new beta beta = beta * beta; alpha = 1.0 / (d - beta); // calculate new alpha p = z + beta * p; // update search direction } q = A * p; x += alpha * p; // update approximation vector r -= alpha * q; // compute residual if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = i; return 0; // convergence } } tol = resid; return 1; // no convergence }