bicg.h

Go to the documentation of this file.

/* $Id: bicg.h,v 1.8.2.6 2019/05/28 11:13:02 garloff Exp $ */

#ifndef TBCI_SOLVER_BICG_H
#define TBCI_SOLVER_BICG_H

#include "tbci/basics.h"

#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/solver/bicg_gd.h"
#endif

NAMESPACE_TBCI

INST3(template <typename T, Matrix<T>, Vector<T> > class NN friend \
int BiCG(const Matrix<T>&, Vector<T>&, const Vector<T>&,\
         const Preconditioner_Sig<T, Matrix<T> >&, unsigned int&, double&);)
INST3(template <typename T, BdMatrix<T>, Vector<T> > class NN friend \
int BiCG(const BdMatrix<T>&, Vector<T>&, const Vector<T>&,\
         const Preconditioner_Sig<T, BdMatrix<T> >&, unsigned int&, double&);)

//Aenderung Marco:
template <typename T, typename SysMatrix, typename SysVector>
int BiCG(const SysMatrix& A, SysVector& x, const SysVector& b,
         const Preconditioner_Sig<T, SysMatrix>& M, unsigned int& max_iter, double& tol)
{
        unsigned int dim = A.rows();
        typename SysVector::value_type  rho_1, rho_2(0), alpha, beta;
        SysVector z(dim), ztilde(dim), p(dim), ptilde(dim), q(dim), qtilde(dim);

        double residsqr;
        const double tolsqr = sqr(tol);
        double normbsqr = b.fabssqr();
                
        if (normbsqr < 1e-32)  
                normbsqr = 1e32;
        else 
                normbsqr = 1.0 / normbsqr;

        SysVector r(b - A * x);
        SysVector rtilde(conj(r));

        if ((residsqr = r.fabssqr() * normbsqr) <= tolsqr) {
                tol = MATH__ sqrt(residsqr);
                max_iter = 0;
                return 0;
        }

        for (unsigned int i = 1; i <= max_iter; ++i) {
                z = M./*trans*/solve(r);
                ztilde = M.transSolve(rtilde);
                rho_1 = z * rtilde;
                if (rho_1 == (typename SysVector::value_type)0) {
                        tol = MATH__ sqrt(r.fabssqr() * normbsqr);
                        max_iter = i;
                        return 2;
                }
                if (UNLIKELY(i == 1)) {
                        p = z;
                        ptilde = ztilde;
                } else {
                        beta = rho_1 / rho_2;
                        p = z + beta * p;
                        ptilde = ztilde + beta * ptilde;
                }

                q = A * p;
                qtilde = A.transMult(ptilde);
                alpha = rho_1 / (q * ptilde);
                x += alpha * p;
                r -= alpha * q;
                rtilde -= alpha * qtilde;
                
                rho_2 = rho_1;
                if ((residsqr = r.fabssqr() * normbsqr) < tolsqr) {
                        tol = MATH__ sqrt(residsqr);
                        max_iter = i;
                        return 0;
                }
        }

        tol = MATH__ sqrt(residsqr);
        return 1;
}

NAMESPACE_END

#endif  /* TBCI_SOLVER_BICG_H */