Actual source code: tcqmr.c

      2: /*
      3:     This file contains an implementation of Tony Chan's transpose-free QMR.

      5:     Note: The vector dot products in the code have not been checked for the
      6:     complex numbers version, so most probably some are incorrect.
      7: */

      9: #include <../src/ksp/ksp/impls/tcqmr/tcqmrimpl.h>

     11: static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
     12: {
     13:   PetscReal      rnorm0,rnorm,dp1,Gamma;
     14:   PetscScalar    theta,ep,cl1,sl1,cl,sl,sprod,tau_n1,f;
     15:   PetscScalar    deltmp,rho,beta,eptmp,ta,s,c,tau_n,delta;
     16:   PetscScalar    dp11,dp2,rhom1,alpha,tmp;

     20:   ksp->its = 0;

     22:   KSPInitialResidual(ksp,x,u,v,r,b);
     23:   VecNorm(r,NORM_2,&rnorm0);          /*  rnorm0 = ||r|| */
     24:   KSPCheckNorm(ksp,rnorm0);
     25:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm0;
     26:   else ksp->rnorm = 0;
     27:   (*ksp->converged)(ksp,0,ksp->rnorm,&ksp->reason,ksp->cnvP);
     28:   if (ksp->reason) return(0);

     30:   VecSet(um1,0.0);
     31:   VecCopy(r,u);
     32:   rnorm = rnorm0;
     33:   tmp   = 1.0/rnorm; VecScale(u,tmp);
     34:   VecSet(vm1,0.0);
     35:   VecCopy(u,v);
     36:   VecCopy(u,v0);
     37:   VecSet(pvec1,0.0);
     38:   VecSet(pvec2,0.0);
     39:   VecSet(p,0.0);
     40:   theta = 0.0;
     41:   ep    = 0.0;
     42:   cl1   = 0.0;
     43:   sl1   = 0.0;
     44:   cl    = 0.0;
     45:   sl    = 0.0;
     46:   sprod = 1.0;
     47:   tau_n1= rnorm0;
     48:   f     = 1.0;
     49:   Gamma = 1.0;
     50:   rhom1 = 1.0;

     52:   /*
     53:    CALCULATE SQUARED LANCZOS  vectors
     54:    */
     55:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
     56:   else ksp->rnorm = 0;
     57:   (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
     58:   while (!ksp->reason) {
     59:     KSPMonitor(ksp,ksp->its,ksp->rnorm);
     60:     ksp->its++;

     62:     KSP_PCApplyBAorAB(ksp,u,y,vtmp); /* y = A*u */
     63:     VecDot(y,v0,&dp11);
     64:     KSPCheckDot(ksp,dp11);
     65:     VecDot(u,v0,&dp2);
     66:     alpha  = dp11 / dp2;                          /* alpha = v0'*y/v0'*u */
     67:     deltmp = alpha;
     68:     VecCopy(y,z);
     69:     VecAXPY(z,-alpha,u); /* z = y - alpha u */
     70:     VecDot(u,v0,&rho);
     71:     beta   = rho / (f*rhom1);
     72:     rhom1  = rho;
     73:     VecCopy(z,utmp);    /* up1 = (A-alpha*I)*
     74:                                                 (z-2*beta*p) + f*beta*
     75:                                                 beta*um1 */
     76:     VecAXPY(utmp,-2.0*beta,p);
     77:     KSP_PCApplyBAorAB(ksp,utmp,up1,vtmp);
     78:     VecAXPY(up1,-alpha,utmp);
     79:     VecAXPY(up1,f*beta*beta,um1);
     80:     VecNorm(up1,NORM_2,&dp1);
     81:     KSPCheckNorm(ksp,dp1);
     82:     f      = 1.0 / dp1;
     83:     VecScale(up1,f);
     84:     VecAYPX(p,-beta,z);   /* p = f*(z-beta*p) */
     85:     VecScale(p,f);
     86:     VecCopy(u,um1);
     87:     VecCopy(up1,u);
     88:     beta   = beta/Gamma;
     89:     eptmp  = beta;
     90:     KSP_PCApplyBAorAB(ksp,v,vp1,vtmp);
     91:     VecAXPY(vp1,-alpha,v);
     92:     VecAXPY(vp1,-beta,vm1);
     93:     VecNorm(vp1,NORM_2,&Gamma);
     94:     KSPCheckNorm(ksp,Gamma);
     95:     VecScale(vp1,1.0/Gamma);
     96:     VecCopy(v,vm1);
     97:     VecCopy(vp1,v);

     99:     /*
    100:        SOLVE  Ax = b
    101:      */
    102:     /* Apply last two Given's (Gl-1 and Gl) rotations to (beta,alpha,Gamma) */
    103:     if (ksp->its > 2) {
    104:       theta =  sl1*beta;
    105:       eptmp = -cl1*beta;
    106:     }
    107:     if (ksp->its > 1) {
    108:       ep     = -cl*eptmp + sl*alpha;
    109:       deltmp = -sl*eptmp - cl*alpha;
    110:     }
    111:     if (PetscAbsReal(Gamma) > PetscAbsScalar(deltmp)) {
    112:       ta = -deltmp / Gamma;
    113:       s  = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
    114:       c  = s*ta;
    115:     } else {
    116:       ta = -Gamma/deltmp;
    117:       c  = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
    118:       s  = c*ta;
    119:     }

    121:     delta = -c*deltmp + s*Gamma;
    122:     tau_n = -c*tau_n1; tau_n1 = -s*tau_n1;
    123:     VecCopy(vm1,pvec);
    124:     VecAXPY(pvec,-theta,pvec2);
    125:     VecAXPY(pvec,-ep,pvec1);
    126:     VecScale(pvec,1.0/delta);
    127:     VecAXPY(x,tau_n,pvec);
    128:     cl1   = cl; sl1 = sl; cl = c; sl = s;

    130:     VecCopy(pvec1,pvec2);
    131:     VecCopy(pvec,pvec1);

    133:     /* Compute the upper bound on the residual norm r (See QMR paper p. 13) */
    134:     sprod = sprod*PetscAbsScalar(s);
    135:     rnorm = rnorm0 * PetscSqrtReal((PetscReal)ksp->its+2.0) * PetscRealPart(sprod);
    136:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
    137:     else ksp->rnorm = 0;
    138:     (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
    139:     if (ksp->its >= ksp->max_it) {
    140:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
    141:       break;
    142:     }
    143:   }
    144:   KSPMonitor(ksp,ksp->its,ksp->rnorm);
    145:   KSPUnwindPreconditioner(ksp,x,vtmp);
    146:   return(0);
    147: }

    149: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
    150: {

    154:   if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no symmetric preconditioning for KSPTCQMR");
    155:   KSPSetWorkVecs(ksp,TCQMR_VECS);
    156:   return(0);
    157: }

    159: /*MC
    160:      KSPTCQMR - A variant of QMR (quasi minimal residual) developed by Tony Chan

    162:    Options Database Keys:
    163: .   see KSPSolve()

    165:    Level: beginner

    167:   Notes:
    168:     Supports either left or right preconditioning, but not symmetric

    170:           The "residual norm" computed in this algorithm is actually just an upper bound on the actual residual norm.
    171:           That is for left preconditioning it is a bound on the preconditioned residual and for right preconditioning
    172:           it is a bound on the true residual.

    174:   References:
    175: .  1. - Tony F. Chan, Lisette de Pillis, and Henk van der Vorst, Transpose free formulations of Lanczos type methods for nonsymmetric linear systems,
    176:   Numerical Algorithms, Volume 17, 1998.

    178: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPTFQMR

    180: M*/

    182: PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
    183: {

    187:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
    188:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
    189:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

    191:   ksp->data                = (void*)0;
    192:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    193:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    194:   ksp->ops->setup          = KSPSetUp_TCQMR;
    195:   ksp->ops->solve          = KSPSolve_TCQMR;
    196:   ksp->ops->destroy        = KSPDestroyDefault;
    197:   ksp->ops->setfromoptions = NULL;
    198:   ksp->ops->view           = NULL;
    199:   return(0);
    200: }