Actual source code: stcg.c

  1: #include <../src/ksp/ksp/impls/cg/stcg/stcgimpl.h>

  3: #define STCG_PRECONDITIONED_DIRECTION   0
  4: #define STCG_UNPRECONDITIONED_DIRECTION 1
  5: #define STCG_DIRECTION_TYPES            2

  7: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};

  9: static PetscErrorCode KSPCGSolve_STCG(KSP ksp)
 10: {
 11: #if defined(PETSC_USE_COMPLEX)
 12:   SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "STCG is not available for complex systems");
 13: #else
 14:   KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
 15:   Mat         Qmat, Mmat;
 16:   Vec         r, z, p, d;
 17:   PC          pc;
 18:   PetscReal   norm_r, norm_d, norm_dp1, norm_p, dMp;
 19:   PetscReal   alpha, beta, kappa, rz, rzm1;
 20:   PetscReal   rr, r2, step;
 21:   PetscInt    max_cg_its;
 22:   PetscBool   diagonalscale;

 24:   /***************************************************************************/
 25:   /* Check the arguments and parameters.                                     */
 26:   /***************************************************************************/

 28:   PetscFunctionBegin;
 29:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 30:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
 31:   PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");

 33:   /***************************************************************************/
 34:   /* Get the workspace vectors and initialize variables                      */
 35:   /***************************************************************************/

 37:   r2 = cg->radius * cg->radius;
 38:   r  = ksp->work[0];
 39:   z  = ksp->work[1];
 40:   p  = ksp->work[2];
 41:   d  = ksp->vec_sol;
 42:   pc = ksp->pc;

 44:   PetscCall(PCGetOperators(pc, &Qmat, &Mmat));

 46:   PetscCall(VecGetSize(d, &max_cg_its));
 47:   max_cg_its = PetscMin(max_cg_its, ksp->max_it);
 48:   ksp->its   = 0;

 50:   /***************************************************************************/
 51:   /* Initialize objective function and direction.                            */
 52:   /***************************************************************************/

 54:   cg->o_fcn = 0.0;

 56:   PetscCall(VecSet(d, 0.0)); /* d = 0             */
 57:   cg->norm_d = 0.0;

 59:   /***************************************************************************/
 60:   /* Begin the conjugate gradient method.  Check the right-hand side for     */
 61:   /* numerical problems.  The check for not-a-number and infinite values     */
 62:   /* need be performed only once.                                            */
 63:   /***************************************************************************/

 65:   PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad         */
 66:   PetscCall(VecDot(r, r, &rr));        /* rr = r^T r        */
 67:   KSPCheckDot(ksp, rr);

 69:   /***************************************************************************/
 70:   /* Check the preconditioner for numerical problems and for positive        */
 71:   /* definiteness.  The check for not-a-number and infinite values need be   */
 72:   /* performed only once.                                                    */
 73:   /***************************************************************************/

 75:   PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */
 76:   PetscCall(VecDot(r, z, &rz));      /* rz = r^T inv(M) r */
 77:   if (PetscIsInfOrNanScalar(rz)) {
 78:     /*************************************************************************/
 79:     /* The preconditioner contains not-a-number or an infinite value.        */
 80:     /* Return the gradient direction intersected with the trust region.      */
 81:     /*************************************************************************/

 83:     ksp->reason = KSP_DIVERGED_NANORINF;
 84:     PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad preconditioner: rz=%g\n", (double)rz));

 86:     if (cg->radius != 0) {
 87:       if (r2 >= rr) {
 88:         alpha      = 1.0;
 89:         cg->norm_d = PetscSqrtReal(rr);
 90:       } else {
 91:         alpha      = PetscSqrtReal(r2 / rr);
 92:         cg->norm_d = cg->radius;
 93:       }

 95:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

 97:       /***********************************************************************/
 98:       /* Compute objective function.                                         */
 99:       /***********************************************************************/

101:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
102:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
103:       PetscCall(VecDot(d, z, &cg->o_fcn));
104:       cg->o_fcn = -cg->o_fcn;
105:       ++ksp->its;
106:     }
107:     PetscFunctionReturn(PETSC_SUCCESS);
108:   }

