gmres.h

Go to the documentation of this file.

/* $Id: gmres.h,v 1.10.2.23 2022/11/03 17:28:12 garloff Exp $ */

#ifndef TBCI_SOLVER_GMRES_H
#define TBCI_SOLVER_GMRES_H

#include "tbci/basics.h"

NAMESPACE_TBCI

INST2(template <Matrix<T>, Vector<T> > class NN friend \
void Update (Vector<T>&, int, Matrix<T>&, Vector<T>&, Vector<T>*);)
INST2(template <BdMatrix<T>, Vector<T> > class NN friend \
void Update (Vector<T>&, int, BdMatrix<T>&, Vector<T>&, Vector<T>*);)

                        
//Aenderung Marco:
template < typename SysMatrix, typename SysVector >
void Update(SysVector& x, int k, SysMatrix& h, SysVector& s, SysVector *V)
{
        SysVector y(s);

        // Backsolve:  
        for (int i = k; i >= 0; --i) {
                y(i) /= h(i,i);
                for (int j = i - 1; j >= 0; --j)
                        y(j) -= h(j,i) * y(i);            
        }

        for (int j = 0; j <= k; ++j)
                //Aenderung Marco:
                //x += V[j] * y(j);
                x +=  /*CPLX__ conj*/ (y(j)) * V[j];
                /*statt:
                x += v[j] * y(j);
                */
}

INST(template <typename T> class NN friend \
     void GeneratePlaneRotation(const T &, const double &, double &, T &);)

template<typename T> 
inline void GeneratePlaneRotation(const T &dx, const double &dy, double &cs, T &sn)
{
        if (UNLIKELY(dy == 0.0)) {
                cs = 1.0;
                sn = (T)0.0;
        } else if (UNLIKELY (dx == (T)0.0)) {
                cs = 0.0;
                sn = (T)1.0;
        }
#if 1
        else {  /* dy is always larger or equal than dx */
                const double rel = fabssqr(dx/dy);
                const double ynorm = 1.0 / MATH__ sqrt (1.0 + rel);
                cs = MATH__ sqrt(rel) * ynorm;
                sn = cs * dy / dx;
        }
#else
        else {
                const double rnorm = 1.0/MATH__ sqrt (fabssqr(dx) + fabssqr(dy));
                cs = MATH__ fabs(dx) * rnorm /* * sign(dx)*/;
                sn = cs * dy / dx;
        }
#endif  
}


INST(template <typename T> class NN friend \
     void ApplyPlaneRotation(T &, T &, const double &, const T &);)

template<typename T> 
inline void ApplyPlaneRotation(T &dx, T &dy, const double &cs, const T &sn)
{
        REGISTER T temp  =  /*CPLX__ conj*/(cs) * dx + CPLX__ conj (sn) * dy;
        dy = cs * dy - sn * dx;
        dx = temp;
}


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


template < typename T, typename SysMatrix, typename SysVector >
int GMRES (const SysMatrix &A, SysVector &x, const SysVector &b,
           const Preconditioner_Sig<T, SysMatrix> &M, int &m,
           unsigned int &max_iter, double &tol)
{
        const unsigned int dim = A.rows();
        if (m == 0)
                m = dim/8+max_iter/4;
        if ((unsigned)m > max_iter)
                m = max_iter;
        SysVector s(m+1), sn(m+1), w(dim), r(dim);
        Vector<double> cs(m+1);
        Matrix<T> H(m+1, m+1);
        int restart = 0;

        double normb = MATH__ sqrt(fabssqr(M.solve(b)));
        if (UNLIKELY(normb < 1e-16))
                normb = 1e-16;

        r = M.solve(b - A * x);
        
        const double rel_tol = tol*normb;

        double beta = r.fabs();

        if (beta <= rel_tol) {
                tol = beta / normb;
                max_iter = 0;
                //m = 0;
                return 0;
        }

        unsigned int j;
        int err = 0;
        
        //SysVector *V = new SysVector[m+1](dim);
        SysVector *V = new SysVector[m+1];
        for (j = 0; j <= (unsigned)m; ++j)
                V[j].resize(dim);
        
        j = 1;
        while (j <= max_iter) {
                
                V[0] = r /(T)beta;
                s.fill ((T)0.0); s(0) = (T) beta;
                //cs.fill (0.0); sn.fill((T)0.0);
                
                for (int i = 0; (i < m) && (j <= max_iter); ++i, ++j) {
                        int k;
                        
                        w = M.solve (A * V[i]);
                        
                        for (k = 0; k <= i; k++) {
                                H(k, i) = dot (V[k], w);
                                //H(k, i) = dot (w, V[k]);
                                //H(k, i) = w * V[k];
                                //H (k,i) = CPLX__ conj (w * V[k]);
                                w -= H(k, i) * V[k];
                                //w -= CPLX__ conj(H(k, i)) * V[k];
                        }
                        const REGISTER double normw = w.fabs();
                        H(i+1, i) = normw;
                        // This breakdown is a good one ...
                        //if (normw == 0.0) 
                        //      goto GMR_out;
                        if (UNLIKELY(normw == 0))
                                V[i+1] = w;
                        else    
                                V[i+1] = w / normw;
                        //for (k = 0; k < dim; ++k)
                        //      V[i+1](k) = CPLX__ conj(w(k) /normw);
                        //V[i+1] = CPLX__ conj(w) / H(i+1, i);
                        
                        for (k = 0; k < i; k++)
                                ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));

                        GeneratePlaneRotation(H(i,i), normw, cs(i), sn(i));
                        ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
                        ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
                        
#ifdef _MSC_VER
                        using ::fabs;
                        using ::std::fabs;
                        beta = fabs(s(i+1));
#else
                        beta = MATH__ fabs(s(i+1));
#endif
                        
                        if (beta <= rel_tol) {
                                Update(x, i, H, s, V);
                                goto GMR_out;
                        }
                }
                Update(x, m - 1, H, s, V);
                r = M.solve(b - A * x);
                beta = r.fabs();
                if (beta < rel_tol)
                        goto GMR_out;
                ++restart;
        }
        
        err = 1;
        
    GMR_out:
        tol = beta / normb;
        max_iter = j;
        //m = restart;
        delete[] V;
        return err;
}

NAMESPACE_END

#endif  /* TBCI_SOLVER_GMRES_H */