qmr.h

Go to the documentation of this file.

#ifndef TBCI_SOLVER_QMR_H
#define TBCI_SOLVER_QMR_H
 
#include "tbci/basics.h"

NAMESPACE_TBCI


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


//Aenderung Marco:
template < typename T,typename SysMatrix,typename SysVector >
int QMR(const SysMatrix &A, SysVector &x, const SysVector &b,
        const Preconditioner_Sig<T, SysMatrix> &M1, const Preconditioner_Sig<T, SysMatrix> &M2,
        unsigned int &max_iter, double &tol)
{
  double rho, rho_1, xi;
  // KG, 2001-06-27: Probably some of these should be double instead of complex
  //  resulting in fabs to be used instead of abs. Maybe fabssqr() as optimization?
  // KG, 2002-07-24: OK, converted things ...
  //  The solver now seems to work for complex variables.
  double gamma, gamma_1;
  //typename SysVector::value_type gamma, gamma_1;
  typename SysVector::value_type theta, theta_1;
  //double theta, theta_1;
  typename SysVector::value_type eta, delta, beta;
  typename SysVector::value_type ep = (typename SysVector::value_type)0;

  const unsigned int dim = A.rows();
  // KG, 2001-06-27: A few of those should be optimised away ...
  SysVector r(dim), v_tld(dim), y(dim), w_tld(dim), z(dim);
  SysVector v(dim), w(dim), y_tld(dim), z_tld(dim);
  SysVector p(dim), q(dim), p_tld(dim), d(dim), s(dim);

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

  r = b - A * x;

#define QMR_OUT(r) do { err = r; goto qmr_out; } while(0)
        
  unsigned i = 0;
  if ((residsqr = r.fabssqr() * normbsqr) <= tolsqr)
        QMR_OUT (0);

  //M1.update(A);
  //M2.update(A);
        
  v_tld = r;
  y = M1.solve(v_tld);
  rho = /*(typename SysVector::value_type)*/ y.fabs();

  w_tld = r;
  z = M2.transSolve(w_tld);
  xi = /*(typename SysVector::value_type)*/ z.fabs();

  gamma = /*(typename SysVector::value_type)*/ 1.0;
  eta   = (typename SysVector::value_type)-1.0;
  theta = (typename SysVector::value_type) 0.0;

  for (i = 1; i <= max_iter; i++)
  {
    if (rho == 0.0)
        QMR_OUT(2);                     // return on breakdown

    if (xi == 0.0)
        QMR_OUT(7);                     // return on breakdown

    v = v_tld * (1. / rho);
    y = y     * (1. / rho);

    w = w_tld * (1. / xi);
    z = z     * (1. / xi);

    //delta = dot(z, y);
    delta = z * y;
    if (delta == 0.0)
        QMR_OUT(5);                     // return on breakdown

    y_tld = M2.solve(y);                // apply preconditioners
    z_tld = M1.transSolve(z);

    if (i > 1) {
        p = y_tld - (typename SysVector::value_type)(xi * delta / ep) * p;
        q = z_tld - (typename SysVector::value_type)(rho* delta / ep) * q;
    } else {
        p = y_tld;
        q = z_tld;
    }

    p_tld = A * p;
    //ep = dot(q, p_tld);
    ep = q * p_tld;
    if (ep == 0.0) 
        QMR_OUT(6);                     // return on breakdown

    beta = ep / delta;
    if (beta == 0.0) 
        QMR_OUT(3);                     // return on breakdown

    v_tld = p_tld - beta * v;
    //v_tld = p_tld - CPLX__ conj(beta) * v;
    y = M1.solve(v_tld);

    rho_1 = rho;
    rho   = /*(typename SysVector::value_type)*/ y.fabs();
    w_tld = A.transMult(q) - beta * w; 
    //w_tld = A.transMult(q) - CPLX__ conj(beta) * w; 
    z = M2.transSolve(w_tld);

    xi = /*(typename SysVector::value_type)*/ z.fabs();

    gamma_1 = gamma;
    theta_1 = theta;
    theta = rho / (gamma_1 * beta);
    //theta = rho / dot(gamma_1, beta);
          
#ifdef _MSC_VER         /*  STD__ != MATH__ */
//    typename SysVector::value_type (1.0 + theta*theta); // dot (theta, theta) ??
      double c(1.0 + fabssqr(theta));
        { 
                using std::sqrt; using ::sqrt; 
                gamma = sqrt(c); // dot (theta, theta) ??
        }
#else
# ifdef HAVE_LIBC_NEQ_CPLX_BUG
#  if !defined(_SGI_SOURCE) || !defined(_COMPILER_VERSION)
    using STD__ sqrt; 
#  endif    
    using MATH__ sqrt;
    gamma = 1.0 / sqrt(1.0 + fabssqr(theta)); // dot (theta, theta) ??
# else
    gamma = 1.0 / MATH__ sqrt(1.0 + fabssqr(theta)); // dot (theta, theta) ??
# endif
#endif

    if (gamma == 0.0) 
        QMR_OUT(4);                     // return on breakdown

#if defined(USE_BUILTIN_CPLX) && defined(TYPE)
    eta = -hcplx(eta)* rho_1 * gamma * gamma /
          (beta * gamma_1 * gamma_1);
#else     
    eta = -eta * rho_1 * gamma * gamma / 
          (beta * gamma_1 * gamma_1);
#endif    

    if (i > 1) {
        const typename SysVector::value_type _tmp = fabssqr(theta_1) * fabssqr(gamma);
        d = eta * p     + _tmp * d;
        s = eta * p_tld + _tmp * s;
    } else {
        d = eta * p;
        s = eta * p_tld;
    }

    x += d;                            // update approximation vector
    r -= s;                            // compute residual

    if ((residsqr = r.fabssqr() * normbsqr) <= tolsqr) 
          QMR_OUT(0);
  }

 qmr_out:
  max_iter = i;
  tol = MATH__ sqrt(residsqr);
  return err;                            // no convergence
}

NAMESPACE_END

#endif  /* TBCI_SOLVER_QMR_H */