110:   if (rz < 0.0) {
111:     /*************************************************************************/
112:     /* The preconditioner is indefinite.  Because this is the first          */
113:     /* and we do not have a direction yet, we use the gradient step.  Note   */
114:     /* that we cannot use the preconditioned norm when computing the step    */
115:     /* because the matrix is indefinite.                                     */
116:     /*************************************************************************/

118:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
119:     PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz));

121:     if (cg->radius != 0.0) {
122:       if (r2 >= rr) {
123:         alpha      = 1.0;
124:         cg->norm_d = PetscSqrtReal(rr);
125:       } else {
126:         alpha      = PetscSqrtReal(r2 / rr);
127:         cg->norm_d = cg->radius;
128:       }

130:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

132:       /***********************************************************************/
133:       /* Compute objective function.                                         */
134:       /***********************************************************************/

136:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
137:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
138:       PetscCall(VecDot(d, z, &cg->o_fcn));
139:       cg->o_fcn = -cg->o_fcn;
140:       ++ksp->its;
141:     }
142:     PetscFunctionReturn(PETSC_SUCCESS);
143:   }

145:   /***************************************************************************/
146:   /* As far as we know, the preconditioner is positive semidefinite.         */
147:   /* Compute and log the residual.  Check convergence because this           */
148:   /* initializes things, but do not terminate until at least one conjugate   */
149:   /* gradient iteration has been performed.                                  */
150:   /***************************************************************************/

152:   switch (ksp->normtype) {
153:   case KSP_NORM_PRECONDITIONED:
154:     PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
155:     break;

157:   case KSP_NORM_UNPRECONDITIONED:
158:     norm_r = PetscSqrtReal(rr); /* norm_r = |r|      */
159:     break;

161:   case KSP_NORM_NATURAL:
162:     norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M    */
163:     break;

165:   default:
166:     norm_r = 0.0;
167:     break;
168:   }

170:   PetscCall(KSPLogResidualHistory(ksp, norm_r));
171:   PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
172:   ksp->rnorm = norm_r;

174:   PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));

176:   /***************************************************************************/
177:   /* Compute the first direction and update the iteration.                   */
178:   /***************************************************************************/

180:   PetscCall(VecCopy(z, p));                /* p = z             */
181:   PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
182:   ++ksp->its;

184:   /***************************************************************************/
185:   /* Check the matrix for numerical problems.                                */
186:   /***************************************************************************/

188:   PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p   */
189:   if (PetscIsInfOrNanScalar(kappa)) {
190:     /*************************************************************************/
191:     /* The matrix produced not-a-number or an infinite value.  In this case, */
192:     /* we must stop and use the gradient direction.  This condition need     */
193:     /* only be checked once.                                                 */
194:     /*************************************************************************/

196:     ksp->reason = KSP_DIVERGED_NANORINF;
197:     PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad matrix: kappa=%g\n", (double)kappa));

199:     if (cg->radius) {
200:       if (r2 >= rr) {
201:         alpha      = 1.0;
202:         cg->norm_d = PetscSqrtReal(rr);
203:       } else {
204:         alpha      = PetscSqrtReal(r2 / rr);
205:         cg->norm_d = cg->radius;
206:       }

208:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

210:       /***********************************************************************/
211:       /* Compute objective function.                                         */
212:       /***********************************************************************/

214:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
215:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
216:       PetscCall(VecDot(d, z, &cg->o_fcn));
217:       cg->o_fcn = -cg->o_fcn;
218:       ++ksp->its;
219:     }
220:     PetscFunctionReturn(PETSC_SUCCESS);
221:   }

223:   /***************************************************************************/
224:   /* Initialize variables for calculating the norm of the direction.         */
225:   /***************************************************************************/

227:   dMp    = 0.0;
228:   norm_d = 0.0;
229:   switch (cg->dtype) {
230:   case STCG_PRECONDITIONED_DIRECTION:
231:     norm_p = rz;
232:     break;

234:   default:
235:     PetscCall(VecDot(p, p, &norm_p));
236:     break;
237:   }

239:   /***************************************************************************/
240:   /* Check for negative curvature.                                           */
241:   /***************************************************************************/

