cheby.h

Go to the documentation of this file.

template < class Matrix, class Vector, class Preconditioner, class Real,
class Type >
int 
CHEBY(const Matrix &A, Vector &x, const Vector &b,
      const Preconditioner &M, int64_t &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 (int64_t 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
}