Actual source code: nash.c
1: #include <../src/ksp/ksp/impls/cg/nash/nashimpl.h>
3: #define NASH_PRECONDITIONED_DIRECTION 0
4: #define NASH_UNPRECONDITIONED_DIRECTION 1
5: #define NASH_DIRECTION_TYPES 2
7: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
9: static PetscErrorCode KSPCGSolve_NASH(KSP ksp)
10: {
11: #if PetscDefined(USE_COMPLEX)
12: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "NASH is not available for complex systems");
13: #else
14: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
15: Mat Qmat, Mmat;
16: Vec r, z, p, d;
17: 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;
23: PetscInt max_cg_its;
25: PetscFunctionBegin;
26: /***************************************************************************/
27: /* Check the arguments and parameters. */
28: /***************************************************************************/
30: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
32: /***************************************************************************/
33: /* Get the workspace vectors and initialize variables */
34: /***************************************************************************/
36: r2 = cg->radius * cg->radius;
37: r = ksp->work[0];
38: z = ksp->work[1];
39: p = ksp->work[2];
40: d = ksp->vec_sol;
41: pc = ksp->pc;
43: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
45: PetscCall(VecGetSize(d, &max_cg_its));
46: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
47: ksp->its = 0;
49: /***************************************************************************/
50: /* Initialize objective function and direction. */
51: /***************************************************************************/
53: cg->o_fcn = 0.0;
55: PetscCall(VecSet(d, 0.0)); /* d = 0 */
56: cg->norm_d = 0.0;
58: /***************************************************************************/
59: /* Begin the conjugate gradient method. Check the right-hand side for */
60: /* numerical problems. The check for not-a-number and infinite values */
61: /* need be performed only once. */
62: /***************************************************************************/
64: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
65: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
66: KSPCheckDot(ksp, rr);
68: /***************************************************************************/
69: /* Check the preconditioner for numerical problems and for positive */
70: /* definiteness. The check for not-a-number and infinite values need be */
71: /* performed only once. */
72: /***************************************************************************/
74: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
75: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
76: if (PetscIsInfOrNanScalar(rz)) {
77: /*************************************************************************/
78: /* The preconditioner contains not-a-number or an infinite value. */
79: /* Return the gradient direction intersected with the trust region. */
80: /*************************************************************************/
82: ksp->reason = KSP_DIVERGED_NANORINF;
83: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: bad preconditioner: rz=%g\n", (double)rz));
85: if (cg->radius) {
86: if (r2 >= rr) {
87: alpha = 1.0;
88: cg->norm_d = PetscSqrtReal(rr);
89: } else {
90: alpha = PetscSqrtReal(r2 / rr);
91: cg->norm_d = cg->radius;
92: }
94: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
96: /***********************************************************************/
97: /* Compute objective function. */
98: /***********************************************************************/
100: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
101: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
102: PetscCall(VecDot(d, z, &cg->o_fcn));
103: cg->o_fcn = -cg->o_fcn;
104: ++ksp->its;
105: }
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: if (rz < 0.0) {
110: /*************************************************************************/
111: /* The preconditioner is indefinite. Because this is the first */
112: /* and we do not have a direction yet, we use the gradient step. Note */
113: /* that we cannot use the preconditioned norm when computing the step */
114: /* because the matrix is indefinite. */
115: /*************************************************************************/
117: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
118: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz));
120: if (cg->radius) {
121: if (r2 >= rr) {
122: alpha = 1.0;
123: cg->norm_d = PetscSqrtReal(rr);
124: } else {
125: alpha = PetscSqrtReal(r2 / rr);
126: cg->norm_d = cg->radius;
127: }
129: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
131: /***********************************************************************/
132: /* Compute objective function. */
133: /***********************************************************************/
135: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
136: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
137: PetscCall(VecDot(d, z, &cg->o_fcn));
138: cg->o_fcn = -cg->o_fcn;
139: ++ksp->its;
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /***************************************************************************/
145: /* As far as we know, the preconditioner is positive semidefinite. */
146: /* Compute and log the residual. Check convergence because this */
147: /* initializes things, but do not terminate until at least one conjugate */
148: /* gradient iteration has been performed. */
149: /***************************************************************************/
151: switch (ksp->normtype) {
152: case KSP_NORM_PRECONDITIONED:
153: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
154: break;
156: case KSP_NORM_UNPRECONDITIONED:
157: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
158: break;
160: case KSP_NORM_NATURAL:
161: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
162: break;
164: default:
165: norm_r = 0.0;
166: break;
167: }
169: PetscCall(KSPLogResidualHistory(ksp, norm_r));
170: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
171: ksp->rnorm = norm_r;
173: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
175: /***************************************************************************/
176: /* Compute the first direction and update the iteration. */
177: /***************************************************************************/
179: PetscCall(VecCopy(z, p)); /* p = z */
180: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
181: ++ksp->its;
183: /***************************************************************************/
184: /* Check the matrix for numerical problems. */
185: /***************************************************************************/
187: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
188: if (PetscIsInfOrNanScalar(kappa)) {
189: /*************************************************************************/
190: /* The matrix produced not-a-number or an infinite value. In this case, */
191: /* we must stop and use the gradient direction. This condition need */
192: /* only be checked once. */
193: /*************************************************************************/
195: ksp->reason = KSP_DIVERGED_NANORINF;
196: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: bad matrix: kappa=%g\n", (double)kappa));
198: if (cg->radius) {
199: if (r2 >= rr) {
200: alpha = 1.0;
201: cg->norm_d = PetscSqrtReal(rr);
202: } else {
203: alpha = PetscSqrtReal(r2 / rr);
204: cg->norm_d = cg->radius;
205: }
207: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
209: /***********************************************************************/
210: /* Compute objective function. */
211: /***********************************************************************/
213: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
214: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
215: PetscCall(VecDot(d, z, &cg->o_fcn));
216: cg->o_fcn = -cg->o_fcn;
217: ++ksp->its;
218: }
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /***************************************************************************/
223: /* Initialize variables for calculating the norm of the direction. */
224: /***************************************************************************/
226: dMp = 0.0;
227: norm_d = 0.0;
228: switch (cg->dtype) {
229: case NASH_PRECONDITIONED_DIRECTION:
230: norm_p = rz;
231: break;
233: default:
234: PetscCall(VecDot(p, p, &norm_p));
235: break;
236: }
238: /***************************************************************************/
239: /* Check for negative curvature. */
240: /***************************************************************************/
242: if (kappa <= 0.0) {
243: /*************************************************************************/
244: /* In this case, the matrix is indefinite and we have encountered a */
245: /* direction of negative curvature. Because negative curvature occurs */
246: /* during the first step, we must follow a direction. */
247: /*************************************************************************/
249: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
250: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa));
252: if (cg->radius && norm_p > 0.0) {
253: /***********************************************************************/
254: /* Follow direction of negative curvature to the boundary of the */
255: /* trust region. */
256: /***********************************************************************/
258: step = PetscSqrtReal(r2 / norm_p);
259: cg->norm_d = cg->radius;
261: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
263: /***********************************************************************/
264: /* Update objective function. */
265: /***********************************************************************/
267: cg->o_fcn += step * (0.5 * step * kappa - rz);
268: } else if (cg->radius) {
269: /***********************************************************************/
270: /* The norm of the preconditioned direction is zero; use the gradient */
271: /* step. */
272: /***********************************************************************/
274: if (r2 >= rr) {
275: alpha = 1.0;
276: cg->norm_d = PetscSqrtReal(rr);
277: } else {
278: alpha = PetscSqrtReal(r2 / rr);
279: cg->norm_d = cg->radius;
280: }
282: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
284: /***********************************************************************/
285: /* Compute objective function. */
286: /***********************************************************************/
288: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
289: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
290: PetscCall(VecDot(d, z, &cg->o_fcn));
291: cg->o_fcn = -cg->o_fcn;
292: ++ksp->its;
293: }
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /***************************************************************************/
298: /* Run the conjugate gradient method until either the problem is solved, */
299: /* we encounter the boundary of the trust region, or the conjugate */
300: /* gradient method breaks down. */
301: /***************************************************************************/
303: while (1) {
304: /*************************************************************************/
305: /* Know that kappa is nonzero, because we have not broken down, so we */
306: /* can compute the steplength. */
307: /*************************************************************************/
309: alpha = rz / kappa;
311: /*************************************************************************/
312: /* Compute the steplength and check for intersection with the trust */
313: /* region. */
314: /*************************************************************************/
316: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
317: if (cg->radius && norm_dp1 >= r2) {
318: /***********************************************************************/
319: /* In this case, the matrix is positive definite as far as we know. */
320: /* However, the full step goes beyond the trust region. */
321: /***********************************************************************/
323: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
324: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: constrained step: radius=%g\n", (double)cg->radius));
326: if (norm_p > 0.0) {
327: /*********************************************************************/
328: /* Follow the direction to the boundary of the trust region. */
329: /*********************************************************************/
331: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
332: cg->norm_d = cg->radius;
334: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
336: /*********************************************************************/
337: /* Update objective function. */
338: /*********************************************************************/
340: cg->o_fcn += step * (0.5 * step * kappa - rz);
341: } else {
342: /*********************************************************************/
343: /* The norm of the direction is zero; there is nothing to follow. */
344: /*********************************************************************/
345: }
346: break;
347: }
349: /*************************************************************************/
350: /* Now we can update the direction and residual. */
351: /*************************************************************************/
353: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
354: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
355: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
357: switch (cg->dtype) {
358: case NASH_PRECONDITIONED_DIRECTION:
359: norm_d = norm_dp1;
360: break;
362: default:
363: PetscCall(VecDot(d, d, &norm_d));
364: break;
365: }
366: cg->norm_d = PetscSqrtReal(norm_d);
368: /*************************************************************************/
369: /* Update objective function. */
370: /*************************************************************************/
372: cg->o_fcn -= 0.5 * alpha * rz;
374: /*************************************************************************/
375: /* Check that the preconditioner appears positive semidefinite. */
376: /*************************************************************************/
378: rzm1 = rz;
379: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
380: if (rz < 0.0) {
381: /***********************************************************************/
382: /* The preconditioner is indefinite. */
383: /***********************************************************************/
385: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
386: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz));
387: break;
388: }
390: /*************************************************************************/
391: /* As far as we know, the preconditioner is positive semidefinite. */
392: /* Compute the residual and check for convergence. */
393: /*************************************************************************/
395: switch (ksp->normtype) {
396: case KSP_NORM_PRECONDITIONED:
397: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
398: break;
400: case KSP_NORM_UNPRECONDITIONED:
401: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
402: break;
404: case KSP_NORM_NATURAL:
405: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
406: break;
408: default:
409: norm_r = 0.;
410: break;
411: }
413: PetscCall(KSPLogResidualHistory(ksp, norm_r));
414: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
415: ksp->rnorm = norm_r;
417: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
418: if (ksp->reason) {
419: /***********************************************************************/
420: /* The method has converged. */
421: /***********************************************************************/
423: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
424: break;
425: }
427: /*************************************************************************/
428: /* We have not converged yet. Check for breakdown. */
429: /*************************************************************************/
431: beta = rz / rzm1;
432: if (PetscAbsReal(beta) <= 0.0) {
433: /***********************************************************************/
434: /* Conjugate gradients has broken down. */
435: /***********************************************************************/
437: ksp->reason = KSP_DIVERGED_BREAKDOWN;
438: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: breakdown: beta=%g\n", (double)beta));
439: break;
440: }
442: /*************************************************************************/
443: /* Check iteration limit. */
444: /*************************************************************************/
446: if (ksp->its >= max_cg_its) {
447: ksp->reason = KSP_DIVERGED_ITS;
448: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
449: break;
450: }
452: /*************************************************************************/
453: /* Update p and the norms. */
454: /*************************************************************************/
456: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
458: switch (cg->dtype) {
459: case NASH_PRECONDITIONED_DIRECTION:
460: dMp = beta * (dMp + alpha * norm_p);
461: norm_p = beta * (rzm1 + beta * norm_p);
462: break;
464: default:
465: PetscCall(VecDot(d, p, &dMp));
466: PetscCall(VecDot(p, p, &norm_p));
467: break;
468: }
470: /*************************************************************************/
471: /* Compute the new direction and update the iteration. */
472: /*************************************************************************/
474: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
475: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
476: ++ksp->its;
478: /*************************************************************************/
479: /* Check for negative curvature. */
480: /*************************************************************************/
482: if (kappa <= 0.0) {
483: /***********************************************************************/
484: /* In this case, the matrix is indefinite and we have encountered */
485: /* a direction of negative curvature. Stop at the base. */
486: /***********************************************************************/
488: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
489: PetscCall(PetscInfo(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa));
490: break;
491: }
492: }
493: PetscFunctionReturn(PETSC_SUCCESS);
494: #endif
495: }
497: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
498: {
499: /***************************************************************************/
500: /* Set work vectors needed by conjugate gradient method and allocate */
501: /***************************************************************************/
503: PetscFunctionBegin;
504: PetscCall(KSPSetWorkVecs(ksp, 3));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
509: {
510: PetscFunctionBegin;
511: /***************************************************************************/
512: /* Clear composed functions */
513: /***************************************************************************/
515: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
516: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
517: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
519: /***************************************************************************/
520: /* Destroy KSP object. */
521: /***************************************************************************/
523: PetscCall(KSPDestroyDefault(ksp));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
528: {
529: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
531: PetscFunctionBegin;
532: cg->radius = radius;
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: static PetscErrorCode KSPCGGetNormD_NASH(KSP ksp, PetscReal *norm_d)
537: {
538: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
540: PetscFunctionBegin;
541: *norm_d = cg->norm_d;
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: static PetscErrorCode KSPCGGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
546: {
547: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
549: PetscFunctionBegin;
550: *o_fcn = cg->o_fcn;
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: static PetscErrorCode KSPCGSetFromOptions_NASH(KSP ksp, PetscOptionItems PetscOptionsObject)
555: {
556: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
558: PetscFunctionBegin;
559: PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG NASH options");
561: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
563: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
565: PetscOptionsHeadEnd();
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: /*MC
570: KSPNASH - Code to run conjugate gradient method subject to a constraint on the solution norm in a trust region method {cite}`nash1984newton`
572: Options Database Keys:
573: . -ksp_cg_radius radius - Trust Region Radius
575: Level: developer
577: Notes:
578: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
580: Uses preconditioned conjugate gradient to compute
581: an approximate minimizer of the quadratic function
583: $$
584: q(s) = g^T * s + 0.5 * s^T * H * s
585: $$
587: subject to the trust region constraint
589: $$
590: || s || \le delta,
591: $$
593: where
594: .vb
595: delta is the trust region radius,
596: g is the gradient vector,
597: H is the Hessian approximation, and
598: .ve
600: `KSPConvergedReason` may include
601: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
602: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,
604: The preconditioner supplied must be symmetric and positive definite.
606: .seealso: [](ch_ksp), `KSPQCG`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`
607: M*/
609: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
610: {
611: KSPCG_NASH *cg;
613: PetscFunctionBegin;
614: PetscCall(PetscNew(&cg));
615: cg->radius = 0.0;
616: cg->dtype = NASH_UNPRECONDITIONED_DIRECTION;
618: ksp->data = (void *)cg;
619: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
620: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
621: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
622: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
623: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
625: /***************************************************************************/
626: /* Sets the functions that are associated with this data structure */
627: /* (in C++ this is the same as defining virtual functions). */
628: /***************************************************************************/
630: ksp->ops->setup = KSPCGSetUp_NASH;
631: ksp->ops->solve = KSPCGSolve_NASH;
632: ksp->ops->destroy = KSPCGDestroy_NASH;
633: ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
634: ksp->ops->buildsolution = KSPBuildSolutionDefault;
635: ksp->ops->buildresidual = KSPBuildResidualDefault;
636: ksp->ops->view = NULL;
638: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_NASH));
639: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_NASH));
640: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_NASH));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }