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;

 18:   ksp->its = 0;

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

 28:   VecSet(um1, 0.0);
 29:   VecCopy(r, u);
 30:   rnorm = rnorm0;
 31:   tmp   = 1.0 / rnorm;
 32:   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:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
 55:   else ksp->rnorm = 0;
 56:   (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
 57:   while (!ksp->reason) {
 58:     KSPMonitor(ksp, ksp->its, ksp->rnorm);
 59:     ksp->its++;

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

133:     VecCopy(pvec1, pvec2);
134:     VecCopy(pvec, pvec1);

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

152: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
153: {
154:   KSPSetWorkVecs(ksp, TCQMR_VECS);
155:   return 0;
156: }

158: /*MC
159:      KSPTCQMR - A variant of QMR (quasi minimal residual) [1]

161:    Level: beginner

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

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

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

174: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPTFQMR`
175: M*/

177: PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
178: {
179:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
180:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2);
181:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);

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