243:   if (kappa <= 0.0) {
244:     /*************************************************************************/
245:     /* In this case, the matrix is indefinite and we have encountered a      */
246:     /* direction of negative curvature.  Because negative curvature occurs   */
247:     /* during the first step, we must follow a direction.                    */
248:     /*************************************************************************/

250:     ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
251:     PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));

253:     if (cg->radius != 0.0 && norm_p > 0.0) {
254:       /***********************************************************************/
255:       /* Follow direction of negative curvature to the boundary of the       */
256:       /* trust region.                                                       */
257:       /***********************************************************************/

259:       step       = PetscSqrtReal(r2 / norm_p);
260:       cg->norm_d = cg->radius;

262:       PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

264:       /***********************************************************************/
265:       /* Update objective function.                                          */
266:       /***********************************************************************/

268:       cg->o_fcn += step * (0.5 * step * kappa - rz);
269:     } else if (cg->radius != 0.0) {
270:       /***********************************************************************/
271:       /* The norm of the preconditioned direction is zero; use the gradient  */
272:       /* step.                                                               */
273:       /***********************************************************************/

275:       if (r2 >= rr) {
276:         alpha      = 1.0;
277:         cg->norm_d = PetscSqrtReal(rr);
278:       } else {
279:         alpha      = PetscSqrtReal(r2 / rr);
280:         cg->norm_d = cg->radius;
281:       }

283:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

285:       /***********************************************************************/
286:       /* Compute objective function.                                         */
287:       /***********************************************************************/

289:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
290:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
291:       PetscCall(VecDot(d, z, &cg->o_fcn));

293:       cg->o_fcn = -cg->o_fcn;
294:       ++ksp->its;
295:     }
296:     PetscFunctionReturn(PETSC_SUCCESS);
297:   }

299:   /***************************************************************************/
300:   /* Run the conjugate gradient method until either the problem is solved,   */
301:   /* we encounter the boundary of the trust region, or the conjugate         */
302:   /* gradient method breaks down.                                            */
303:   /***************************************************************************/

305:   while (1) {
306:     /*************************************************************************/
307:     /* Know that kappa is nonzero, because we have not broken down, so we    */
308:     /* can compute the steplength.                                           */
309:     /*************************************************************************/

311:     alpha = rz / kappa;

313:     /*************************************************************************/
314:     /* Compute the steplength and check for intersection with the trust      */
315:     /* region.                                                               */
316:     /*************************************************************************/

318:     norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
319:     if (cg->radius != 0.0 && norm_dp1 >= r2) {
320:       /***********************************************************************/
321:       /* In this case, the matrix is positive definite as far as we know.    */
322:       /* However, the full step goes beyond the trust region.                */
323:       /***********************************************************************/

325:       ksp->reason = KSP_CONVERGED_STEP_LENGTH;
326:       PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: constrained step: radius=%g\n", (double)cg->radius));

328:       if (norm_p > 0.0) {
329:         /*********************************************************************/
330:         /* Follow the direction to the boundary of the trust region.         */
331:         /*********************************************************************/

333:         step       = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
334:         cg->norm_d = cg->radius;

336:         PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

338:         /*********************************************************************/
339:         /* Update objective function.                                        */
340:         /*********************************************************************/

342:         cg->o_fcn += step * (0.5 * step * kappa - rz);
343:       } else {
344:         /*********************************************************************/
345:         /* The norm of the direction is zero; there is nothing to follow.    */
346:         /*********************************************************************/
347:       }
348:       break;
349:     }

351:     /*************************************************************************/
352:     /* Now we can update the direction and residual.                         */
353:     /*************************************************************************/

355:     PetscCall(VecAXPY(d, alpha, p));   /* d = d + alpha p   */
356:     PetscCall(VecAXPY(r, -alpha, z));  /* r = r - alpha Q p */
357:     PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */

359:     switch (cg->dtype) {
360:     case STCG_PRECONDITIONED_DIRECTION:
361:       norm_d = norm_dp1;
362:       break;

364:     default:
365:       PetscCall(VecDot(d, d, &norm_d));
366:       break;
367:     }
368:     cg->norm_d = PetscSqrtReal(norm_d);

370:     /*************************************************************************/
371:     /* Update objective function.                                            */
372:     /*************************************************************************/

374:     cg->o_fcn -= 0.5 * alpha * rz;

376:     /*************************************************************************/
377:     /* Check that the preconditioner appears positive semidefinite.          */
378:     /*************************************************************************/

380:     rzm1 = rz;
381:     PetscCall(VecDot(r, z, &rz)); /* rz = r^T z        */
382:     if (rz < 0.0) {
383:       /***********************************************************************/
384:       /* The preconditioner is indefinite.                                   */
385:       /***********************************************************************/

387:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
388:       PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz));
389:       break;
390:     }

