Actual source code: cg.c

      1: /*
      2:     This file implements the conjugate gradient method in PETSc as part of
      3:     KSP. You can use this as a starting point for implementing your own
      4:     Krylov method that is not provided with PETSc.

      6:     The following basic routines are required for each Krylov method.
      7:         KSPCreate_XXX()          - Creates the Krylov context
      8:         KSPSetFromOptions_XXX()  - Sets runtime options
      9:         KSPSolve_XXX()           - Runs the Krylov method
     10:         KSPDestroy_XXX()         - Destroys the Krylov context, freeing all
     11:                                    memory it needed
     12:     Here the "_XXX" denotes a particular implementation, in this case
     13:     we use _CG (e.g. KSPCreate_CG, KSPDestroy_CG). These routines
     14:     are actually called via the common user interface routines
     15:     KSPSetType(), KSPSetFromOptions(), KSPSolve(), and KSPDestroy() so the
     16:     application code interface remains identical for all preconditioners.

     18:     Other basic routines for the KSP objects include
     19:         KSPSetUp_XXX()
     20:         KSPView_XXX()            - Prints details of solver being used.

     22:     Detailed Notes:
     23:     By default, this code implements the CG (Conjugate Gradient) method,
     24:     which is valid for real symmetric (and complex Hermitian) positive
     25:     definite matrices. Note that for the complex Hermitian case, the
     26:     VecDot() arguments within the code MUST remain in the order given
     27:     for correct computation of inner products.

     29:     Reference: Hestenes and Steifel, 1952.

     31:     By switching to the indefinite vector inner product, VecTDot(), the
     32:     same code is used for the complex symmetric case as well.  The user
     33:     must call KSPCGSetType(ksp,KSP_CG_SYMMETRIC) or use the option
     34:     -ksp_cg_type symmetric to invoke this variant for the complex case.
     35:     Note, however, that the complex symmetric code is NOT valid for
     36:     all such matrices ... and thus we don't recommend using this method.
     37: */
     38: /*
     39:     cgimpl.h defines the simple data structured used to store information
     40:     related to the type of matrix (e.g. complex symmetric) being solved and
     41:     data used during the optional Lanczos process used to compute eigenvalues
     42: */
     43: #include <../src/ksp/ksp/impls/cg/cgimpl.h>
     44: extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP, PetscReal *, PetscReal *);
     45: extern PetscErrorCode KSPComputeEigenvalues_CG(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);

     47: static PetscErrorCode KSPCGSetObjectiveTarget_CG(KSP ksp, PetscReal obj_min)
     48: {
     49:   KSP_CG *cg = (KSP_CG *)ksp->data;

     51:   PetscFunctionBegin;
     52:   cg->obj_min = obj_min;
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

     56: static PetscErrorCode KSPCGSetRadius_CG(KSP ksp, PetscReal radius)
     57: {
     58:   KSP_CG *cg = (KSP_CG *)ksp->data;

     60:   PetscFunctionBegin;
     61:   cg->radius = radius;
     62:   PetscFunctionReturn(PETSC_SUCCESS);
     63: }

     65: static PetscErrorCode KSPCGGetObjFcn_CG(KSP ksp, PetscReal *obj)
     66: {
     67:   KSP_CG *cg = (KSP_CG *)ksp->data;

     69:   PetscFunctionBegin;
     70:   *obj = cg->obj;
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /*
     75:      KSPSetUp_CG - Sets up the workspace needed by the CG method.

     77:       This is called once, usually automatically by KSPSolve() or KSPSetUp()
     78:      but can be called directly by KSPSetUp()
     79: */
     80: static PetscErrorCode KSPSetUp_CG(KSP ksp)
     81: {
     82:   KSP_CG  *cgP   = (KSP_CG *)ksp->data;
     83:   PetscInt maxit = ksp->max_it, nwork = 3;

     85:   PetscFunctionBegin;
     86:   /* get work vectors needed by CG */
     87:   if (cgP->singlereduction) nwork += 2;
     88:   PetscCall(KSPSetWorkVecs(ksp, nwork));

     90:   /*
     91:      If user requested computations of eigenvalues then allocate
     92:      work space needed
     93:   */
     94:   if (ksp->calc_sings) {
     95:     PetscCall(PetscFree4(cgP->e, cgP->d, cgP->ee, cgP->dd));
     96:     PetscCall(PetscMalloc4(maxit + 1, &cgP->e, maxit, &cgP->d, maxit, &cgP->ee, maxit, &cgP->dd));

     98:     ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
     99:     ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_CG;
    100:   }
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: /*
    105:      A macro used in the following KSPSolve_CG and KSPSolve_CG_SingleReduction routines
    106: */
    107: #define VecXDot(x, y, a) (cg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))

    109: /*
    110:      KSPSolve_CG - This routine actually applies the conjugate gradient method

    112:      Note : this routine can be replaced with another one (see below) which implements
    113:             another variant of CG.

    115:    Input Parameter:
    116: .     ksp - the Krylov space object that was set to use conjugate gradient, by, for
    117:             example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
    118: */
    119: static PetscErrorCode KSPSolve_CG(KSP ksp)
    120: {
    121:   PetscInt    i, stored_max_it, eigs;
    122:   PetscScalar dpi = 0.0, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL, dpiold;
    123:   PetscReal   dp = 0.0;
    124:   PetscReal   r2, norm_p, norm_d, dMp;
    125:   Vec         X, B, Z, R, P, W;
    126:   KSP_CG     *cg;
    127:   Mat         Amat, Pmat;
    128:   PetscBool   diagonalscale, testobj;

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

    134:   cg            = (KSP_CG *)ksp->data;
    135:   eigs          = ksp->calc_sings;
    136:   stored_max_it = ksp->max_it;
    137:   X             = ksp->vec_sol;
    138:   B             = ksp->vec_rhs;
    139:   R             = ksp->work[0];
    140:   Z             = ksp->work[1];
    141:   P             = ksp->work[2];
    142:   W             = Z;
    143:   r2            = PetscSqr(cg->radius);

    145:   if (eigs) {
    146:     e    = cg->e;
    147:     d    = cg->d;
    148:     e[0] = 0.0;
    149:   }
    150:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

    152:   ksp->its = 0;
    153:   if (!ksp->guess_zero) {
    154:     PetscCall(KSP_MatMult(ksp, Amat, X, R)); /*    r <- b - Ax                       */

    156:     PetscCall(VecAYPX(R, -1.0, B));
    157:     if (cg->radius) { /* XXX direction? */
    158:       PetscCall(VecNorm(X, NORM_2, &norm_d));
    159:       norm_d *= norm_d;
    160:     }
    161:   } else {
    162:     PetscCall(VecCopy(B, R)); /*    r <- b (x is 0)                   */
    163:     norm_d = 0.0;
    164:   }
    165:   /* This may be true only on a subset of MPI ranks; setting it here so it will be detected by the first norm computation below */
    166:   PetscCall(VecFlag(R, ksp->reason == KSP_DIVERGED_PC_FAILED));

    168:   switch (ksp->normtype) {
    169:   case KSP_NORM_PRECONDITIONED:
    170:     PetscCall(KSP_PCApply(ksp, R, Z));  /*    z <- Br                           */
    171:     PetscCall(VecNorm(Z, NORM_2, &dp)); /*    dp <- z'*z = e'*A'*B'*B*A*e       */
    172:     KSPCheckNorm(ksp, dp);
    173:     break;
    174:   case KSP_NORM_UNPRECONDITIONED:
    175:     PetscCall(VecNorm(R, NORM_2, &dp)); /*    dp <- r'*r = e'*A'*A*e            */
    176:     KSPCheckNorm(ksp, dp);
    177:     break;
    178:   case KSP_NORM_NATURAL:
    179:     PetscCall(KSP_PCApply(ksp, R, Z)); /*    z <- Br                           */
    180:     PetscCall(VecXDot(Z, R, &beta));   /*    beta <- z'*r                      */
    181:     KSPCheckDot(ksp, beta);
    182:     dp = PetscSqrtReal(PetscAbsScalar(beta)); /*    dp <- r'*z = r'*B*r = e'*A'*B*A*e */
    183:     break;
    184:   case KSP_NORM_NONE:
    185:     dp = 0.0;
    186:     break;
    187:   default:
    188:     SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
    189:   }

    191:   /* Initialize objective function
    192:      obj = 1/2 x^T A x - x^T b */
    193:   testobj = (PetscBool)(cg->obj_min < 0.0);
    194:   PetscCall(VecXDot(R, X, &a));
    195:   cg->obj = 0.5 * PetscRealPart(a);
    196:   PetscCall(VecXDot(B, X, &a));
    197:   cg->obj -= 0.5 * PetscRealPart(a);

    199:   if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", ksp->its, (double)cg->obj));
    200:   PetscCall(KSPLogResidualHistory(ksp, dp));
    201:   PetscCall(KSPMonitor(ksp, ksp->its, dp));
    202:   ksp->rnorm = dp;

    204:   PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */

    206:   if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
    207:     PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
    208:     ksp->reason = KSP_CONVERGED_ATOL;
    209:   }

    211:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

    213:   if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /*     z <- Br                           */
    214:   if (ksp->normtype != KSP_NORM_NATURAL) {
    215:     PetscCall(VecXDot(Z, R, &beta)); /*     beta <- z'*r                      */
    216:     KSPCheckDot(ksp, beta);
    217:   }

    219:   i = 0;
    220:   do {
    221:     ksp->its = i + 1;
    222:     if (beta == 0.0) {
    223:       ksp->reason = KSP_CONVERGED_ATOL;
    224:       PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
    225:       break;
    226: #if !defined(PETSC_USE_COMPLEX)
    227:     } else if ((i > 0) && (beta * betaold < 0.0)) {
    228:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner, beta %g, betaold %g", (double)PetscRealPart(beta), (double)PetscRealPart(betaold));
    229:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    230:       PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
    231:       break;
    232: #endif
    233:     }
    234:     if (!i) {
    235:       PetscCall(VecCopy(Z, P)); /*     p <- z                           */
    236:       if (cg->radius) {
    237:         PetscCall(VecNorm(P, NORM_2, &norm_p));
    238:         norm_p *= norm_p;
    239:         dMp = 0.0;
    240:         if (!ksp->guess_zero) PetscCall(VecDotRealPart(X, P, &dMp));
    241:       }
    242:       b = 0.0;
    243:     } else {
    244:       b = beta / betaold;
    245:       if (eigs) {
    246:         PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
    247:         e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
    248:       }
    249:       PetscCall(VecAYPX(P, b, Z)); /*     p <- z + b* p                    */
    250:       if (cg->radius) {
    251:         PetscCall(VecDotRealPart(X, P, &dMp));
    252:         PetscCall(VecNorm(P, NORM_2, &norm_p));
    253:         norm_p *= norm_p;
    254:       }
    255:     }
    256:     dpiold = dpi;
    257:     PetscCall(KSP_MatMult(ksp, Amat, P, W)); /*     w <- Ap                          */
    258:     PetscCall(VecXDot(P, W, &dpi));          /*     dpi <- p'w                       */
    259:     KSPCheckDot(ksp, dpi);
    260:     betaold = beta;

    262:     if ((dpi == 0.0) || ((i > 0) && ((PetscSign(PetscRealPart(dpi)) * PetscSign(PetscRealPart(dpiold))) < 0.0))) {
    263:       if (cg->radius) {
    264:         a = 0.0;
    265:         if (i == 0) {
    266:           if (norm_p > 0.0) {
    267:             a = PetscSqrtReal(r2 / norm_p);
    268:           } else {
    269:             PetscCall(VecNorm(R, NORM_2, &dp));
    270:             a = cg->radius > dp ? 1.0 : cg->radius / dp;
    271:           }
    272:         } else if (norm_p > 0.0) {
    273:           a = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
    274:         }
    275:         PetscCall(VecAXPY(X, a, P)); /*     x <- x + ap                      */
    276:         cg->obj += PetscRealPart(a * (0.5 * a * dpi - betaold));
    277:       }
    278:       if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " N obj %g\n", i + 1, (double)cg->obj));
    279:       if (ksp->converged_neg_curve) {
    280:         PetscCall(PetscInfo(ksp, "converged due to negative curvature: %g\n", (double)(PetscRealPart(dpi))));
    281:         ksp->reason = KSP_CONVERGED_NEG_CURVE;
    282:       } else {
    283:         PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix, dpi %g, dpiold %g", (double)PetscRealPart(dpi), (double)PetscRealPart(dpiold));
    284:         ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
    285:         PetscCall(PetscInfo(ksp, "diverging due to indefinite matrix\n"));
    286:       }
    287:       break;
    288:     }
    289:     a = beta / dpi; /*     a = beta/p'w                     */
    290:     if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
    291:     if (cg->radius) { /* Steihaugh-Toint */
    292:       PetscReal norm_dp1 = norm_d + PetscRealPart(a) * (2.0 * dMp + PetscRealPart(a) * norm_p);
    293:       if (norm_dp1 > r2) {
    294:         ksp->reason = KSP_CONVERGED_STEP_LENGTH;
    295:         PetscCall(PetscInfo(ksp, "converged to the trust region radius %g\n", (double)cg->radius));
    296:         if (norm_p > 0.0) {
    297:           dp = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
    298:           PetscCall(VecAXPY(X, dp, P)); /*     x <- x + ap                      */
    299:           cg->obj += PetscRealPart(dp * (0.5 * dp * dpi - beta));
    300:         }
    301:         if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " R obj %g\n", i + 1, (double)cg->obj));
    302:         break;
    303:       }
    304:     }
    305:     PetscCall(VecAXPY(X, a, P));  /*     x <- x + ap                      */
    306:     PetscCall(VecAXPY(R, -a, W)); /*     r <- r - aw                      */
    307:     if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
    308:       PetscCall(KSP_PCApply(ksp, R, Z));  /*     z <- Br                          */
    309:       PetscCall(VecNorm(Z, NORM_2, &dp)); /*     dp <- z'*z                       */
    310:       KSPCheckNorm(ksp, dp);
    311:     } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
    312:       PetscCall(VecNorm(R, NORM_2, &dp)); /*     dp <- r'*r                       */
    313:       KSPCheckNorm(ksp, dp);
    314:     } else if (ksp->normtype == KSP_NORM_NATURAL) {
    315:       PetscCall(KSP_PCApply(ksp, R, Z)); /*     z <- Br                          */
    316:       PetscCall(VecXDot(Z, R, &beta));   /*     beta <- r'*z                     */
    317:       KSPCheckDot(ksp, beta);
    318:       dp = PetscSqrtReal(PetscAbsScalar(beta));
    319:     } else {
    320:       dp = 0.0;
    321:     }
    322:     cg->obj -= PetscRealPart(0.5 * a * betaold);
    323:     if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", i + 1, (double)cg->obj));

    325:     ksp->rnorm = dp;
    326:     PetscCall(KSPLogResidualHistory(ksp, dp));
    327:     PetscCall(KSPMonitor(ksp, i + 1, dp));
    328:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));

    330:     if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
    331:       PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
    332:       ksp->reason = KSP_CONVERGED_ATOL;
    333:     }

    335:     if (ksp->reason) break;

    337:     if (cg->radius) {
    338:       PetscCall(VecNorm(X, NORM_2, &norm_d));
    339:       norm_d *= norm_d;
    340:     }

    342:     if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) PetscCall(KSP_PCApply(ksp, R, Z)); /*     z <- Br                          */
    343:     if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
    344:       PetscCall(VecXDot(Z, R, &beta)); /*     beta <- z'*r                     */
    345:       KSPCheckDot(ksp, beta);
    346:     }

    348:     i++;
    349:   } while (i < ksp->max_it);
    350:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: /*
    355:        KSPSolve_CG_SingleReduction

    357:        This variant of CG is identical in exact arithmetic to the standard algorithm,
    358:        but is rearranged to use only a single reduction stage per iteration, using additional
    359:        intermediate vectors.

    361:        See KSPCGUseSingleReduction_CG()

    363: */
    364: static PetscErrorCode KSPSolve_CG_SingleReduction(KSP ksp)
    365: {
    366:   PetscInt    i, stored_max_it, eigs;
    367:   PetscScalar dpi = 0.0, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL, delta, dpiold, tmp[2];
    368:   PetscReal   dp = 0.0;
    369:   Vec         X, B, Z, R, P, S, W, tmpvecs[2];
    370:   KSP_CG     *cg;
    371:   Mat         Amat, Pmat;
    372:   PetscBool   diagonalscale;

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

    378:   cg            = (KSP_CG *)ksp->data;
    379:   eigs          = ksp->calc_sings;
    380:   stored_max_it = ksp->max_it;
    381:   X             = ksp->vec_sol;
    382:   B             = ksp->vec_rhs;
    383:   R             = ksp->work[0];
    384:   Z             = ksp->work[1];
    385:   P             = ksp->work[2];
    386:   S             = ksp->work[3];
    387:   W             = ksp->work[4];

    389:   if (eigs) {
    390:     e    = cg->e;
    391:     d    = cg->d;
    392:     e[0] = 0.0;
    393:   }
    394:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

    396:   ksp->its = 0;
    397:   if (!ksp->guess_zero) {
    398:     PetscCall(KSP_MatMult(ksp, Amat, X, R)); /*    r <- b - Ax                       */
    399:     PetscCall(VecAYPX(R, -1.0, B));
    400:   } else {
    401:     PetscCall(VecCopy(B, R)); /*    r <- b (x is 0)                   */
    402:   }

    404:   switch (ksp->normtype) {
    405:   case KSP_NORM_PRECONDITIONED:
    406:     PetscCall(KSP_PCApply(ksp, R, Z));  /*    z <- Br                           */
    407:     PetscCall(VecNorm(Z, NORM_2, &dp)); /*    dp <- z'*z = e'*A'*B'*B*A'*e'     */
    408:     KSPCheckNorm(ksp, dp);
    409:     break;
    410:   case KSP_NORM_UNPRECONDITIONED:
    411:     PetscCall(VecNorm(R, NORM_2, &dp)); /*    dp <- r'*r = e'*A'*A*e            */
    412:     KSPCheckNorm(ksp, dp);
    413:     break;
    414:   case KSP_NORM_NATURAL:
    415:     PetscCall(KSP_PCApply(ksp, R, Z)); /*    z <- Br                           */
    416:     PetscCall(KSP_MatMult(ksp, Amat, Z, S));
    417:     PetscCall(VecXDot(Z, S, &delta)); /*    delta <- z'*A*z = r'*B*A*B*r      */
    418:     PetscCall(VecXDot(Z, R, &beta));  /*    beta <- z'*r                      */
    419:     KSPCheckDot(ksp, beta);
    420:     dp = PetscSqrtReal(PetscAbsScalar(beta)); /*    dp <- r'*z = r'*B*r = e'*A'*B*A*e */
    421:     break;
    422:   case KSP_NORM_NONE:
    423:     dp = 0.0;
    424:     break;
    425:   default:
    426:     SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
    427:   }
    428:   PetscCall(KSPLogResidualHistory(ksp, dp));
    429:   PetscCall(KSPMonitor(ksp, 0, dp));
    430:   ksp->rnorm = dp;

    432:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
    433:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

    435:   if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /*    z <- Br                           */
    436:   if (ksp->normtype != KSP_NORM_NATURAL) {
    437:     PetscCall(KSP_MatMult(ksp, Amat, Z, S));
    438:     PetscCall(VecXDot(Z, S, &delta)); /*    delta <- z'*A*z = r'*B*A*B*r      */
    439:     PetscCall(VecXDot(Z, R, &beta));  /*    beta <- z'*r                      */
    440:     KSPCheckDot(ksp, beta);
    441:   }

    443:   i = 0;
    444:   do {
    445:     ksp->its = i + 1;
    446:     if (beta == 0.0) {
    447:       ksp->reason = KSP_CONVERGED_ATOL;
    448:       PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
    449:       break;
    450: #if !defined(PETSC_USE_COMPLEX)
    451:     } else if ((i > 0) && (beta * betaold < 0.0)) {
    452:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner");
    453:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    454:       PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
    455:       break;
    456: #endif
    457:     }
    458:     if (!i) {
    459:       PetscCall(VecCopy(Z, P)); /*    p <- z                           */
    460:       b = 0.0;
    461:     } else {
    462:       b = beta / betaold;
    463:       if (eigs) {
    464:         PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
    465:         e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
    466:       }
    467:       PetscCall(VecAYPX(P, b, Z)); /*    p <- z + b* p                     */
    468:     }
    469:     dpiold = dpi;
    470:     if (!i) {
    471:       PetscCall(KSP_MatMult(ksp, Amat, P, W)); /*    w <- Ap                           */
    472:       PetscCall(VecXDot(P, W, &dpi));          /*    dpi <- p'w                        */
    473:     } else {
    474:       PetscCall(VecAYPX(W, beta / betaold, S));                 /*    w <- Ap                           */
    475:       dpi = delta - beta * beta * dpiold / (betaold * betaold); /*    dpi <- p'w                        */
    476:     }
    477:     betaold = beta;
    478:     KSPCheckDot(ksp, beta);

    480:     if ((dpi == 0.0) || ((i > 0) && (PetscRealPart(dpi * dpiold) <= 0.0))) {
    481:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix");
    482:       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
    483:       PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n"));
    484:       break;
    485:     }
    486:     a = beta / dpi; /*    a = beta/p'w                      */
    487:     if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
    488:     PetscCall(VecAXPY(X, a, P));  /*    x <- x + ap                       */
    489:     PetscCall(VecAXPY(R, -a, W)); /*    r <- r - aw                       */
    490:     if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
    491:       PetscCall(KSP_PCApply(ksp, R, Z)); /*    z <- Br                           */
    492:       PetscCall(KSP_MatMult(ksp, Amat, Z, S));
    493:       PetscCall(VecNorm(Z, NORM_2, &dp)); /*    dp <- z'*z                        */
    494:       KSPCheckNorm(ksp, dp);
    495:     } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
    496:       PetscCall(VecNorm(R, NORM_2, &dp)); /*    dp <- r'*r                        */
    497:       KSPCheckNorm(ksp, dp);
    498:     } else if (ksp->normtype == KSP_NORM_NATURAL) {
    499:       PetscCall(KSP_PCApply(ksp, R, Z)); /*    z <- Br                           */
    500:       tmpvecs[0] = S;
    501:       tmpvecs[1] = R;
    502:       PetscCall(KSP_MatMult(ksp, Amat, Z, S));
    503:       PetscCall(VecMDot(Z, 2, tmpvecs, tmp)); /*    delta <- z'*A*z = r'*B*A*B*r      */
    504:       delta = tmp[0];
    505:       beta  = tmp[1]; /*    beta <- z'*r                      */
    506:       KSPCheckDot(ksp, beta);
    507:       dp = PetscSqrtReal(PetscAbsScalar(beta)); /*    dp <- r'*z = r'*B*r = e'*A'*B*A*e */
    508:     } else {
    509:       dp = 0.0;
    510:     }
    511:     ksp->rnorm = dp;
    512:     PetscCall(KSPLogResidualHistory(ksp, dp));
    513:     PetscCall(KSPMonitor(ksp, i + 1, dp));
    514:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
    515:     if (ksp->reason) break;

    517:     if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) {
    518:       PetscCall(KSP_PCApply(ksp, R, Z)); /*    z <- Br                           */
    519:       PetscCall(KSP_MatMult(ksp, Amat, Z, S));
    520:     }
    521:     if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
    522:       tmpvecs[0] = S;
    523:       tmpvecs[1] = R;
    524:       PetscCall(VecMDot(Z, 2, tmpvecs, tmp));
    525:       delta = tmp[0];
    526:       beta  = tmp[1];         /*    delta <- z'*A*z = r'*B'*A*B*r     */
    527:       KSPCheckDot(ksp, beta); /*    beta <- z'*r                      */
    528:     }

    530:     i++;
    531:   } while (i < ksp->max_it);
    532:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    533:   PetscFunctionReturn(PETSC_SUCCESS);
    534: }

    536: /*
    537:      KSPDestroy_CG - Frees resources allocated in KSPSetup_CG and clears function
    538:                      compositions from KSPCreate_CG. If adding your own KSP implementation,
    539:                      you must be sure to free all allocated resources here to prevent
    540:                      leaks.
    541: */
    542: PetscErrorCode KSPDestroy_CG(KSP ksp)
    543: {
    544:   KSP_CG *cg = (KSP_CG *)ksp->data;

    546:   PetscFunctionBegin;
    547:   PetscCall(PetscFree4(cg->e, cg->d, cg->ee, cg->dd));
    548:   PetscCall(KSPDestroyDefault(ksp));
    549:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", NULL));
    550:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
    551:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", NULL));
    552:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", NULL));
    553:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
    554:   PetscFunctionReturn(PETSC_SUCCESS);
    555: }

    557: /*
    558:      KSPView_CG - Prints information about the current Krylov method being used.
    559:                   If your Krylov method has special options or flags that information
    560:                   should be printed here.
    561: */
    562: PetscErrorCode KSPView_CG(KSP ksp, PetscViewer viewer)
    563: {
    564:   KSP_CG   *cg = (KSP_CG *)ksp->data;
    565:   PetscBool isascii;

    567:   PetscFunctionBegin;
    568:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    569:   if (isascii) {
    570: #if defined(PETSC_USE_COMPLEX)
    571:     PetscCall(PetscViewerASCIIPrintf(viewer, "  variant %s\n", KSPCGTypes[cg->type]));
    572: #endif
    573:     if (cg->singlereduction) PetscCall(PetscViewerASCIIPrintf(viewer, "  using single-reduction variant\n"));
    574:   }
    575:   PetscFunctionReturn(PETSC_SUCCESS);
    576: }

    578: /*
    579:     KSPSetFromOptions_CG - Checks the options database for options related to the
    580:                            conjugate gradient method.
    581: */
    582: PetscErrorCode KSPSetFromOptions_CG(KSP ksp, PetscOptionItems PetscOptionsObject)
    583: {
    584:   KSP_CG   *cg = (KSP_CG *)ksp->data;
    585:   PetscBool flg;

    587:   PetscFunctionBegin;
    588:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP CG and CGNE options");
    589: #if defined(PETSC_USE_COMPLEX)
    590:   PetscCall(PetscOptionsEnum("-ksp_cg_type", "Matrix is Hermitian or complex symmetric", "KSPCGSetType", KSPCGTypes, (PetscEnum)cg->type, (PetscEnum *)&cg->type, NULL));
    591: #endif
    592:   PetscCall(PetscOptionsBool("-ksp_cg_single_reduction", "Merge inner products into single MPI_Allreduce()", "KSPCGUseSingleReduction", cg->singlereduction, &cg->singlereduction, &flg));
    593:   if (flg) PetscCall(KSPCGUseSingleReduction(ksp, cg->singlereduction));
    594:   PetscOptionsHeadEnd();
    595:   PetscFunctionReturn(PETSC_SUCCESS);
    596: }

    598: /*
    599:     KSPCGSetType_CG - This is an option that is SPECIFIC to this particular Krylov method.
    600:                       This routine is registered below in KSPCreate_CG() and called from the
    601:                       routine KSPCGSetType() (see the file cgtype.c).
    602: */
    603: PetscErrorCode KSPCGSetType_CG(KSP ksp, KSPCGType type)
    604: {
    605:   KSP_CG *cg = (KSP_CG *)ksp->data;

    607:   PetscFunctionBegin;
    608:   cg->type = type;
    609:   PetscFunctionReturn(PETSC_SUCCESS);
    610: }

    612: /*
    613:     KSPCGUseSingleReduction_CG

    615:     This routine sets a flag to use a variant of CG. Note that (in somewhat
    616:     atypical fashion) it also swaps out the routine called when KSPSolve()
    617:     is invoked.
    618: */
    619: static PetscErrorCode KSPCGUseSingleReduction_CG(KSP ksp, PetscBool flg)
    620: {
    621:   KSP_CG *cg = (KSP_CG *)ksp->data;

    623:   PetscFunctionBegin;
    624:   cg->singlereduction = flg;
    625:   if (cg->singlereduction) {
    626:     ksp->ops->solve = KSPSolve_CG_SingleReduction;
    627:     PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
    628:   } else {
    629:     ksp->ops->solve = KSPSolve_CG;
    630:     PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
    631:   }
    632:   PetscFunctionReturn(PETSC_SUCCESS);
    633: }

    635: PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP ksp, Vec t, Vec v, Vec *V)
    636: {
    637:   PetscFunctionBegin;
    638:   PetscCall(VecCopy(ksp->work[0], v));
    639:   *V = v;
    640:   PetscFunctionReturn(PETSC_SUCCESS);
    641: }

    643: /*MC
    644:    KSPCG - The Preconditioned Conjugate Gradient (PCG) iterative method {cite}`hs:52` and {cite}`malek2014preconditioning` for solving linear systems using `KSP`.

    646:    Options Database Keys:
    647: +   -ksp_cg_type (hermitian|symmetric) - (for complex matrices only) indicates the matrix is Hermitian or symmetric, see `KSPCGSetType()`
    648: -   -ksp_cg_single_reduction           - performs both inner products needed in the algorithm with a single `MPI_Allreduce()` call, see `KSPCGUseSingleReduction()`

    650:    Level: beginner

    652:    Notes:
    653:    The `KSPCG` method requires both the matrix and preconditioner to be symmetric positive (or negative) (semi) definite.

    655:    `KSPCG` is the best Krylov method, `KSPType`, when the matrix and preconditioner are symmetric positive definite (SPD).

    657:    Only left preconditioning is supported with `KSPCG`; there are several ways to motivate preconditioned CG, but they all produce the same algorithm.
    658:    One can interpret preconditioning $A$ with $B$ to mean any of the following\:
    659: .vb
    660:    (1) Solve a left-preconditioned system $BAx = Bb $, using $ B^{-1}$ to define an inner product in the algorithm.
    661:    (2) Solve a right-preconditioned system $ABy = b, x = By,$ using $B$ to define an inner product in the algorithm.
    662:    (3) Solve a symmetrically-preconditioned system, $ E^TAEy = E^Tb, x = Ey, $ where $B = EE^T.$
    663:    (4) Solve $Ax=b$ with CG, but use the inner product defined by $B$ to define the method.
    664:    In all cases, the resulting algorithm only requires application of $B$ to vectors, the other inner-product does not appear explicitly in the code
    665: .ve

    667:    For complex numbers there are two different CG methods, one for Hermitian symmetric matrices and one for non-Hermitian symmetric matrices. Use
    668:    `KSPCGSetType()` to indicate which type you are using.

    670:    One can use `KSPSetComputeEigenvalues()` and `KSPComputeEigenvalues()` to compute the eigenvalues of the (preconditioned) operator

    672:    There are two pipelined implementations of CG in PETSc `KSPPIPECG` and `KSPGROPPCG`. These may perform better for very large
    673:    numbers of MPI processes since they overlap communication and computation so the reduction operations in CG, that is inner products and norms,
    674:    do not dominate the compute time.

    676:    Developer Note:
    677:    KSPSolve_CG() should actually query the matrix to determine if it is Hermitian or symmetric and NOT require the user to
    678:    indicate it to the `KSP` object.

    680: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSetComputeEigenvalues()`, `KSPComputeEigenvalues()`,
    681:           `KSPCGSetType()`, `KSPCGUseSingleReduction()`, `KSPPIPECG`, `KSPGROPPCG`
    682: M*/

    684: /*
    685:     KSPCreate_CG - Creates the data structure for the Krylov method CG and sets the
    686:        function pointers for all the routines it needs to call (KSPSolve_CG() etc)

    688:     It must be labeled as PETSC_EXTERN to be dynamically linkable in C++
    689: */
    690: PETSC_EXTERN PetscErrorCode KSPCreate_CG(KSP ksp)
    691: {
    692:   KSP_CG *cg;

    694:   PetscFunctionBegin;
    695:   PetscCall(PetscNew(&cg));
    696:   cg->type    = !PetscDefined(USE_COMPLEX) ? KSP_CG_SYMMETRIC : KSP_CG_HERMITIAN;
    697:   cg->obj_min = 0.0;
    698:   ksp->data   = (void *)cg;

    700:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
    701:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
    702:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
    703:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

    705:   /*
    706:        Sets the functions that are associated with this data structure
    707:        (in C++ this is the same as defining virtual functions)
    708:   */
    709:   ksp->ops->setup          = KSPSetUp_CG;
    710:   ksp->ops->solve          = KSPSolve_CG;
    711:   ksp->ops->destroy        = KSPDestroy_CG;
    712:   ksp->ops->view           = KSPView_CG;
    713:   ksp->ops->setfromoptions = KSPSetFromOptions_CG;
    714:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    715:   ksp->ops->buildresidual  = KSPBuildResidual_CG;

    717:   /*
    718:       Attach the function KSPCGSetType_CG() to this object. The routine
    719:       KSPCGSetType() checks for this attached function and calls it if it finds
    720:       it. (Sort of like a dynamic member function that can be added at run time
    721:   */
    722:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CG));
    723:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", KSPCGUseSingleReduction_CG));
    724:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_CG));
    725:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", KSPCGSetObjectiveTarget_CG));
    726:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
    727:   PetscFunctionReturn(PETSC_SUCCESS);
    728: }