Actual source code: symmlq.c

      1: #include <petsc/private/kspimpl.h>

      3: typedef struct {
      4:   PetscReal haptol;
      5: } KSP_SYMMLQ;

      7: static PetscErrorCode KSPSetUp_SYMMLQ(KSP ksp)
      8: {
      9:   PetscFunctionBegin;
     10:   PetscCall(KSPSetWorkVecs(ksp, 9));
     11:   PetscFunctionReturn(PETSC_SUCCESS);
     12: }

     14: static PetscErrorCode KSPSolve_SYMMLQ(KSP ksp)
     15: {
     16:   PetscInt    i;
     17:   PetscScalar alpha, beta, ibeta, betaold, beta1, ceta = 0, ceta_oold = 0.0, ceta_old = 0.0, ceta_bar;
     18:   PetscScalar c = 1.0, cold = 1.0, s = 0.0, sold = 0.0, coold, soold, rho0, rho1, rho2, rho3;
     19:   PetscScalar dp = 0.0;
     20:   PetscReal   np = 0.0, s_prod;
     21:   Vec         X, B, R, Z, U, V, W, UOLD, VOLD, Wbar;
     22:   Mat         Amat, Pmat;
     23:   KSP_SYMMLQ *symmlq = (KSP_SYMMLQ *)ksp->data;
     24:   PetscBool   diagonalscale;

     26:   PetscFunctionBegin;
     27:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
     28:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

     30:   X    = ksp->vec_sol;
     31:   B    = ksp->vec_rhs;
     32:   R    = ksp->work[0];
     33:   Z    = ksp->work[1];
     34:   U    = ksp->work[2];
     35:   V    = ksp->work[3];
     36:   W    = ksp->work[4];
     37:   UOLD = ksp->work[5];
     38:   VOLD = ksp->work[6];
     39:   Wbar = ksp->work[7];

     41:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

     43:   ksp->its = 0;

     45:   PetscCall(VecSet(UOLD, 0.0));   /* u_old <- zeros;  */
     46:   PetscCall(VecCopy(UOLD, VOLD)); /* v_old <- u_old;  */
     47:   PetscCall(VecCopy(UOLD, W));    /* w     <- u_old;  */
     48:   PetscCall(VecCopy(UOLD, Wbar)); /* w_bar <- u_old;  */
     49:   if (!ksp->guess_zero) {
     50:     PetscCall(KSP_MatMult(ksp, Amat, X, R)); /*     r <- b - A*x */
     51:     PetscCall(VecAYPX(R, -1.0, B));
     52:   } else {
     53:     PetscCall(VecCopy(B, R)); /*     r <- b (x is 0) */
     54:   }

     56:   PetscCall(KSP_PCApply(ksp, R, Z)); /* z  <- B*r       */
     57:   PetscCall(VecDot(R, Z, &dp));      /* dp = r'*z;      */
     58:   KSPCheckDot(ksp, dp);
     59:   if (PetscAbsScalar(dp) < symmlq->haptol) {
     60:     PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol));
     61:     ksp->rnorm  = 0.0;                           /* what should we really put here? */
     62:     ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN; /* bugfix proposed by Lourens (lourens.vanzanen@shell.com) */
     63:     PetscFunctionReturn(PETSC_SUCCESS);
     64:   }

     66: #if !defined(PETSC_USE_COMPLEX)
     67:   if (dp < 0.0) {
     68:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
     69:     PetscFunctionReturn(PETSC_SUCCESS);
     70:   }
     71: #endif
     72:   dp     = PetscSqrtScalar(dp);
     73:   beta   = dp; /*  beta <- sqrt(r'*z)  */
     74:   beta1  = beta;
     75:   s_prod = PetscAbsScalar(beta1);

     77:   PetscCall(VecCopy(R, V)); /* v <- r; */
     78:   PetscCall(VecCopy(Z, U)); /* u <- z; */
     79:   ibeta = 1.0 / beta;
     80:   PetscCall(VecScale(V, ibeta)); /* v <- ibeta*v; */
     81:   PetscCall(VecScale(U, ibeta)); /* u <- ibeta*u; */
     82:   PetscCall(VecCopy(U, Wbar));   /* w_bar <- u;   */
     83:   if (ksp->normtype != KSP_NORM_NONE) {
     84:     PetscCall(VecNorm(Z, NORM_2, &np)); /*   np <- ||z||        */
     85:     KSPCheckNorm(ksp, np);
     86:   }
     87:   PetscCall(KSPLogResidualHistory(ksp, np));
     88:   PetscCall(KSPMonitor(ksp, 0, np));
     89:   ksp->rnorm = np;
     90:   PetscCall((*ksp->converged)(ksp, 0, np, &ksp->reason, ksp->cnvP)); /* test for convergence */
     91:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

     93:   i    = 0;
     94:   ceta = 0.;
     95:   do {
     96:     ksp->its = i + 1;

     98:     /*    Update    */
     99:     if (ksp->its > 1) {
    100:       PetscCall(VecCopy(V, VOLD)); /* v_old <- v; */
    101:       PetscCall(VecCopy(U, UOLD)); /* u_old <- u; */

    103:       PetscCall(VecCopy(R, V));
    104:       PetscCall(VecScale(V, 1.0 / beta)); /* v <- ibeta*r; */
    105:       PetscCall(VecCopy(Z, U));
    106:       PetscCall(VecScale(U, 1.0 / beta)); /* u <- ibeta*z; */

    108:       PetscCall(VecCopy(Wbar, W));
    109:       PetscCall(VecScale(W, c));
    110:       PetscCall(VecAXPY(W, s, U)); /* w  <- c*w_bar + s*u;    (w_k) */
    111:       PetscCall(VecScale(Wbar, -s));
    112:       PetscCall(VecAXPY(Wbar, c, U)); /* w_bar <- -s*w_bar + c*u; (w_bar_(k+1)) */
    113:       PetscCall(VecAXPY(X, ceta, W)); /* x <- x + ceta * w;       (xL_k)  */

    115:       ceta_oold = ceta_old;
    116:       ceta_old  = ceta;
    117:     }

    119:     /*   Lanczos  */
    120:     PetscCall(KSP_MatMult(ksp, Amat, U, R)); /*  r     <- Amat*u; */
    121:     PetscCall(VecDot(U, R, &alpha));         /*  alpha <- u'*r;   */
    122:     PetscCall(KSP_PCApply(ksp, R, Z));       /*      z <- B*r;    */

    124:     PetscCall(VecAXPY(R, -alpha, V));   /*  r <- r - alpha* v;  */
    125:     PetscCall(VecAXPY(Z, -alpha, U));   /*  z <- z - alpha* u;  */
    126:     PetscCall(VecAXPY(R, -beta, VOLD)); /*  r <- r - beta * v_old; */
    127:     PetscCall(VecAXPY(Z, -beta, UOLD)); /*  z <- z - beta * u_old; */
    128:     betaold = beta;                     /* beta_k                  */
    129:     PetscCall(VecDot(R, Z, &dp));       /* dp <- r'*z;             */
    130:     KSPCheckDot(ksp, dp);
    131:     if (PetscAbsScalar(dp) < symmlq->haptol) {
    132:       PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol));
    133:       dp = 0.0;
    134:     }

    136: #if !defined(PETSC_USE_COMPLEX)
    137:     if (dp < 0.0) {
    138:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    139:       break;
    140:     }
    141: #endif
    142:     beta = PetscSqrtScalar(dp); /*  beta = sqrt(dp); */

    144:     /*    QR factorization    */
    145:     coold = cold;
    146:     cold  = c;
    147:     soold = sold;
    148:     sold  = s;
    149:     rho0  = cold * alpha - coold * sold * betaold;      /* gamma_bar */
    150:     rho1  = PetscSqrtScalar(rho0 * rho0 + beta * beta); /* gamma     */
    151:     rho2  = sold * alpha + coold * cold * betaold;      /* delta     */
    152:     rho3  = soold * betaold;                            /* epsilon   */

    154:     /* Givens rotation: [c -s; s c] (different from the Reference!) */
    155:     c = rho0 / rho1;
    156:     s = beta / rho1;

    158:     if (ksp->its == 1) ceta = beta1 / rho1;
    159:     else ceta = -(rho2 * ceta_old + rho3 * ceta_oold) / rho1;

    161:     s_prod = s_prod * PetscAbsScalar(s);
    162:     if (c == 0.0) np = s_prod * 1.e16;
    163:     else np = s_prod / PetscAbsScalar(c); /* residual norm for xc_k (CGNORM) */

    165:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
    166:     else ksp->rnorm = 0.0;
    167:     PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
    168:     PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
    169:     PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
    170:     if (ksp->reason) break;
    171:     i++;
    172:   } while (i < ksp->max_it);

    174:   /* move to the CG point: xc_(k+1) */
    175:   if (c == 0.0) ceta_bar = ceta * 1.e15;
    176:   else ceta_bar = ceta / c;

    178:   PetscCall(VecAXPY(X, ceta_bar, Wbar)); /* x <- x + ceta_bar*w_bar */

    180:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: /*MC
    185:    KSPSYMMLQ - Implements the SYMMLQ method {cite}`paige.saunders:solution`.

    187:    Level: beginner

    189:    Notes:
    190:    The operator and the preconditioner must be symmetric for this method.

    192:    The preconditioner must be POSITIVE-DEFINITE.

    194:    Supports only left preconditioning.

    196: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`
    197: M*/
    198: PETSC_EXTERN PetscErrorCode KSPCreate_SYMMLQ(KSP ksp)
    199: {
    200:   KSP_SYMMLQ *symmlq;

    202:   PetscFunctionBegin;
    203:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
    204:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

    206:   PetscCall(PetscNew(&symmlq));
    207:   symmlq->haptol = 1.e-18;
    208:   ksp->data      = (void *)symmlq;

    210:   /*
    211:        Sets the functions that are associated with this data structure
    212:        (in C++ this is the same as defining virtual functions)
    213:   */
    214:   ksp->ops->setup          = KSPSetUp_SYMMLQ;
    215:   ksp->ops->solve          = KSPSolve_SYMMLQ;
    216:   ksp->ops->destroy        = KSPDestroyDefault;
    217:   ksp->ops->setfromoptions = NULL;
    218:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    219:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    220:   PetscFunctionReturn(PETSC_SUCCESS);
    221: }