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: /***************************************************************************/
29: /* Check the arguments and parameters. */
30: /***************************************************************************/
32: PCGetDiagonalScale(ksp->pc, &diagonalscale);
36: /***************************************************************************/
37: /* Get the workspace vectors and initialize variables */
38: /***************************************************************************/
40: r2 = cg->radius * cg->radius;
41: r = ksp->work[0];
42: z = ksp->work[1];
43: p = ksp->work[2];
44: d = ksp->vec_sol;
45: pc = ksp->pc;
47: PCGetOperators(pc, &Qmat, &Mmat);
49: VecGetSize(d, &max_cg_its);
50: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
51: ksp->its = 0;
53: /***************************************************************************/
54: /* Initialize objective function and direction. */
55: /***************************************************************************/
57: cg->o_fcn = 0.0;
59: VecSet(d, 0.0); /* d = 0 */
60: cg->norm_d = 0.0;
62: /***************************************************************************/
63: /* Begin the conjugate gradient method. Check the right-hand side for */
64: /* numerical problems. The check for not-a-number and infinite values */
65: /* need be performed only once. */
66: /***************************************************************************/
68: VecCopy(ksp->vec_rhs, r); /* r = -grad */
69: VecDot(r, r, &rr); /* rr = r^T r */
70: KSPCheckDot(ksp, rr);
72: /***************************************************************************/
73: /* Check the preconditioner for numerical problems and for positive */
74: /* definiteness. The check for not-a-number and infinite values need be */
75: /* performed only once. */
76: /***************************************************************************/
78: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
79: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
80: if (PetscIsInfOrNanScalar(rz)) {
81: /*************************************************************************/
82: /* The preconditioner contains not-a-number or an infinite value. */
83: /* Return the gradient direction intersected with the trust region. */
84: /*************************************************************************/
86: ksp->reason = KSP_DIVERGED_NANORINF;
87: PetscInfo(ksp, "KSPCGSolve_NASH: bad preconditioner: rz=%g\n", (double)rz);
89: if (cg->radius) {
90: if (r2 >= rr) {
91: alpha = 1.0;
92: cg->norm_d = PetscSqrtReal(rr);
93: } else {
94: alpha = PetscSqrtReal(r2 / rr);
95: cg->norm_d = cg->radius;
96: }
98: VecAXPY(d, alpha, r); /* d = d + alpha r */
100: /***********************************************************************/
101: /* Compute objective function. */
102: /***********************************************************************/
104: KSP_MatMult(ksp, Qmat, d, z);
105: VecAYPX(z, -0.5, ksp->vec_rhs);
106: VecDot(d, z, &cg->o_fcn);
107: cg->o_fcn = -cg->o_fcn;
108: ++ksp->its;
109: }
110: return 0;
111: }
113: if (rz < 0.0) {
114: /*************************************************************************/
115: /* The preconditioner is indefinite. Because this is the first */
116: /* and we do not have a direction yet, we use the gradient step. Note */
117: /* that we cannot use the preconditioned norm when computing the step */
118: /* because the matrix is indefinite. */
119: /*************************************************************************/
121: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
122: PetscInfo(ksp, "KSPCGSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz);
124: if (cg->radius) {
125: if (r2 >= rr) {
126: alpha = 1.0;
127: cg->norm_d = PetscSqrtReal(rr);
128: } else {
129: alpha = PetscSqrtReal(r2 / rr);
130: cg->norm_d = cg->radius;
131: }
133: VecAXPY(d, alpha, r); /* d = d + alpha r */
135: /***********************************************************************/
136: /* Compute objective function. */
137: /***********************************************************************/
139: KSP_MatMult(ksp, Qmat, d, z);
140: VecAYPX(z, -0.5, ksp->vec_rhs);
141: VecDot(d, z, &cg->o_fcn);
142: cg->o_fcn = -cg->o_fcn;
143: ++ksp->its;
144: }
145: return 0;
146: }
148: /***************************************************************************/
149: /* As far as we know, the preconditioner is positive semidefinite. */
150: /* Compute and log the residual. Check convergence because this */
151: /* initializes things, but do not terminate until at least one conjugate */
152: /* gradient iteration has been performed. */
153: /***************************************************************************/
155: switch (ksp->normtype) {
156: case KSP_NORM_PRECONDITIONED:
157: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
158: break;
160: case KSP_NORM_UNPRECONDITIONED:
161: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
162: break;
164: case KSP_NORM_NATURAL:
165: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
166: break;
168: default:
169: norm_r = 0.0;
170: break;
171: }
173: KSPLogResidualHistory(ksp, norm_r);
174: KSPMonitor(ksp, ksp->its, norm_r);
175: ksp->rnorm = norm_r;
177: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
179: /***************************************************************************/
180: /* Compute the first direction and update the iteration. */
181: /***************************************************************************/
183: VecCopy(z, p); /* p = z */
184: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
185: ++ksp->its;
187: /***************************************************************************/
188: /* Check the matrix for numerical problems. */
189: /***************************************************************************/
191: VecDot(p, z, &kappa); /* kappa = p^T Q p */
192: if (PetscIsInfOrNanScalar(kappa)) {
193: /*************************************************************************/
194: /* The matrix produced not-a-number or an infinite value. In this case, */
195: /* we must stop and use the gradient direction. This condition need */
196: /* only be checked once. */
197: /*************************************************************************/
199: ksp->reason = KSP_DIVERGED_NANORINF;
200: PetscInfo(ksp, "KSPCGSolve_NASH: bad matrix: kappa=%g\n", (double)kappa);
202: if (cg->radius) {
203: if (r2 >= rr) {
204: alpha = 1.0;
205: cg->norm_d = PetscSqrtReal(rr);
206: } else {
207: alpha = PetscSqrtReal(r2 / rr);
208: cg->norm_d = cg->radius;
209: }
211: VecAXPY(d, alpha, r); /* d = d + alpha r */
213: /***********************************************************************/
214: /* Compute objective function. */
215: /***********************************************************************/
217: KSP_MatMult(ksp, Qmat, d, z);
218: VecAYPX(z, -0.5, ksp->vec_rhs);
219: VecDot(d, z, &cg->o_fcn);
220: cg->o_fcn = -cg->o_fcn;
221: ++ksp->its;
222: }
223: return 0;
224: }
226: /***************************************************************************/
227: /* Initialize variables for calculating the norm of the direction. */
228: /***************************************************************************/
230: dMp = 0.0;
231: norm_d = 0.0;
232: switch (cg->dtype) {
233: case NASH_PRECONDITIONED_DIRECTION:
234: norm_p = rz;
235: break;
237: default:
238: VecDot(p, p, &norm_p);
239: break;
240: }
242: /***************************************************************************/
243: /* Check for negative curvature. */
244: /***************************************************************************/
246: if (kappa <= 0.0) {
247: /*************************************************************************/
248: /* In this case, the matrix is indefinite and we have encountered a */
249: /* direction of negative curvature. Because negative curvature occurs */
250: /* during the first step, we must follow a direction. */
251: /*************************************************************************/
253: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
254: PetscInfo(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
256: if (cg->radius && norm_p > 0.0) {
257: /***********************************************************************/
258: /* Follow direction of negative curvature to the boundary of the */
259: /* trust region. */
260: /***********************************************************************/
262: step = PetscSqrtReal(r2 / norm_p);
263: cg->norm_d = cg->radius;
265: VecAXPY(d, step, p); /* d = d + step p */
267: /***********************************************************************/
268: /* Update objective function. */
269: /***********************************************************************/
271: cg->o_fcn += step * (0.5 * step * kappa - rz);
272: } else if (cg->radius) {
273: /***********************************************************************/
274: /* The norm of the preconditioned direction is zero; use the gradient */
275: /* step. */
276: /***********************************************************************/
278: if (r2 >= rr) {
279: alpha = 1.0;
280: cg->norm_d = PetscSqrtReal(rr);
281: } else {
282: alpha = PetscSqrtReal(r2 / rr);
283: cg->norm_d = cg->radius;
284: }
286: VecAXPY(d, alpha, r); /* d = d + alpha r */
288: /***********************************************************************/
289: /* Compute objective function. */
290: /***********************************************************************/
292: KSP_MatMult(ksp, Qmat, d, z);
293: VecAYPX(z, -0.5, ksp->vec_rhs);
294: VecDot(d, z, &cg->o_fcn);
295: cg->o_fcn = -cg->o_fcn;
296: ++ksp->its;
297: }
298: return 0;
299: }
301: /***************************************************************************/
302: /* Run the conjugate gradient method until either the problem is solved, */
303: /* we encounter the boundary of the trust region, or the conjugate */
304: /* gradient method breaks down. */
305: /***************************************************************************/
307: while (1) {
308: /*************************************************************************/
309: /* Know that kappa is nonzero, because we have not broken down, so we */
310: /* can compute the steplength. */
311: /*************************************************************************/
313: alpha = rz / kappa;
315: /*************************************************************************/
316: /* Compute the steplength and check for intersection with the trust */
317: /* region. */
318: /*************************************************************************/
320: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
321: if (cg->radius && norm_dp1 >= r2) {
322: /***********************************************************************/
323: /* In this case, the matrix is positive definite as far as we know. */
324: /* However, the full step goes beyond the trust region. */
325: /***********************************************************************/
327: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
328: PetscInfo(ksp, "KSPCGSolve_NASH: constrained step: radius=%g\n", (double)cg->radius);
330: if (norm_p > 0.0) {
331: /*********************************************************************/
332: /* Follow the direction to the boundary of the trust region. */
333: /*********************************************************************/
335: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
336: cg->norm_d = cg->radius;
338: VecAXPY(d, step, p); /* d = d + step p */
340: /*********************************************************************/
341: /* Update objective function. */
342: /*********************************************************************/
344: cg->o_fcn += step * (0.5 * step * kappa - rz);
345: } else {
346: /*********************************************************************/
347: /* The norm of the direction is zero; there is nothing to follow. */
348: /*********************************************************************/
349: }
350: break;
351: }
353: /*************************************************************************/
354: /* Now we can update the direction and residual. */
355: /*************************************************************************/
357: VecAXPY(d, alpha, p); /* d = d + alpha p */
358: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
359: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
361: switch (cg->dtype) {
362: case NASH_PRECONDITIONED_DIRECTION:
363: norm_d = norm_dp1;
364: break;
366: default:
367: VecDot(d, d, &norm_d);
368: break;
369: }
370: cg->norm_d = PetscSqrtReal(norm_d);
372: /*************************************************************************/
373: /* Update objective function. */
374: /*************************************************************************/
376: cg->o_fcn -= 0.5 * alpha * rz;
378: /*************************************************************************/
379: /* Check that the preconditioner appears positive semidefinite. */
380: /*************************************************************************/
382: rzm1 = rz;
383: VecDot(r, z, &rz); /* rz = r^T z */
384: if (rz < 0.0) {
385: /***********************************************************************/
386: /* The preconditioner is indefinite. */
387: /***********************************************************************/
389: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
390: PetscInfo(ksp, "KSPCGSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz);
391: break;
392: }
394: /*************************************************************************/
395: /* As far as we know, the preconditioner is positive semidefinite. */
396: /* Compute the residual and check for convergence. */
397: /*************************************************************************/
399: switch (ksp->normtype) {
400: case KSP_NORM_PRECONDITIONED:
401: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
402: break;
404: case KSP_NORM_UNPRECONDITIONED:
405: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
406: break;
408: case KSP_NORM_NATURAL:
409: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
410: break;
412: default:
413: norm_r = 0.;
414: break;
415: }
417: KSPLogResidualHistory(ksp, norm_r);
418: KSPMonitor(ksp, ksp->its, norm_r);
419: ksp->rnorm = norm_r;
421: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
422: if (ksp->reason) {
423: /***********************************************************************/
424: /* The method has converged. */
425: /***********************************************************************/
427: PetscInfo(ksp, "KSPCGSolve_NASH: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
428: break;
429: }
431: /*************************************************************************/
432: /* We have not converged yet. Check for breakdown. */
433: /*************************************************************************/
435: beta = rz / rzm1;
436: if (PetscAbsReal(beta) <= 0.0) {
437: /***********************************************************************/
438: /* Conjugate gradients has broken down. */
439: /***********************************************************************/
441: ksp->reason = KSP_DIVERGED_BREAKDOWN;
442: PetscInfo(ksp, "KSPCGSolve_NASH: breakdown: beta=%g\n", (double)beta);
443: break;
444: }
446: /*************************************************************************/
447: /* Check iteration limit. */
448: /*************************************************************************/
450: if (ksp->its >= max_cg_its) {
451: ksp->reason = KSP_DIVERGED_ITS;
452: PetscInfo(ksp, "KSPCGSolve_NASH: iterlim: its=%" PetscInt_FMT "\n", ksp->its);
453: break;
454: }
456: /*************************************************************************/
457: /* Update p and the norms. */
458: /*************************************************************************/
460: VecAYPX(p, beta, z); /* p = z + beta p */
462: switch (cg->dtype) {
463: case NASH_PRECONDITIONED_DIRECTION:
464: dMp = beta * (dMp + alpha * norm_p);
465: norm_p = beta * (rzm1 + beta * norm_p);
466: break;
468: default:
469: VecDot(d, p, &dMp);
470: VecDot(p, p, &norm_p);
471: break;
472: }
474: /*************************************************************************/
475: /* Compute the new direction and update the iteration. */
476: /*************************************************************************/
478: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
479: VecDot(p, z, &kappa); /* kappa = p^T Q p */
480: ++ksp->its;
482: /*************************************************************************/
483: /* Check for negative curvature. */
484: /*************************************************************************/
486: if (kappa <= 0.0) {
487: /***********************************************************************/
488: /* In this case, the matrix is indefinite and we have encountered */
489: /* a direction of negative curvature. Stop at the base. */
490: /***********************************************************************/
492: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
493: PetscInfo(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
494: break;
495: }
496: }
497: return 0;
498: #endif
499: }
501: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
502: {
503: /***************************************************************************/
504: /* Set work vectors needed by conjugate gradient method and allocate */
505: /***************************************************************************/
507: KSPSetWorkVecs(ksp, 3);
508: return 0;
509: }
511: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
512: {
513: /***************************************************************************/
514: /* Clear composed functions */
515: /***************************************************************************/
517: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL);
518: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL);
519: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL);
521: /***************************************************************************/
522: /* Destroy KSP object. */
523: /***************************************************************************/
525: KSPDestroyDefault(ksp);
526: return 0;
527: }
529: static PetscErrorCode KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
530: {
531: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
533: cg->radius = radius;
534: return 0;
535: }
537: static PetscErrorCode KSPCGGetNormD_NASH(KSP ksp, PetscReal *norm_d)
538: {
539: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
541: *norm_d = cg->norm_d;
542: return 0;
543: }
545: static PetscErrorCode KSPCGGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
546: {
547: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
549: *o_fcn = cg->o_fcn;
550: return 0;
551: }
553: static PetscErrorCode KSPCGSetFromOptions_NASH(KSP ksp, PetscOptionItems *PetscOptionsObject)
554: {
555: KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;
557: PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG NASH options");
559: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
561: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
563: PetscOptionsHeadEnd();
564: return 0;
565: }
567: /*MC
568: KSPNASH - Code to run conjugate gradient method subject to a constraint on the solution norm.
570: Options Database Keys:
571: . -ksp_cg_radius <r> - Trust Region Radius
573: Level: developer
575: Notes:
576: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
578: Uses preconditioned conjugate gradient to compute
579: an approximate minimizer of the quadratic function
581: q(s) = g^T * s + 0.5 * s^T * H * s
583: subject to the trust region constraint
585: || s || <= delta,
587: where
589: delta is the trust region radius,
590: g is the gradient vector,
591: H is the Hessian approximation, and
592: M is the positive definite preconditioner matrix.
594: `KSPConvergedReason` may be
595: .vb
596: KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
597: KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
598: .ve
599: other `KSP` converged/diverged reasons
601: The preconditioner supplied should be symmetric and positive definite.
603: Reference:
604: Nash, Stephen G. Newton-type minimization via the Lanczos method. SIAM Journal on Numerical Analysis 21, no. 4 (1984): 770-788.
606: .seealso: [](chapter_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: PetscNew(&cg);
614: cg->radius = 0.0;
615: cg->dtype = NASH_UNPRECONDITIONED_DIRECTION;
617: ksp->data = (void *)cg;
618: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
619: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
620: KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
621: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
623: /***************************************************************************/
624: /* Sets the functions that are associated with this data structure */
625: /* (in C++ this is the same as defining virtual functions). */
626: /***************************************************************************/
628: ksp->ops->setup = KSPCGSetUp_NASH;
629: ksp->ops->solve = KSPCGSolve_NASH;
630: ksp->ops->destroy = KSPCGDestroy_NASH;
631: ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
632: ksp->ops->buildsolution = KSPBuildSolutionDefault;
633: ksp->ops->buildresidual = KSPBuildResidualDefault;
634: ksp->ops->view = NULL;
636: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_NASH);
637: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_NASH);
638: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_NASH);
639: return 0;
640: }