Actual source code: tcqmr.c

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

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

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

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

     17:   PetscFunctionBegin;
     18:   ksp->its = 0;

     20:   PetscCall(KSPInitialResidual(ksp, x, u, v, r, b));
     21:   PetscCall(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:   PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
     26:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

     28:   PetscCall(VecSet(um1, 0.0));
     29:   PetscCall(VecCopy(r, u));
     30:   rnorm = rnorm0;
     31:   tmp   = 1.0 / rnorm;
     32:   PetscCall(VecScale(u, tmp));
     33:   PetscCall(VecSet(vm1, 0.0));
     34:   PetscCall(VecCopy(u, v));
     35:   PetscCall(VecCopy(u, v0));
     36:   PetscCall(VecSet(pvec1, 0.0));
     37:   PetscCall(VecSet(pvec2, 0.0));
     38:   PetscCall(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:   PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
     57:   while (!ksp->reason) {
     58:     PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
     59:     ksp->its++;

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

    133:     PetscCall(VecCopy(pvec1, pvec2));
    134:     PetscCall(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:     PetscCall((*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:   PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
    148:   PetscCall(KSPUnwindPreconditioner(ksp, x, vtmp));
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
    153: {
    154:   PetscFunctionBegin;
    155:   PetscCall(KSPSetWorkVecs(ksp, TCQMR_VECS));
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: /*MC
    160:    KSPTCQMR - A variant of a tranpose-free QMR (quasi minimal residual) {cite}`chan1998transpose` algorithm

    162:    Level: beginner

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

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

    171:    This algorithm is related to the standard transpose-free QMR implemented by `KSPTFQMR`.

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

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

    183:   ksp->data                = NULL;
    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:   PetscFunctionReturn(PETSC_SUCCESS);
    192: }