Actual source code: stcg.c


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

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

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

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

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

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

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

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

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

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

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

 55:   cg->o_fcn = 0.0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

312:     alpha = rz / kappa;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

591: /*MC
592:      KSPSTCG -   Code to run conjugate gradient method subject to a constraint on the solution norm.

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

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

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

603:             q(s) = g^T * s + 0.5 * s^T * H * s

605:    subject to the trust region constraint

607:             || s || <= delta,

609:    where

611:      delta is the trust region radius,
612:      g is the gradient vector,
613:      H is the Hessian approximation, and
614:      M is the positive definite preconditioner matrix.

616:    `KSPConvergedReason` may be
617: .vb
618:    KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
619:    KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step,
620: .ve
621:    other `KSP` converged/diverged reasons

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

625:    References:
626: +  * - Steihaug, T. (1983): The conjugate gradient method and trust regions in large scale optimization. SIAM J. Numer. Anal. 20, 626--637
627: -  * - Toint, Ph.L. (1981): Towards an efficient sparsity exploiting Newton method for minimization. In: Duff, I., ed., Sparse Matrices and Their Uses, pp. 57--88. Academic Press

629:    Level: developer

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

634: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
635: {
636:   KSPCG_STCG *cg;

638:   PetscFunctionBegin;
639:   PetscCall(PetscNew(&cg));

641:   cg->radius = 0.0;
642:   cg->dtype  = STCG_UNPRECONDITIONED_DIRECTION;

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

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

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

664:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_STCG));
665:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_STCG));
666:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_STCG));
667:   PetscFunctionReturn(PETSC_SUCCESS);
668: }