392:     /*************************************************************************/
393:     /* As far as we know, the preconditioner is positive semidefinite.       */
394:     /* Compute the residual and check for convergence.                       */
395:     /*************************************************************************/

397:     switch (ksp->normtype) {
398:     case KSP_NORM_PRECONDITIONED:
399:       PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
400:       break;

402:     case KSP_NORM_UNPRECONDITIONED:
403:       PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r|      */
404:       break;

406:     case KSP_NORM_NATURAL:
407:       norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M    */
408:       break;

410:     default:
411:       norm_r = 0.0;
412:       break;
413:     }

415:     PetscCall(KSPLogResidualHistory(ksp, norm_r));
416:     PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
417:     ksp->rnorm = norm_r;

419:     PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
420:     if (ksp->reason) {
421:       /***********************************************************************/
422:       /* The method has converged.                                           */
423:       /***********************************************************************/

425:       PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
426:       break;
427:     }

429:     /*************************************************************************/
430:     /* We have not converged yet.  Check for breakdown.                      */
431:     /*************************************************************************/

433:     beta = rz / rzm1;
434:     if (PetscAbsScalar(beta) <= 0.0) {
435:       /***********************************************************************/
436:       /* Conjugate gradients has broken down.                                */
437:       /***********************************************************************/

439:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
440:       PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: breakdown: beta=%g\n", (double)beta));
441:       break;
442:     }

444:     /*************************************************************************/
445:     /* Check iteration limit.                                                */
446:     /*************************************************************************/

448:     if (ksp->its >= max_cg_its) {
449:       ksp->reason = KSP_DIVERGED_ITS;
450:       PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
451:       break;
452:     }

454:     /*************************************************************************/
455:     /* Update p and the norms.                                               */
456:     /*************************************************************************/

458:     PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p    */

460:     switch (cg->dtype) {
461:     case STCG_PRECONDITIONED_DIRECTION:
462:       dMp    = beta * (dMp + alpha * norm_p);
463:       norm_p = beta * (rzm1 + beta * norm_p);
464:       break;

466:     default:
467:       PetscCall(VecDot(d, p, &dMp));
468:       PetscCall(VecDot(p, p, &norm_p));
469:       break;
470:     }

472:     /*************************************************************************/
473:     /* Compute the new direction and update the iteration.                   */
474:     /*************************************************************************/

476:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
477:     PetscCall(VecDot(p, z, &kappa));         /* kappa = p^T Q p   */
478:     ++ksp->its;

480:     /*************************************************************************/
481:     /* Check for negative curvature.                                         */
482:     /*************************************************************************/

484:     if (kappa <= 0.0) {
485:       /***********************************************************************/
486:       /* In this case, the matrix is indefinite and we have encountered      */
487:       /* a direction of negative curvature.  Follow the direction to the     */
488:       /* boundary of the trust region.                                       */
489:       /***********************************************************************/

491:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
492:       PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));

494:       if (cg->radius != 0.0 && norm_p > 0.0) {
495:         /*********************************************************************/
496:         /* Follow direction of negative curvature to boundary.               */
497:         /*********************************************************************/

499:         step       = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
500:         cg->norm_d = cg->radius;

502:         PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

504:         /*********************************************************************/
505:         /* Update objective function.                                        */
506:         /*********************************************************************/

508:         cg->o_fcn += step * (0.5 * step * kappa - rz);
509:       } else if (cg->radius != 0.0) {
510:         /*********************************************************************/
511:         /* The norm of the direction is zero; there is nothing to follow.    */
512:         /*********************************************************************/
513:       }
514:       break;
515:     }
516:   }
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: #endif
519: }

