Actual source code: tcqmr.c

petsc-master 2019-08-19
Report Typos and Errors

  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);

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

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

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

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

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

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

127:     VecCopy(pvec1,pvec2);
128:     VecCopy(pvec,pvec1);

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

144: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
145: {

149:   if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no symmetric preconditioning for KSPTCQMR");
150:   KSPSetWorkVecs(ksp,TCQMR_VECS);
151:   return(0);
152: }

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

157:    Options Database Keys:
158: .   see KSPSolve()

160:    Level: beginner

162:   Notes:
163:     Supports either left or right preconditioning, but not symmetric

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

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

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

175: M*/

177: PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
178: {

182:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
183:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);

185:   ksp->data                = (void*)0;
186:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
187:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
188:   ksp->ops->setup          = KSPSetUp_TCQMR;
189:   ksp->ops->solve          = KSPSolve_TCQMR;
190:   ksp->ops->destroy        = KSPDestroyDefault;
191:   ksp->ops->setfromoptions = 0;
192:   ksp->ops->view           = 0;
193:   return(0);
194: }