Actual source code: qcg.c

      1: #include <../src/ksp/ksp/impls/qcg/qcgimpl.h>

      3: /*
      4:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
      5:          ||s + step*p|| - delta = 0
      6:    such that step1 >= 0 >= step2.
      7:    where
      8:       delta:
      9:         On entry delta must contain scalar delta.
     10:         On exit delta is unchanged.
     11:       step1:
     12:         On entry step1 need not be specified.
     13:         On exit step1 contains the non-negative root.
     14:       step2:
     15:         On entry step2 need not be specified.
     16:         On exit step2 contains the non-positive root.
     17:    C code is translated from the Fortran version of the MINPACK-2 Project,
     18:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
     19: */
     20: static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
     21: {
     22:   PetscReal dsq, ptp, pts, rad, sts;

     24:   PetscFunctionBegin;
     25:   PetscCall(VecDotRealPart(p, s, &pts));
     26:   PetscCall(VecDotRealPart(p, p, &ptp));
     27:   PetscCall(VecDotRealPart(s, s, &sts));
     28:   dsq = delta * delta;
     29:   rad = PetscSqrtReal((pts * pts) - ptp * (sts - dsq));
     30:   if (pts > 0.0) {
     31:     *step2 = -(pts + rad) / ptp;
     32:     *step1 = (sts - dsq) / (ptp * *step2);
     33:   } else {
     34:     *step1 = -(pts - rad) / ptp;
     35:     *step2 = (sts - dsq) / (ptp * *step1);
     36:   }
     37:   PetscFunctionReturn(PETSC_SUCCESS);
     38: }

     40: /*@
     41:   KSPQCGSetTrustRegionRadius - Sets the radius of the trust region for `KSPQCG`

     43:   Logically Collective

     45:   Input Parameters:
     46: + ksp   - the iterative context
     47: - delta - the trust region radius (Infinity is the default)

     49:   Options Database Key:
     50: . -ksp_qcg_trustregionradius delta - trust region radius

     52:   Level: advanced

     54:   Developer Note:
     55:   `KSPMINRESSetRadius()`, for example, does not have TrustRegion in the name

     57: .seealso: [](ch_ksp), `KSPQCG`, `KSPQCGGetTrialStepNorm()`
     58: @*/
     59: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
     60: {
     61:   PetscFunctionBegin;
     63:   PetscCheck(delta >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
     64:   PetscTryMethod(ksp, "KSPQCGSetTrustRegionRadius_C", (KSP, PetscReal), (ksp, delta));
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: /*@
     69:   KSPQCGGetTrialStepNorm - Gets the norm of a trial step vector in `KSPQCG`.  The WCG step may be
     70:   constrained, so this is not necessarily the length of the ultimate step taken in `KSPQCG`.

     72:   Not Collective

     74:   Input Parameter:
     75: . ksp - the iterative context

     77:   Output Parameter:
     78: . tsnorm - the norm

     80:   Level: advanced

     82: .seealso: [](ch_ksp), `KSPQCG`, `KSPQCGSetTrustRegionRadius()`
     83: @*/
     84: PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
     85: {
     86:   PetscFunctionBegin;
     88:   PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: /*@
     93:   KSPQCGGetQuadratic - Gets the value of the quadratic function, evaluated at the new iterate

     95:   Collective

     97:   Input Parameter:
     98: . ksp - the iterative context

    100:   Output Parameter:
    101: . quadratic - the quadratic function evaluated at the new iterate

    103:   Level: advanced

    105:   Note:
    106:   The quadratic function is

    108:   $$
    109:   q(s) = g^T * s + 0.5 * s^T * H * s
    110:   $$

    112:   which satisfies the Euclidean Norm trust region constraint

    114:   $$
    115:   || D * s || \le delta,
    116:   $$

    118:   where
    119: .vb
    120:   delta is the trust region radius,
    121:   g is the gradient vector, and
    122:   H is Hessian matrix,
    123:   D is a scaling matrix.
    124: .ve

    126: .seealso: [](ch_ksp), `KSPQCG`
    127: @*/
    128: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
    129: {
    130:   PetscFunctionBegin;
    132:   PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
    133:   PetscFunctionReturn(PETSC_SUCCESS);
    134: }

    136: static PetscErrorCode KSPSolve_QCG(KSP ksp)
    137: {
    138:   /*
    139:    Correspondence with documentation above:
    140:       B = g = gradient,
    141:       X = s = step
    142:    Note:  This is not coded correctly for complex arithmetic!
    143:  */

    145:   KSP_QCG    *pcgP = (KSP_QCG *)ksp->data;
    146:   Mat         Amat, Pmat;
    147:   Vec         W, WA, WA2, R, P, ASP, BS, X, B;
    148:   PetscScalar scal, beta, rntrn, step;
    149:   PetscReal   q1, q2, xnorm, step1, step2, rnrm = 0.0, btx, xtax;
    150:   PetscReal   ptasp, rtr, wtasp, bstp;
    151:   PetscReal   dzero = 0.0, bsnrm = 0.0;
    152:   PetscInt    i, maxit;
    153:   PC          pc = ksp->pc;
    154:   PetscBool   diagonalscale;

    156:   PetscFunctionBegin;
    157:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
    158:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
    159:   PetscCheck(!ksp->transpose_solve, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Currently does not support transpose solve");

    161:   ksp->its = 0;
    162:   maxit    = ksp->max_it;
    163:   WA       = ksp->work[0];
    164:   R        = ksp->work[1];
    165:   P        = ksp->work[2];
    166:   ASP      = ksp->work[3];
    167:   BS       = ksp->work[4];
    168:   W        = ksp->work[5];
    169:   WA2      = ksp->work[6];
    170:   X        = ksp->vec_sol;
    171:   B        = ksp->vec_rhs;

    173:   PetscCheck(pcgP->delta > dzero, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: delta <= 0");

    175:   /* Initialize variables */
    176:   PetscCall(VecSet(W, 0.0)); /* W = 0 */
    177:   PetscCall(VecSet(X, 0.0)); /* X = 0 */
    178:   PetscCall(PCGetOperators(pc, &Amat, &Pmat));

    180:   /* Compute:  BS = D^{-1} B */
    181:   PetscCall(PCApplySymmetricLeft(pc, B, BS));

    183:   if (ksp->normtype != KSP_NORM_NONE) {
    184:     PetscCall(VecNorm(BS, NORM_2, &bsnrm));
    185:     KSPCheckNorm(ksp, bsnrm);
    186:   }
    187:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
    188:   ksp->its   = 0;
    189:   ksp->rnorm = bsnrm;
    190:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
    191:   PetscCall(KSPLogResidualHistory(ksp, bsnrm));
    192:   PetscCall(KSPMonitor(ksp, 0, bsnrm));
    193:   PetscCall((*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP));
    194:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

    196:   /* Compute the initial scaled direction and scaled residual */
    197:   PetscCall(VecCopy(BS, R));
    198:   PetscCall(VecScale(R, -1.0));
    199:   PetscCall(VecCopy(R, P));
    200:   PetscCall(VecDotRealPart(R, R, &rtr));

    202:   for (i = 0; i <= maxit; i++) {
    203:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
    204:     ksp->its++;
    205:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

    207:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
    208:     PetscCall(PCApplySymmetricRight(pc, P, WA));
    209:     PetscCall(KSP_MatMult(ksp, Amat, WA, WA2));
    210:     PetscCall(PCApplySymmetricLeft(pc, WA2, ASP));

    212:     /* Check for negative curvature */
    213:     PetscCall(VecDotRealPart(P, ASP, &ptasp));
    214:     if (ptasp <= dzero) {
    215:       /* Scaled negative curvature direction:  Compute a step so that
    216:         ||w + step*p|| = delta and QS(w + step*p) is least */

    218:       if (!i) {
    219:         PetscCall(VecCopy(P, X));
    220:         PetscCall(VecNorm(X, NORM_2, &xnorm));
    221:         KSPCheckNorm(ksp, xnorm);
    222:         scal = pcgP->delta / xnorm;
    223:         PetscCall(VecScale(X, scal));
    224:       } else {
    225:         /* Compute roots of quadratic */
    226:         PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
    227:         PetscCall(VecDotRealPart(W, ASP, &wtasp));
    228:         PetscCall(VecDotRealPart(BS, P, &bstp));
    229:         PetscCall(VecCopy(W, X));
    230:         q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
    231:         q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
    232:         if (q1 <= q2) PetscCall(VecAXPY(X, step1, P));
    233:         else PetscCall(VecAXPY(X, step2, P));
    234:       }
    235:       pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
    236:       ksp->reason  = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    237:       if (!i) {
    238:         PetscCall(PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta));
    239:       } else {
    240:         PetscCall(PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
    241:       }

    243:     } else {
    244:       /* Compute step along p */
    245:       step = rtr / ptasp;
    246:       PetscCall(VecCopy(W, X));       /*  x = w  */
    247:       PetscCall(VecAXPY(X, step, P)); /*  x <- step*p + x  */
    248:       PetscCall(VecNorm(X, NORM_2, &pcgP->ltsnrm));
    249:       KSPCheckNorm(ksp, pcgP->ltsnrm);

    251:       if (pcgP->ltsnrm > pcgP->delta) {
    252:         /* Since the trial iterate is outside the trust region,
    253:             evaluate a constrained step along p so that
    254:                     ||w + step*p|| = delta
    255:           The positive step is always better in this case. */
    256:         if (!i) {
    257:           scal = pcgP->delta / pcgP->ltsnrm;
    258:           PetscCall(VecScale(X, scal));
    259:         } else {
    260:           /* Compute roots of quadratic */
    261:           PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
    262:           PetscCall(VecCopy(W, X));
    263:           PetscCall(VecAXPY(X, step1, P)); /*  x <- step1*p + x  */
    264:         }
    265:         pcgP->ltsnrm = pcgP->delta;
    266:         ksp->reason  = KSP_CONVERGED_STEP_LENGTH; /* convergence along constrained step */
    267:         if (!i) {
    268:           PetscCall(PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta));
    269:         } else {
    270:           PetscCall(PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
    271:         }

    273:       } else {
    274:         /* Evaluate the current step */
    275:         PetscCall(VecCopy(X, W));          /* update interior iterate */
    276:         PetscCall(VecAXPY(R, -step, ASP)); /* r <- -step*asp + r */
    277:         if (ksp->normtype != KSP_NORM_NONE) {
    278:           PetscCall(VecNorm(R, NORM_2, &rnrm));
    279:           KSPCheckNorm(ksp, rnrm);
    280:         }
    281:         PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
    282:         ksp->rnorm = rnrm;
    283:         PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
    284:         PetscCall(KSPLogResidualHistory(ksp, rnrm));
    285:         PetscCall(KSPMonitor(ksp, i + 1, rnrm));
    286:         PetscCall((*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP));
    287:         if (ksp->reason) { /* convergence for */
    288:           PetscCall(PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta));
    289:         }
    290:       }
    291:     }
    292:     if (ksp->reason) break; /* Convergence has been attained */
    293:     else { /* Compute a new AS-orthogonal direction */ PetscCall(VecDot(R, R, &rntrn));
    294:       beta = rntrn / rtr;
    295:       PetscCall(VecAYPX(P, beta, R)); /*  p <- r + beta*p  */
    296:       rtr = PetscRealPart(rntrn);
    297:     }
    298:   }
    299:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

    301:   /* Unscale x */
    302:   PetscCall(VecCopy(X, WA2));
    303:   PetscCall(PCApplySymmetricRight(pc, WA2, X));

    305:   PetscCall(KSP_MatMult(ksp, Amat, X, WA));
    306:   PetscCall(VecDotRealPart(B, X, &btx));
    307:   PetscCall(VecDotRealPart(X, WA, &xtax));

    309:   pcgP->quadratic = btx + .5 * xtax;
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: static PetscErrorCode KSPSetUp_QCG(KSP ksp)
    314: {
    315:   PetscFunctionBegin;
    316:   /* Get work vectors from user code */
    317:   PetscCall(KSPSetWorkVecs(ksp, 7));
    318:   PetscFunctionReturn(PETSC_SUCCESS);
    319: }

    321: static PetscErrorCode KSPDestroy_QCG(KSP ksp)
    322: {
    323:   PetscFunctionBegin;
    324:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL));
    325:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL));
    326:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL));
    327:   PetscCall(KSPDestroyDefault(ksp));
    328:   PetscFunctionReturn(PETSC_SUCCESS);
    329: }

    331: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
    332: {
    333:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

    335:   PetscFunctionBegin;
    336:   cgP->delta = delta;
    337:   PetscFunctionReturn(PETSC_SUCCESS);
    338: }

    340: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
    341: {
    342:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

    344:   PetscFunctionBegin;
    345:   *ltsnrm = cgP->ltsnrm;
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
    350: {
    351:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

    353:   PetscFunctionBegin;
    354:   *quadratic = cgP->quadratic;
    355:   PetscFunctionReturn(PETSC_SUCCESS);
    356: }

    358: static PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems PetscOptionsObject)
    359: {
    360:   PetscReal delta;
    361:   KSP_QCG  *cgP = (KSP_QCG *)ksp->data;
    362:   PetscBool flg;

    364:   PetscFunctionBegin;
    365:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
    366:   PetscCall(PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg));
    367:   if (flg) PetscCall(KSPQCGSetTrustRegionRadius(ksp, delta));
    368:   PetscOptionsHeadEnd();
    369:   PetscFunctionReturn(PETSC_SUCCESS);
    370: }

    372: /*MC
    373:    KSPQCG - Code to run conjugate gradient method subject to a constraint on the solution norm  {cite}`steihaug:83`.

    375:    Options Database Key:
    376: .  -ksp_qcg_trustregionradius r - Trust Region Radius

    378:    Level: developer

    380:    Notes:
    381:    This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`

    383:    Uses preconditioned conjugate gradient to compute
    384:    an approximate minimizer of the quadratic function $ q(s) = g^T * s + .5 * s^T * H * s $   subject to the Euclidean norm trust region constraint
    385:    $ || D * s || \le delta$, where
    386: .vb
    387:      delta is the trust region radius,
    388:      g is the gradient vector, and
    389:      H is Hessian matrix,
    390:      D is a scaling matrix.
    391: .ve

    393:    `KSPConvergedReason` may include
    394: +  `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
    395: -  `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,

    397:   Note:
    398:   Allows symmetric preconditioning with the following scaling matrices:
    399: .vb
    400:       `PCNONE`:   D = Identity matrix
    401:       `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
    402:       `PCICC`:    D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
    403: .ve

    405: .seealso: [](ch_ksp), `KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`,
    406:           `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
    407: M*/

    409: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
    410: {
    411:   KSP_QCG *cgP;

    413:   PetscFunctionBegin;
    414:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3));
    415:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1));
    416:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
    417:   PetscCall(PetscNew(&cgP));

    419:   ksp->data                = (void *)cgP;
    420:   ksp->ops->setup          = KSPSetUp_QCG;
    421:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
    422:   ksp->ops->solve          = KSPSolve_QCG;
    423:   ksp->ops->destroy        = KSPDestroy_QCG;
    424:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    425:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    426:   ksp->ops->view           = NULL;

    428:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG));
    429:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG));
    430:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG));
    431:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
    432:   PetscFunctionReturn(PETSC_SUCCESS);
    433: }