Actual source code: nash.c


  2: #include <../src/ksp/ksp/impls/cg/nash/nashimpl.h>

  4: #define NASH_PRECONDITIONED_DIRECTION   0
  5: #define NASH_UNPRECONDITIONED_DIRECTION 1
  6: #define NASH_DIRECTION_TYPES            2

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

 10: static PetscErrorCode KSPCGSolve_NASH(KSP ksp)
 11: {
 12: #if defined(PETSC_USE_COMPLEX)
 13:   SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "NASH is not available for complex systems");
 14: #else
 15:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
 16:   Mat         Qmat, Mmat;
 17:   Vec         r, z, p, d;
 18:   PC          pc;

 20:   PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
 21:   PetscReal alpha, beta, kappa, rz, rzm1;
 22:   PetscReal rr, r2, step;

 24:   PetscInt max_cg_its;

 26:   PetscBool diagonalscale;

 28:   PetscFunctionBegin;
 29:   /***************************************************************************/
 30:   /* Check the arguments and parameters.                                     */
 31:   /***************************************************************************/

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

 37:   /***************************************************************************/
 38:   /* Get the workspace vectors and initialize variables                      */
 39:   /***************************************************************************/

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

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

 50:   PetscCall(VecGetSize(d, &max_cg_its));
 51:   max_cg_its = PetscMin(max_cg_its, ksp->max_it);
 52:   ksp->its   = 0;

 54:   /***************************************************************************/
 55:   /* Initialize objective function and direction.                            */
 56:   /***************************************************************************/

 58:   cg->o_fcn = 0.0;

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

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

 69:   PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad         */
 70:   PetscCall(VecDot(r, r, &rr));        /* rr = r^T r        */
 71:   KSPCheckDot(ksp, rr);

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

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

 87:     ksp->reason = KSP_DIVERGED_NANORINF;
 88:     PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: bad preconditioner: rz=%g\n", (double)rz));

 90:     if (cg->radius) {
 91:       if (r2 >= rr) {
 92:         alpha      = 1.0;
 93:         cg->norm_d = PetscSqrtReal(rr);
 94:       } else {
 95:         alpha      = PetscSqrtReal(r2 / rr);
 96:         cg->norm_d = cg->radius;
 97:       }

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

101:       /***********************************************************************/
102:       /* Compute objective function.                                         */
103:       /***********************************************************************/

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

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

122:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
123:     PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz));

125:     if (cg->radius) {
126:       if (r2 >= rr) {
127:         alpha      = 1.0;
128:         cg->norm_d = PetscSqrtReal(rr);
129:       } else {
130:         alpha      = PetscSqrtReal(r2 / rr);
131:         cg->norm_d = cg->radius;
132:       }

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

136:       /***********************************************************************/
137:       /* Compute objective function.                                         */
138:       /***********************************************************************/

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

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

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

161:   case KSP_NORM_UNPRECONDITIONED:
162:     norm_r = PetscSqrtReal(rr); /* norm_r = |r|      */
163:     break;

165:   case KSP_NORM_NATURAL:
166:     norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M    */
167:     break;

169:   default:
170:     norm_r = 0.0;
171:     break;
172:   }

174:   PetscCall(KSPLogResidualHistory(ksp, norm_r));
175:   PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
176:   ksp->rnorm = norm_r;

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

180:   /***************************************************************************/
181:   /* Compute the first direction and update the iteration.                   */
182:   /***************************************************************************/

184:   PetscCall(VecCopy(z, p));                /* p = z             */
185:   PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
186:   ++ksp->its;

188:   /***************************************************************************/
189:   /* Check the matrix for numerical problems.                                */
190:   /***************************************************************************/

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

200:     ksp->reason = KSP_DIVERGED_NANORINF;
201:     PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: bad matrix: kappa=%g\n", (double)kappa));

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

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

214:       /***********************************************************************/
215:       /* Compute objective function.                                         */
216:       /***********************************************************************/

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

227:   /***************************************************************************/
228:   /* Initialize variables for calculating the norm of the direction.         */
229:   /***************************************************************************/

231:   dMp    = 0.0;
232:   norm_d = 0.0;
233:   switch (cg->dtype) {
234:   case NASH_PRECONDITIONED_DIRECTION:
235:     norm_p = rz;
236:     break;

238:   default:
239:     PetscCall(VecDot(p, p, &norm_p));
240:     break;
241:   }

243:   /***************************************************************************/
244:   /* Check for negative curvature.                                           */
245:   /***************************************************************************/

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

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

257:     if (cg->radius && norm_p > 0.0) {
258:       /***********************************************************************/
259:       /* Follow direction of negative curvature to the boundary of the       */
260:       /* trust region.                                                       */
261:       /***********************************************************************/

263:       step       = PetscSqrtReal(r2 / norm_p);
264:       cg->norm_d = cg->radius;

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

268:       /***********************************************************************/
269:       /* Update objective function.                                          */
270:       /***********************************************************************/

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

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

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

289:       /***********************************************************************/
290:       /* Compute objective function.                                         */
291:       /***********************************************************************/

293:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
294:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
295:       PetscCall(VecDot(d, z, &cg->o_fcn));
296:       cg->o_fcn = -cg->o_fcn;
297:       ++ksp->its;
298:     }
299:     PetscFunctionReturn(PETSC_SUCCESS);
300:   }

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

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

314:     alpha = rz / kappa;

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

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

328:       ksp->reason = KSP_CONVERGED_STEP_LENGTH;
329:       PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: constrained step: radius=%g\n", (double)cg->radius));

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

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

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

341:         /*********************************************************************/
342:         /* Update objective function.                                        */
343:         /*********************************************************************/

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

354:     /*************************************************************************/
355:     /* Now we can update the direction and residual.                         */
356:     /*************************************************************************/

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

