ir.h

Go to the documentation of this file.

template < class Matrix, class Vector, class Preconditioner, class Real >
int IR( Matrix &A, Vector &x,const Vector &b,
   Preconditioner &M, Vector *residual, int64_t &max_iter, Real &tol,const int FromCurrent)
{
                Real resid;
                Vector z;
                
    Real normb = TPZExtractVal::val(Norm(b));
                Vector resbackup;
                Vector *res = residual;
                if(!res) res = &resbackup;
                Vector &r = *res;

                if(FromCurrent) A.MultAdd(x,b,r,-1.,1.);
                else {
                                x.Zero();
                                r = b;
                }
                
                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);
                                x += z;
                                A.MultAdd(x,b,r,-1.,1.);
                                
                                if ((resid = TPZExtractVal::val(Norm(r)) / normb) <= tol) {
                                                tol = resid;
                                                max_iter = i;
                                                return 0;
                                }
                }
                
                tol = resid;
                return 1;
}