Actual source code: tcqmr.c

      1: #define PETSCKSP_DLL

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

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

     10:  #include src/ksp/ksp/kspimpl.h
     11:  #include src/ksp/ksp/impls/tcqmr/tcqmrp.h

     15: static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
     16: {
     17:   PetscReal      rnorm0,rnorm,dp1,Gamma;
     18:   PetscScalar    theta,ep,cl1,sl1,cl,sl,sprod,tau_n1,f;
     19:   PetscScalar    deltmp,rho,beta,eptmp,ta,s,c,tau_n,delta;
     20:   PetscScalar    dp11,dp2,rhom1,alpha,tmp,zero = 0.0;

     24:   ksp->its = 0;

     26:   KSPInitialResidual(ksp,x,u,v,r,b);
     27:   VecNorm(r,NORM_2,&rnorm0);         /*  rnorm0 = ||r|| */

     29:   (*ksp->converged)(ksp,0,rnorm0,&ksp->reason,ksp->cnvP);
     30:   if (ksp->reason) return(0);

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

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

     62:     KSP_PCApplyBAorAB(ksp,u,y,vtmp); /* y = A*u */
     63:     VecDot(v0,y,&dp11);
     64:     VecDot(v0,u,&dp2);
     65:     alpha  = dp11 / dp2;                          /* alpha = v0'*y/v0'*u */
     66:     deltmp = alpha;
     67:     VecCopy(y,z);
     68:     tmp    = -alpha;
     69:     VecAXPY(z,tmp,u); /* z = y - alpha u */
     70:     VecDot(v0,u,&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:     tmp    = -2.0*beta;VecAXPY(utmp,tmp,p);
     77:     KSP_PCApplyBAorAB(ksp,utmp,up1,vtmp);
     78:     tmp    = -alpha; VecAXPY(up1,tmp,utmp);
     79:     tmp    = f*beta*beta; VecAXPY(up1,tmp,um1);
     80:     VecNorm(up1,NORM_2,&dp1);
     81:     f      = 1.0 / dp1;
     82:     VecScale(up1,f);
     83:     tmp    = -beta;
     84:     VecAYPX(p,tmp,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:     tmp    = -alpha; VecAXPY(vp1,tmp,v);
     92:     tmp    = -beta; VecAXPY(vp1,tmp,vm1);
     93:     VecNorm(vp1,NORM_2,&Gamma);
     94:     tmp    = 1.0/Gamma; VecScale(vp1,tmp);
     95:     VecCopy(v,vm1);
     96:     VecCopy(vp1,v);

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

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

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

    132:     /* Compute the upper bound on the residual norm r (See QMR paper p. 13) */
    133:     sprod = sprod*PetscAbsScalar(s);
    134: #if defined(PETSC_USE_COMPLEX)
    135:     rnorm = rnorm0 * sqrt((double)ksp->its+2.0) * PetscRealPart(sprod);
    136: #else
    137:     rnorm = rnorm0 * sqrt((double)ksp->its+2.0) * sprod;
    138: #endif
    139:     if (ksp->its >= ksp->max_it) {ksp->reason = KSP_DIVERGED_ITS; break;}
    140:     (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
    141:   }
    142:   KSPMonitor(ksp,ksp->its,rnorm);
    143:   KSPUnwindPreconditioner(ksp,x,vtmp);

    145:   return(0);
    146: }

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

    155:   if (ksp->pc_side == PC_SYMMETRIC){
    156:     SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPTCQMR");
    157:   }
    158:   KSPDefaultGetWork(ksp,TCQMR_VECS);
    159:   return(0);
    160: }

    162: /*MC
    163:      KSPRTCQMR - A variant of QMR (quasi minimal residual) developed by Tony Chan

    165:    Options Database Keys:
    166: .   see KSPSolve()

    168:    Level: beginner

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

    172: M*/

    177: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_TCQMR(KSP ksp)
    178: {
    180:   ksp->data                = (void*)0;
    181:   ksp->pc_side             = PC_LEFT;
    182:   ksp->ops->buildsolution  = KSPDefaultBuildSolution;
    183:   ksp->ops->buildresidual  = KSPDefaultBuildResidual;
    184:   ksp->ops->setup          = KSPSetUp_TCQMR;
    185:   ksp->ops->solve          = KSPSolve_TCQMR;
    186:   ksp->ops->destroy        = KSPDefaultDestroy;
    187:   ksp->ops->setfromoptions = 0;
    188:   ksp->ops->view           = 0;
    189:   return(0);
    190: }