Actual source code: nash.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/
3: #include <../src/ksp/ksp/impls/cg/nash/nashimpl.h>
5: #define NASH_PRECONDITIONED_DIRECTION 0
6: #define NASH_UNPRECONDITIONED_DIRECTION 1
7: #define NASH_DIRECTION_TYPES 2
9: static const char *DType_Table[64] = { "preconditioned", "unpreconditioned"};
13: /*@
14: KSPNASHSetRadius - Sets the radius of the trust region.
16: Logically Collective on KSP
18: Input Parameters:
19: + ksp - the iterative context
20: - radius - the trust region radius (Infinity is the default)
22: Options Database Key:
23: . -ksp_nash_radius <r>
25: Level: advanced
27: .keywords: KSP, NASH, set, trust region radius
28: @*/
29: PetscErrorCode KSPNASHSetRadius(KSP ksp, PetscReal radius)
30: {
35: if (radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Radius negative");
37: PetscTryMethod(ksp,"KSPNASHSetRadius_C",(KSP,PetscReal),(ksp,radius));
38: return(0);
39: }
43: /*@
44: KSPNASHGetNormD - Got norm of the direction.
46: Collective on KSP
48: Input Parameters:
49: + ksp - the iterative context
50: - norm_d - the norm of the direction
52: Level: advanced
54: .keywords: KSP, NASH, get, norm direction
55: @*/
56: PetscErrorCode KSPNASHGetNormD(KSP ksp, PetscReal *norm_d)
57: {
62: PetscUseMethod(ksp,"KSPNASHGetNormD_C",(KSP,PetscReal*),(ksp,norm_d));
63: return(0);
64: }
68: /*@
69: KSPNASHGetObjFcn - Get objective function value.
71: Collective on KSP
73: Input Parameters:
74: + ksp - the iterative context
75: - o_fcn - the objective function value
77: Level: advanced
79: .keywords: KSP, NASH, get, objective function
80: @*/
81: PetscErrorCode KSPNASHGetObjFcn(KSP ksp, PetscReal *o_fcn)
82: {
87: PetscUseMethod(ksp,"KSPNASHGetObjFcn_C",(KSP,PetscReal*),(ksp,o_fcn));
88: return(0);
89: }
94: PetscErrorCode KSPSolve_NASH(KSP ksp)
95: {
96: #if defined(PETSC_USE_COMPLEX)
97: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "NASH is not available for complex systems");
98: #else
99: KSP_NASH *cg = (KSP_NASH*)ksp->data;
101: MatStructure pflag;
102: Mat Qmat, Mmat;
103: Vec r, z, p, d;
104: PC pc;
106: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
107: PetscReal alpha, beta, kappa, rz, rzm1;
108: PetscReal rr, r2, step;
110: PetscInt max_cg_its;
112: PetscBool diagonalscale;
115: /***************************************************************************/
116: /* Check the arguments and parameters. */
117: /***************************************************************************/
119: PCGetDiagonalScale(ksp->pc, &diagonalscale);
120: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
121: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
123: /***************************************************************************/
124: /* Get the workspace vectors and initialize variables */
125: /***************************************************************************/
127: r2 = cg->radius * cg->radius;
128: r = ksp->work[0];
129: z = ksp->work[1];
130: p = ksp->work[2];
131: d = ksp->vec_sol;
132: pc = ksp->pc;
134: PCGetOperators(pc, &Qmat, &Mmat, &pflag);
136: VecGetSize(d, &max_cg_its);
137: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
138: ksp->its = 0;
140: /***************************************************************************/
141: /* Initialize objective function and direction. */
142: /***************************************************************************/
144: cg->o_fcn = 0.0;
146: VecSet(d, 0.0); /* d = 0 */
147: cg->norm_d = 0.0;
149: /***************************************************************************/
150: /* Begin the conjugate gradient method. Check the right-hand side for */
151: /* numerical problems. The check for not-a-number and infinite values */
152: /* need be performed only once. */
153: /***************************************************************************/
155: VecCopy(ksp->vec_rhs, r); /* r = -grad */
156: VecDot(r, r, &rr); /* rr = r^T r */
157: if (PetscIsInfOrNanScalar(rr)) {
158: /*************************************************************************/
159: /* The right-hand side contains not-a-number or an infinite value. */
160: /* The gradient step does not work; return a zero value for the step. */
161: /*************************************************************************/
163: ksp->reason = KSP_DIVERGED_NANORINF;
164: PetscInfo1(ksp, "KSPSolve_NASH: bad right-hand side: rr=%g\n", rr);
165: return(0);
166: }
168: /***************************************************************************/
169: /* Check the preconditioner for numerical problems and for positive */
170: /* definiteness. The check for not-a-number and infinite values need be */
171: /* performed only once. */
172: /***************************************************************************/
174: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
175: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
176: if (PetscIsInfOrNanScalar(rz)) {
177: /*************************************************************************/
178: /* The preconditioner contains not-a-number or an infinite value. */
179: /* Return the gradient direction intersected with the trust region. */
180: /*************************************************************************/
182: ksp->reason = KSP_DIVERGED_NANORINF;
183: PetscInfo1(ksp, "KSPSolve_NASH: bad preconditioner: rz=%g\n", rz);
185: if (cg->radius) {
186: if (r2 >= rr) {
187: alpha = 1.0;
188: cg->norm_d = PetscSqrtReal(rr);
189: } else {
190: alpha = PetscSqrtReal(r2 / rr);
191: cg->norm_d = cg->radius;
192: }
194: VecAXPY(d, alpha, r); /* d = d + alpha r */
196: /***********************************************************************/
197: /* Compute objective function. */
198: /***********************************************************************/
200: KSP_MatMult(ksp, Qmat, d, z);
201: VecAYPX(z, -0.5, ksp->vec_rhs);
202: VecDot(d, z, &cg->o_fcn);
203: cg->o_fcn = -cg->o_fcn;
204: ++ksp->its;
205: }
206: return(0);
207: }
209: if (rz < 0.0) {
210: /*************************************************************************/
211: /* The preconditioner is indefinite. Because this is the first */
212: /* and we do not have a direction yet, we use the gradient step. Note */
213: /* that we cannot use the preconditioned norm when computing the step */
214: /* because the matrix is indefinite. */
215: /*************************************************************************/
217: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
218: PetscInfo1(ksp, "KSPSolve_NASH: indefinite preconditioner: rz=%g\n", rz);
220: if (cg->radius) {
221: if (r2 >= rr) {
222: alpha = 1.0;
223: cg->norm_d = PetscSqrtReal(rr);
224: } else {
225: alpha = PetscSqrtReal(r2 / rr);
226: cg->norm_d = cg->radius;
227: }
229: VecAXPY(d, alpha, r); /* d = d + alpha r */
231: /***********************************************************************/
232: /* Compute objective function. */
233: /***********************************************************************/
235: KSP_MatMult(ksp, Qmat, d, z);
236: VecAYPX(z, -0.5, ksp->vec_rhs);
237: VecDot(d, z, &cg->o_fcn);
238: cg->o_fcn = -cg->o_fcn;
239: ++ksp->its;
240: }
241: return(0);
242: }
244: /***************************************************************************/
245: /* As far as we know, the preconditioner is positive semidefinite. */
246: /* Compute and log the residual. Check convergence because this */
247: /* initializes things, but do not terminate until at least one conjugate */
248: /* gradient iteration has been performed. */
249: /***************************************************************************/
251: switch (ksp->normtype) {
252: case KSP_NORM_PRECONDITIONED:
253: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
254: break;
256: case KSP_NORM_UNPRECONDITIONED:
257: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
258: break;
260: case KSP_NORM_NATURAL:
261: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
262: break;
264: default:
265: norm_r = 0.0;
266: break;
267: }
269: KSPLogResidualHistory(ksp, norm_r);
270: KSPMonitor(ksp, ksp->its, norm_r);
271: ksp->rnorm = norm_r;
273: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
275: /***************************************************************************/
276: /* Compute the first direction and update the iteration. */
277: /***************************************************************************/
279: VecCopy(z, p); /* p = z */
280: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
281: ++ksp->its;
283: /***************************************************************************/
284: /* Check the matrix for numerical problems. */
285: /***************************************************************************/
287: VecDot(p, z, &kappa); /* kappa = p^T Q p */
288: if (PetscIsInfOrNanScalar(kappa)) {
289: /*************************************************************************/
290: /* The matrix produced not-a-number or an infinite value. In this case, */
291: /* we must stop and use the gradient direction. This condition need */
292: /* only be checked once. */
293: /*************************************************************************/
295: ksp->reason = KSP_DIVERGED_NANORINF;
296: PetscInfo1(ksp, "KSPSolve_NASH: bad matrix: kappa=%g\n", kappa);
298: if (cg->radius) {
299: if (r2 >= rr) {
300: alpha = 1.0;
301: cg->norm_d = PetscSqrtReal(rr);
302: } else {
303: alpha = PetscSqrtReal(r2 / rr);
304: cg->norm_d = cg->radius;
305: }
307: VecAXPY(d, alpha, r); /* d = d + alpha r */
309: /***********************************************************************/
310: /* Compute objective function. */
311: /***********************************************************************/
313: KSP_MatMult(ksp, Qmat, d, z);
314: VecAYPX(z, -0.5, ksp->vec_rhs);
315: VecDot(d, z, &cg->o_fcn);
316: cg->o_fcn = -cg->o_fcn;
317: ++ksp->its;
318: }
319: return(0);
320: }
322: /***************************************************************************/
323: /* Initialize variables for calculating the norm of the direction. */
324: /***************************************************************************/
326: dMp = 0.0;
327: norm_d = 0.0;
328: switch (cg->dtype) {
329: case NASH_PRECONDITIONED_DIRECTION:
330: norm_p = rz;
331: break;
333: default:
334: VecDot(p, p, &norm_p);
335: break;
336: }
338: /***************************************************************************/
339: /* Check for negative curvature. */
340: /***************************************************************************/
342: if (kappa <= 0.0) {
343: /*************************************************************************/
344: /* In this case, the matrix is indefinite and we have encountered a */
345: /* direction of negative curvature. Because negative curvature occurs */
346: /* during the first step, we must follow a direction. */
347: /*************************************************************************/
349: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
350: PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", kappa);
352: if (cg->radius && norm_p > 0.0) {
353: /***********************************************************************/
354: /* Follow direction of negative curvature to the boundary of the */
355: /* trust region. */
356: /***********************************************************************/
358: step = PetscSqrtReal(r2 / norm_p);
359: cg->norm_d = cg->radius;
361: VecAXPY(d, step, p); /* d = d + step p */
363: /***********************************************************************/
364: /* Update objective function. */
365: /***********************************************************************/
367: cg->o_fcn += step * (0.5 * step * kappa - rz);
368: } else if (cg->radius) {
369: /***********************************************************************/
370: /* The norm of the preconditioned direction is zero; use the gradient */
371: /* step. */
372: /***********************************************************************/
374: if (r2 >= rr) {
375: alpha = 1.0;
376: cg->norm_d = PetscSqrtReal(rr);
377: } else {
378: alpha = PetscSqrtReal(r2 / rr);
379: cg->norm_d = cg->radius;
380: }
382: VecAXPY(d, alpha, r); /* d = d + alpha r */
384: /***********************************************************************/
385: /* Compute objective function. */
386: /***********************************************************************/
388: KSP_MatMult(ksp, Qmat, d, z);
389: VecAYPX(z, -0.5, ksp->vec_rhs);
390: VecDot(d, z, &cg->o_fcn);
391: cg->o_fcn = -cg->o_fcn;
392: ++ksp->its;
393: }
394: return(0);
395: }
397: /***************************************************************************/
398: /* Run the conjugate gradient method until either the problem is solved, */
399: /* we encounter the boundary of the trust region, or the conjugate */
400: /* gradient method breaks down. */
401: /***************************************************************************/
403: while (1) {
404: /*************************************************************************/
405: /* Know that kappa is nonzero, because we have not broken down, so we */
406: /* can compute the steplength. */
407: /*************************************************************************/
409: alpha = rz / kappa;
411: /*************************************************************************/
412: /* Compute the steplength and check for intersection with the trust */
413: /* region. */
414: /*************************************************************************/
416: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
417: if (cg->radius && norm_dp1 >= r2) {
418: /***********************************************************************/
419: /* In this case, the matrix is positive definite as far as we know. */
420: /* However, the full step goes beyond the trust region. */
421: /***********************************************************************/
423: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
424: PetscInfo1(ksp, "KSPSolve_NASH: constrained step: radius=%g\n", cg->radius);
426: if (norm_p > 0.0) {
427: /*********************************************************************/
428: /* Follow the direction to the boundary of the trust region. */
429: /*********************************************************************/
431: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
432: cg->norm_d = cg->radius;
434: VecAXPY(d, step, p); /* d = d + step p */
436: /*********************************************************************/
437: /* Update objective function. */
438: /*********************************************************************/
440: cg->o_fcn += step * (0.5 * step * kappa - rz);
441: } else {
442: /*********************************************************************/
443: /* The norm of the direction is zero; there is nothing to follow. */
444: /*********************************************************************/
445: }
446: break;
447: }
449: /*************************************************************************/
450: /* Now we can update the direction and residual. */
451: /*************************************************************************/
453: VecAXPY(d, alpha, p); /* d = d + alpha p */
454: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
455: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
457: switch (cg->dtype) {
458: case NASH_PRECONDITIONED_DIRECTION:
459: norm_d = norm_dp1;
460: break;
462: default:
463: VecDot(d, d, &norm_d);
464: break;
465: }
466: cg->norm_d = PetscSqrtReal(norm_d);
468: /*************************************************************************/
469: /* Update objective function. */
470: /*************************************************************************/
472: cg->o_fcn -= 0.5 * alpha * rz;
474: /*************************************************************************/
475: /* Check that the preconditioner appears positive semidefinite. */
476: /*************************************************************************/
478: rzm1 = rz;
479: VecDot(r, z, &rz); /* rz = r^T z */
480: if (rz < 0.0) {
481: /***********************************************************************/
482: /* The preconditioner is indefinite. */
483: /***********************************************************************/
485: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
486: PetscInfo1(ksp, "KSPSolve_NASH: cg indefinite preconditioner: rz=%g\n", rz);
487: break;
488: }
490: /*************************************************************************/
491: /* As far as we know, the preconditioner is positive semidefinite. */
492: /* Compute the residual and check for convergence. */
493: /*************************************************************************/
495: switch (ksp->normtype) {
496: case KSP_NORM_PRECONDITIONED:
497: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
498: break;
500: case KSP_NORM_UNPRECONDITIONED:
501: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
502: break;
504: case KSP_NORM_NATURAL:
505: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
506: break;
508: default:
509: norm_r = 0.;
510: break;
511: }
513: KSPLogResidualHistory(ksp, norm_r);
514: KSPMonitor(ksp, ksp->its, norm_r);
515: ksp->rnorm = norm_r;
517: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
518: if (ksp->reason) {
519: /***********************************************************************/
520: /* The method has converged. */
521: /***********************************************************************/
523: PetscInfo2(ksp, "KSPSolve_NASH: truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
524: break;
525: }
527: /*************************************************************************/
528: /* We have not converged yet. Check for breakdown. */
529: /*************************************************************************/
531: beta = rz / rzm1;
532: if (fabs(beta) <= 0.0) {
533: /***********************************************************************/
534: /* Conjugate gradients has broken down. */
535: /***********************************************************************/
537: ksp->reason = KSP_DIVERGED_BREAKDOWN;
538: PetscInfo1(ksp, "KSPSolve_NASH: breakdown: beta=%g\n", beta);
539: break;
540: }
542: /*************************************************************************/
543: /* Check iteration limit. */
544: /*************************************************************************/
546: if (ksp->its >= max_cg_its) {
547: ksp->reason = KSP_DIVERGED_ITS;
548: PetscInfo1(ksp, "KSPSolve_NASH: iterlim: its=%d\n", ksp->its);
549: break;
550: }
552: /*************************************************************************/
553: /* Update p and the norms. */
554: /*************************************************************************/
556: VecAYPX(p, beta, z); /* p = z + beta p */
558: switch (cg->dtype) {
559: case NASH_PRECONDITIONED_DIRECTION:
560: dMp = beta*(dMp + alpha*norm_p);
561: norm_p = beta*(rzm1 + beta*norm_p);
562: break;
564: default:
565: VecDot(d, p, &dMp);
566: VecDot(p, p, &norm_p);
567: break;
568: }
570: /*************************************************************************/
571: /* Compute the new direction and update the iteration. */
572: /*************************************************************************/
574: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
575: VecDot(p, z, &kappa); /* kappa = p^T Q p */
576: ++ksp->its;
578: /*************************************************************************/
579: /* Check for negative curvature. */
580: /*************************************************************************/
582: if (kappa <= 0.0) {
583: /***********************************************************************/
584: /* In this case, the matrix is indefinite and we have encountered */
585: /* a direction of negative curvature. Stop at the base. */
586: /***********************************************************************/
588: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
589: PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", kappa);
590: break;
591: }
592: }
593: return(0);
594: #endif
595: }
599: PetscErrorCode KSPSetUp_NASH(KSP ksp)
600: {
604: /***************************************************************************/
605: /* Set work vectors needed by conjugate gradient method and allocate */
606: /***************************************************************************/
608: KSPSetWorkVecs(ksp, 3);
609: return(0);
610: }
614: PetscErrorCode KSPDestroy_NASH(KSP ksp)
615: {
619: /***************************************************************************/
620: /* Clear composed functions */
621: /***************************************************************************/
623: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",NULL);
624: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",NULL);
625: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",NULL);
627: /***************************************************************************/
628: /* Destroy KSP object. */
629: /***************************************************************************/
631: KSPDestroyDefault(ksp);
632: return(0);
633: }
637: static PetscErrorCode KSPNASHSetRadius_NASH(KSP ksp, PetscReal radius)
638: {
639: KSP_NASH *cg = (KSP_NASH*)ksp->data;
642: cg->radius = radius;
643: return(0);
644: }
648: static PetscErrorCode KSPNASHGetNormD_NASH(KSP ksp, PetscReal *norm_d)
649: {
650: KSP_NASH *cg = (KSP_NASH*)ksp->data;
653: *norm_d = cg->norm_d;
654: return(0);
655: }
659: static PetscErrorCode KSPNASHGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
660: {
661: KSP_NASH *cg = (KSP_NASH*)ksp->data;
664: *o_fcn = cg->o_fcn;
665: return(0);
666: }
670: PetscErrorCode KSPSetFromOptions_NASH(KSP ksp)
671: {
673: KSP_NASH *cg = (KSP_NASH*)ksp->data;
676: PetscOptionsHead("KSP NASH options");
678: PetscOptionsReal("-ksp_nash_radius", "Trust Region Radius", "KSPNASHSetRadius", cg->radius, &cg->radius, NULL);
680: PetscOptionsEList("-ksp_nash_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
682: PetscOptionsTail();
683: return(0);
684: }
686: /*MC
687: KSPNASH - Code to run conjugate gradient method subject to a constraint
688: on the solution norm. This is used in Trust Region methods for
689: nonlinear equations, SNESNEWTONTR
691: Options Database Keys:
692: . -ksp_nash_radius <r> - Trust Region Radius
694: Notes: This is rarely used directly
696: Level: developer
698: Use preconditioned conjugate gradient to compute
699: an approximate minimizer of the quadratic function
701: q(s) = g^T * s + 0.5 * s^T * H * s
703: subject to the trust region constraint
705: || s || <= delta,
707: where
709: delta is the trust region radius,
710: g is the gradient vector,
711: H is the Hessian approximation, and
712: M is the positive definite preconditioner matrix.
714: KSPConvergedReason may be
715: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
716: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
717: $ other KSP converged/diverged reasons
719: Notes:
720: The preconditioner supplied should be symmetric and positive definite.
722: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPNASHSetRadius(), KSPNASHGetNormD(), KSPNASHGetObjFcn()
723: M*/
727: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
728: {
730: KSP_NASH *cg;
733: PetscNewLog(ksp,KSP_NASH, &cg);
734: cg->radius = 0.0;
735: cg->dtype = NASH_UNPRECONDITIONED_DIRECTION;
737: ksp->data = (void*) cg;
738: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
739: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
740: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
742: /***************************************************************************/
743: /* Sets the functions that are associated with this data structure */
744: /* (in C++ this is the same as defining virtual functions). */
745: /***************************************************************************/
747: ksp->ops->setup = KSPSetUp_NASH;
748: ksp->ops->solve = KSPSolve_NASH;
749: ksp->ops->destroy = KSPDestroy_NASH;
750: ksp->ops->setfromoptions = KSPSetFromOptions_NASH;
751: ksp->ops->buildsolution = KSPBuildSolutionDefault;
752: ksp->ops->buildresidual = KSPBuildResidualDefault;
753: ksp->ops->view = 0;
755: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",KSPNASHSetRadius_NASH);
756: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",KSPNASHGetNormD_NASH);
757: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",KSPNASHGetObjFcn_NASH);
758: return(0);
759: }