521: static PetscErrorCode KSPCGSetUp_STCG(KSP ksp)
522: {
523:   PetscFunctionBegin;
524:   /***************************************************************************/
525:   /* Set work vectors needed by conjugate gradient method and allocate       */
526:   /***************************************************************************/

528:   PetscCall(KSPSetWorkVecs(ksp, 3));
529:   PetscFunctionReturn(PETSC_SUCCESS);
530: }

532: static PetscErrorCode KSPCGDestroy_STCG(KSP ksp)
533: {
534:   PetscFunctionBegin;
535:   /***************************************************************************/
536:   /* Clear composed functions                                                */
537:   /***************************************************************************/

539:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
540:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
541:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));

543:   /***************************************************************************/
544:   /* Destroy KSP object.                                                     */
545:   /***************************************************************************/

547:   PetscCall(KSPDestroyDefault(ksp));
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

551: static PetscErrorCode KSPCGSetRadius_STCG(KSP ksp, PetscReal radius)
552: {
553:   KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;

555:   PetscFunctionBegin;
556:   cg->radius = radius;
557:   PetscFunctionReturn(PETSC_SUCCESS);
558: }

560: static PetscErrorCode KSPCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
561: {
562:   KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;

564:   PetscFunctionBegin;
565:   *norm_d = cg->norm_d;
566:   PetscFunctionReturn(PETSC_SUCCESS);
567: }

569: static PetscErrorCode KSPCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
570: {
571:   KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;

573:   PetscFunctionBegin;
574:   *o_fcn = cg->o_fcn;
575:   PetscFunctionReturn(PETSC_SUCCESS);
576: }

578: static PetscErrorCode KSPCGSetFromOptions_STCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
579: {
580:   KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;

582:   PetscFunctionBegin;
583:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG STCG options");
584:   PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
585:   PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
586:   PetscOptionsHeadEnd();
587:   PetscFunctionReturn(PETSC_SUCCESS);
588: }

590: /*MC
591:    KSPSTCG -   Code to run conjugate gradient method subject to a constraint on the solution norm for use in a trust region method
592:    {cite}`steihaug:83`, {cite}`toint1981towards`

594:    Options Database Key:
595: .      -ksp_cg_radius <r> - Trust Region Radius

597:    Level: developer

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

602:    Use preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function

604:    $$
605:    q(s) = g^T * s + 0.5 * s^T * H * s
606:    $$

608:    subject to the trust region constraint

610:    $$
611:             || s || le delta,
612:    $$

614:    where
615: .vb
616:      delta is the trust region radius,
617:      g is the gradient vector,
618:      H is the Hessian approximation, and
619:      M is the positive definite preconditioner matrix.
620: .ve

622:    `KSPConvergedReason` may include
623: +  `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
624: -  `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,

626:   The preconditioner supplied should be symmetric and positive definite.

628: .seealso: [](ch_ksp), `KSPCreate()`, `KSPCGSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPNASH`, `KSPGLTR`, `KSPQCG`
629: M*/

631: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
632: {
633:   KSPCG_STCG *cg;

635:   PetscFunctionBegin;
636:   PetscCall(PetscNew(&cg));

638:   cg->radius = 0.0;
639:   cg->dtype  = STCG_UNPRECONDITIONED_DIRECTION;

641:   ksp->data = (void *)cg;
642:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
643:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
644:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
645:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
646:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));

648:   /***************************************************************************/
649:   /* Sets the functions that are associated with this data structure         */
650:   /* (in C++ this is the same as defining virtual functions).                */
651:   /***************************************************************************/

653:   ksp->ops->setup          = KSPCGSetUp_STCG;
654:   ksp->ops->solve          = KSPCGSolve_STCG;
655:   ksp->ops->destroy        = KSPCGDestroy_STCG;
656:   ksp->ops->setfromoptions = KSPCGSetFromOptions_STCG;
657:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
658:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
659:   ksp->ops->view           = NULL;

661:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_STCG));
662:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_STCG));
663:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_STCG));
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }