bicg.h

Go to the documentation of this file.

/******************************************************************
 * @ingroup solver
 * @brief BiCG solves the unsymmetric linear system \f$ Ax = b \f$
 * using the Preconditioned BiConjugate Gradient method
 * @return The return value indicates convergence within max_iter (input)
 * iterations (0), or no convergence within max_iter iterations (1). \n
 * Upon successful return, output arguments have the following values:
 * @param A  -- matrix of the system
 * @param b  -- vector of the system
 * @param M  -- preconditioner matrix
 * @param x  --  approximate solution to Ax = b
 * @param max_iter  --  the number of iterations performed before the tolerance was reached
 * @param tol  --  the residual after the final iteration
 *****************************************************************/
template < class Matrix, class Vector, class Preconditioner, class Real >
int 
BiCG( Matrix &A, Vector &x, const Vector &b,
     Preconditioner &M, int64_t &max_iter, Real &tol)
{
                Real resid;
                Vector rho_1(1), rho_2(1), alpha(1), beta(1);
                Vector z, ztilde, p, ptilde, q, qtilde;
                
                Real normb = TPZExtractVal::val(Norm(b));
                Vector r = b - A * x;
                Vector rtilde = r;
                
                if (normb == 0.0)
                                normb = 1;
                
                if ((resid = (TPZExtractVal::val(Norm(r)) / normb) <= tol)) {
                                tol = resid;
                                max_iter = 0;
                                return 0;
                }
                
                for (int64_t i = 1; i <= max_iter; i++) {
                                M.Solve(r,z);
                                M.Solve(rtilde,ztilde);
                                rho_1(0) = Dot(z, rtilde);
                                if (TPZExtractVal::val(rho_1(0)) == ((Real)0.)) {
                                                tol = (TPZExtractVal::val(Norm(r))) / normb;
                                                max_iter = i;
                                                return 2;
                                }
                                if (i == 1) {
                                                p = z;
                                                ptilde = ztilde;
                                } else {
                                                beta(0) = rho_1(0) / rho_2(0);
                                                p.TimesBetaPlusZ(beta(0),z);
                                                ptilde.TimesBetaPlusZ(beta(0),ztilde);
                                }
                                q = A * p;
                                A.Multiply(ptilde,qtilde,1);
                                alpha(0) = rho_1(0) / Dot(ptilde, q);
                                x += alpha(0) * p;
                                r -= alpha(0) * q;
                                rtilde -= alpha(0) * qtilde;
                                
                                rho_2(0) = rho_1(0);
                                if ((resid = (TPZExtractVal::val(Norm(r))) / normb) < tol) {
                                                tol = resid;
                                                max_iter = i;
                                                return 0;
                                }
                }
                
                tol = resid;
                return 1;
}