bicgstab.h

Go to the documentation of this file.

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

#ifndef TBCI_SOLVER_BICGSTAB_H
#define TBCI_SOLVER_BICGSTAB_H  

#include "tbci/vector.h"
#include "tbci/solver/precond.h"

#ifdef SOLVER_LOG
# include <fstream>
# include <iomanip>
#endif

#ifndef NO_BICGSTAB_RESTART
# define BICGSTAB_RESTART
#endif

/* There was some rather obscure magic to improve convergence, which is 
 * removed now. (KG, 99/06/15)
 * Unlike the CGS solver, which does weird changes to the vector before 
 * finding a point from which to converge rapidly, the BiCGstab is unlikely
 * to ever deteriorate the Vector. Drawback is that the approach rate is
 * getting lower and lower, when the number of iterations increase.
 * Now the magic comes to play. Whenever we are in such a situation, i.e.
 * the result isn't improving any longer, we try to detect this ,,stall``
 * and recalculate some values resulting in a kick to the solver. This
 * often results in better results.
 * KG, 99/05/26: I've been talking to experts. The loss in the convergence 
 * rate is likely to be caused by the loss of the orthogonality of the two
 * base vectors. This should be detectable and solvable.
 * The magic mentioned above should be removed and replaced by some proper
 * mechanism.
 */

NAMESPACE_TBCI

/* NOTE: Template parameter T and T_tol should really be double (or 
 * similar: float, long double are OK), as fabs() is used, which 
 * returns double. This could be changed to abs(), but will then 
 * cause some trouble for complex numbers.
 */

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


template < typename T, typename SysMatrix, typename T_tol /*= double*/>
int BiCGSTAB (const SysMatrix &A, Vector<T> &x, const Vector<T> &b,
              const Preconditioner_Sig<T, SysMatrix> &M, unsigned int& max_iter,
              T_tol& tol, const unsigned off = 0)
{
        const unsigned int dim = A.rows();
        typename SysMatrix::value_type rho_1, rho_2((typename SysMatrix::value_type)1),
                 alpha((typename SysMatrix::value_type)0), beta, 
                 omega((typename SysMatrix::value_type)0);
        //T rho_1, rho_2, alpha, beta, omega;
        Vector<T> p(dim), phat(dim), s(dim), shat(dim), t(dim), v(dim);

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

        Vector<T> r (b - A * x);
        Vector<T> rtilde (r);
        //M.update (A);

        if (LIKELY((residsqr = r.fabssqr() * normbsqr) <= tolsqr)) {
                tol = MATH__ sqrt(residsqr);
                max_iter = 0;
                return 0;
        }
        alpha = beta = omega = rho_1 = (typename SysMatrix::value_type)0;
#ifdef SOLVER_LOG
        STD__ ofstream cnvg;
        if (off) cnvg.open ("bicgstab_cnvg.gnu", STD__ ios::app);
        else cnvg.open ("bicgstab_cnvg.gnu");
        cnvg << "#  Iter     res     alpha     beta    omega    rho_1     rho_2   rho_1 / sqr(normb)" << STD__ endl;
#endif

        for (unsigned int i = 1; i <= max_iter; ++i) {
#ifdef SOLVER_LOG
                cnvg << i+off << "\t" << MATH__ sqrt(residsqr) << "\t" << MATH__ fabs(alpha) << "\t" << MATH__ fabs(beta) << "\t"
                     << MATH__ fabs (omega) << "\t" << MATH__ fabs(rho_1) << "\t" << MATH__ fabs(rho_2) << "\t" << MATH__ fabs(rho_1)*normbsqr << "\n";
#endif
                rho_1 = dot (rtilde, r);
                if (UNLIKELY(rho_1 == (typename SysMatrix::value_type)0)) {
#ifdef BICGSTAB_RESTART
# ifdef BICGSTAB_REPORT_RESTART
                        STD__ cerr << "BiCGstab failure 2 (iter " << i 
                                   << ") rho1 == " << rho_1 << ". Restart.\n";
# endif
                        // RESTART
                        rtilde = b - A*x; r = rtilde; 
                        alpha = 0;
                        continue;
#else

                        tol = MATH__ sqrt (residsqr);
                        return 2;
#endif
                }
                if (UNLIKELY(i == 1))
                        p = r;
                else {
                        beta = (rho_1 * alpha) / (rho_2 * omega);
                        p = r + beta * (p - omega * v);
                }
                phat = M.solve(p);
                v = A * phat;
                alpha = rho_1 / dot (rtilde, v);
                s = r - alpha * v;
                if ((residsqr = s.fabssqr() * normbsqr) < tolsqr) {
                        x += alpha * phat; 
                        max_iter = i; tol = MATH__ sqrt(residsqr);
                        return 0;
                }
                shat = M.solve(s);
                t = A * shat;
                omega = dot(t,s) / fabssqr(t); //dot(t,t);

                x += alpha * phat + omega * shat;
                r = s - omega * t;

                rho_2 = rho_1;

                if ((residsqr = r.fabssqr() * normbsqr) < tolsqr) {
                        tol = MATH__ sqrt(residsqr); max_iter = i;
#ifdef SOLVER_LOG
                        cnvg << i+off << "\t" << MATH__ sqrt(residsqr) << "\t" << MATH__ fabs(alpha) << "\t" << MATH__ fabs(beta) << "\t"
                             << MATH__ fabs (omega) << "\t" << MATH__ fabs(rho_1) << "\t" << MATH__ fabs(rho_2) << "\t" << MATH__ fabs(rho_1)*normbsqr << STD__ endl;
#endif
                        return 0;
                }
                if (UNLIKELY(omega == (typename SysMatrix::value_type)0)) {
                        tol = MATH__ sqrt(residsqr);
                        STD__ cerr << "BiCGstab failure 3: omega == 0\n";
                        max_iter = i;
                        return 3; 
                }
        }

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

NAMESPACE_END

#endif  /* TBCI_SOLVER_BICGSTAB_H */