362:     switch (cg->dtype) {
363:     case NASH_PRECONDITIONED_DIRECTION:
364:       norm_d = norm_dp1;
365:       break;

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

373:     /*************************************************************************/
374:     /* Update objective function.                                            */
375:     /*************************************************************************/

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

379:     /*************************************************************************/
380:     /* Check that the preconditioner appears positive semidefinite.          */
381:     /*************************************************************************/

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

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

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

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

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

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

413:     default:
414:       norm_r = 0.;
415:       break;
416:     }

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

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

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

432:     /*************************************************************************/
433:     /* We have not converged yet.  Check for breakdown.                      */
434:     /*************************************************************************/

436:     beta = rz / rzm1;
437:     if (PetscAbsReal(beta) <= 0.0) {
438:       /***********************************************************************/
439:       /* Conjugate gradients has broken down.                                */
440:       /***********************************************************************/

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

447:     /*************************************************************************/
448:     /* Check iteration limit.                                                */
449:     /*************************************************************************/

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

457:     /*************************************************************************/
458:     /* Update p and the norms.                                               */
459:     /*************************************************************************/

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

463:     switch (cg->dtype) {
464:     case NASH_PRECONDITIONED_DIRECTION:
465:       dMp    = beta * (dMp + alpha * norm_p);
466:       norm_p = beta * (rzm1 + beta * norm_p);
467:       break;

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

475:     /*************************************************************************/
476:     /* Compute the new direction and update the iteration.                   */
477:     /*************************************************************************/

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

483:     /*************************************************************************/
484:     /* Check for negative curvature.                                         */
485:     /*************************************************************************/

487:     if (kappa <= 0.0) {
488:       /***********************************************************************/
489:       /* In this case, the matrix is indefinite and we have encountered      */
490:       /* a direction of negative curvature.  Stop at the base.               */
491:       /***********************************************************************/

493:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
494:       PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa));
495:       break;
496:     }
497:   }
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: #endif
500: }

502: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
503: {
504:   /***************************************************************************/
505:   /* Set work vectors needed by conjugate gradient method and allocate       */
506:   /***************************************************************************/

508:   PetscFunctionBegin;
509:   PetscCall(KSPSetWorkVecs(ksp, 3));
510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }

513: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
514: {
515:   PetscFunctionBegin;
516:   /***************************************************************************/
517:   /* Clear composed functions                                                */
518:   /***************************************************************************/

520:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
521:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
522:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));

524:   /***************************************************************************/
525:   /* Destroy KSP object.                                                     */
526:   /***************************************************************************/

528:   PetscCall(KSPDestroyDefault(ksp));
529:   PetscFunctionReturn(PETSC_SUCCESS);
530: }

532: static PetscErrorCode KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
533: {
534:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;

536:   PetscFunctionBegin;
537:   cg->radius = radius;
538:   PetscFunctionReturn(PETSC_SUCCESS);
539: }

541: static PetscErrorCode KSPCGGetNormD_NASH(KSP ksp, PetscReal *norm_d)
542: {
543:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;

545:   PetscFunctionBegin;
546:   *norm_d = cg->norm_d;
547:   PetscFunctionReturn(PETSC_SUCCESS);
548: }

550: static PetscErrorCode KSPCGGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
551: {
552:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;

554:   PetscFunctionBegin;
555:   *o_fcn = cg->o_fcn;
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }

559: static PetscErrorCode KSPCGSetFromOptions_NASH(KSP ksp, PetscOptionItems *PetscOptionsObject)
560: {
561:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;

563:   PetscFunctionBegin;
564:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG NASH options");

566:   PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));

568:   PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));

570:   PetscOptionsHeadEnd();
571:   PetscFunctionReturn(PETSC_SUCCESS);
572: }

574: /*MC
575:      KSPNASH -   Code to run conjugate gradient method subject to a constraint on the solution norm.

577:    Options Database Keys:
578: .      -ksp_cg_radius <r> - Trust Region Radius

580:    Level: developer

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

585:   Uses preconditioned conjugate gradient to compute
586:   an approximate minimizer of the quadratic function

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

590:    subject to the trust region constraint

592:             || s || <= delta,

594:    where

596:      delta is the trust region radius,
597:      g is the gradient vector,
598:      H is the Hessian approximation, and
599:      M is the positive definite preconditioner matrix.

601:    `KSPConvergedReason` may be
602: .vb
603:   KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
604:   KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step,
605: .ve
606:   other `KSP` converged/diverged reasons

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

610:   Reference:
611:    Nash, Stephen G. Newton-type minimization via the Lanczos method. SIAM Journal on Numerical Analysis 21, no. 4 (1984): 770-788.

613: .seealso: [](ch_ksp), `KSPQCG`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`
614: M*/

616: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
617: {
618:   KSPCG_NASH *cg;

620:   PetscFunctionBegin;
621:   PetscCall(PetscNew(&cg));
622:   cg->radius = 0.0;
623:   cg->dtype  = NASH_UNPRECONDITIONED_DIRECTION;

625:   ksp->data = (void *)cg;
626:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
627:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
628:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
629:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
630:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));

632:   /***************************************************************************/
633:   /* Sets the functions that are associated with this data structure         */
634:   /* (in C++ this is the same as defining virtual functions).                */
635:   /***************************************************************************/

637:   ksp->ops->setup          = KSPCGSetUp_NASH;
638:   ksp->ops->solve          = KSPCGSolve_NASH;
639:   ksp->ops->destroy        = KSPCGDestroy_NASH;
640:   ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
641:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
642:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
643:   ksp->ops->view           = NULL;

645:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_NASH));
646:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_NASH));
647:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_NASH));
648:   PetscFunctionReturn(PETSC_SUCCESS);
649: }