Actual source code: nash.c
petsc-3.11.4 2019-09-28
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;
17: Mat Qmat, Mmat;
18: Vec r, z, p, d;
19: PC pc;
21: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
22: PetscReal alpha, beta, kappa, rz, rzm1;
23: PetscReal rr, r2, step;
25: PetscInt max_cg_its;
27: PetscBool diagonalscale;
30: /***************************************************************************/
31: /* Check the arguments and parameters. */
32: /***************************************************************************/
34: PCGetDiagonalScale(ksp->pc, &diagonalscale);
35: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
36: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
38: /***************************************************************************/
39: /* Get the workspace vectors and initialize variables */
40: /***************************************************************************/
42: r2 = cg->radius * cg->radius;
43: r = ksp->work[0];
44: z = ksp->work[1];
45: p = ksp->work[2];
46: d = ksp->vec_sol;
47: pc = ksp->pc;
49: PCGetOperators(pc, &Qmat, &Mmat);
51: VecGetSize(d, &max_cg_its);
52: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
53: ksp->its = 0;
55: /***************************************************************************/
56: /* Initialize objective function and direction. */
57: /***************************************************************************/
59: cg->o_fcn = 0.0;
61: VecSet(d, 0.0); /* d = 0 */
62: cg->norm_d = 0.0;
64: /***************************************************************************/
65: /* Begin the conjugate gradient method. Check the right-hand side for */
66: /* numerical problems. The check for not-a-number and infinite values */
67: /* need be performed only once. */
68: /***************************************************************************/
70: VecCopy(ksp->vec_rhs, r); /* r = -grad */
71: VecDot(r, r, &rr); /* rr = r^T r */
72: KSPCheckDot(ksp,rr);
74: /***************************************************************************/
75: /* Check the preconditioner for numerical problems and for positive */
76: /* definiteness. The check for not-a-number and infinite values need be */
77: /* performed only once. */
78: /***************************************************************************/
80: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
81: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
82: if (PetscIsInfOrNanScalar(rz)) {
83: /*************************************************************************/
84: /* The preconditioner contains not-a-number or an infinite value. */
85: /* Return the gradient direction intersected with the trust region. */
86: /*************************************************************************/
88: ksp->reason = KSP_DIVERGED_NANORINF;
89: PetscInfo1(ksp, "KSPCGSolve_NASH: bad preconditioner: rz=%g\n", (double)rz);
91: if (cg->radius) {
92: if (r2 >= rr) {
93: alpha = 1.0;
94: cg->norm_d = PetscSqrtReal(rr);
95: } else {
96: alpha = PetscSqrtReal(r2 / rr);
97: cg->norm_d = cg->radius;
98: }
100: VecAXPY(d, alpha, r); /* d = d + alpha r */
102: /***********************************************************************/
103: /* Compute objective function. */
104: /***********************************************************************/
106: KSP_MatMult(ksp, Qmat, d, z);
107: VecAYPX(z, -0.5, ksp->vec_rhs);
108: VecDot(d, z, &cg->o_fcn);
109: cg->o_fcn = -cg->o_fcn;
110: ++ksp->its;
111: }
112: return(0);
113: }
115: if (rz < 0.0) {
116: /*************************************************************************/
117: /* The preconditioner is indefinite. Because this is the first */
118: /* and we do not have a direction yet, we use the gradient step. Note */
119: /* that we cannot use the preconditioned norm when computing the step */
120: /* because the matrix is indefinite. */
121: /*************************************************************************/
123: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
124: PetscInfo1(ksp, "KSPCGSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz);
126: if (cg->radius) {
127: if (r2 >= rr) {
128: alpha = 1.0;
129: cg->norm_d = PetscSqrtReal(rr);
130: } else {
131: alpha = PetscSqrtReal(r2 / rr);
132: cg->norm_d = cg->radius;
133: }
135: VecAXPY(d, alpha, r); /* d = d + alpha r */
137: /***********************************************************************/
138: /* Compute objective function. */
139: /***********************************************************************/
141: KSP_MatMult(ksp, Qmat, d, z);
142: VecAYPX(z, -0.5, ksp->vec_rhs);
143: VecDot(d, z, &cg->o_fcn);
144: cg->o_fcn = -cg->o_fcn;
145: ++ksp->its;
146: }
147: return(0);
148: }
150: /***************************************************************************/
151: /* As far as we know, the preconditioner is positive semidefinite. */
152: /* Compute and log the residual. Check convergence because this */
153: /* initializes things, but do not terminate until at least one conjugate */
154: /* gradient iteration has been performed. */
155: /***************************************************************************/
157: switch (ksp->normtype) {
158: case KSP_NORM_PRECONDITIONED:
159: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
160: break;
162: case KSP_NORM_UNPRECONDITIONED:
163: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
164: break;
166: case KSP_NORM_NATURAL:
167: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
168: break;
170: default:
171: norm_r = 0.0;
172: break;
173: }
175: KSPLogResidualHistory(ksp, norm_r);
176: KSPMonitor(ksp, ksp->its, norm_r);
177: ksp->rnorm = norm_r;
179: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
181: /***************************************************************************/
182: /* Compute the first direction and update the iteration. */
183: /***************************************************************************/
185: VecCopy(z, p); /* p = z */
186: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
187: ++ksp->its;
189: /***************************************************************************/
190: /* Check the matrix for numerical problems. */
191: /***************************************************************************/
193: VecDot(p, z, &kappa); /* kappa = p^T Q p */
194: if (PetscIsInfOrNanScalar(kappa)) {
195: /*************************************************************************/
196: /* The matrix produced not-a-number or an infinite value. In this case, */
197: /* we must stop and use the gradient direction. This condition need */
198: /* only be checked once. */
199: /*************************************************************************/
201: ksp->reason = KSP_DIVERGED_NANORINF;
202: PetscInfo1(ksp, "KSPCGSolve_NASH: bad matrix: kappa=%g\n", (double)kappa);
204: if (cg->radius) {
205: if (r2 >= rr) {
206: alpha = 1.0;
207: cg->norm_d = PetscSqrtReal(rr);
208: } else {
209: alpha = PetscSqrtReal(r2 / rr);
210: cg->norm_d = cg->radius;
211: }
213: VecAXPY(d, alpha, r); /* d = d + alpha r */
215: /***********************************************************************/
216: /* Compute objective function. */
217: /***********************************************************************/
219: KSP_MatMult(ksp, Qmat, d, z);
220: VecAYPX(z, -0.5, ksp->vec_rhs);
221: VecDot(d, z, &cg->o_fcn);
222: cg->o_fcn = -cg->o_fcn;
223: ++ksp->its;
224: }
225: return(0);
226: }
228: /***************************************************************************/
229: /* Initialize variables for calculating the norm of the direction. */
230: /***************************************************************************/
232: dMp = 0.0;
233: norm_d = 0.0;
234: switch (cg->dtype) {
235: case NASH_PRECONDITIONED_DIRECTION:
236: norm_p = rz;
237: break;
239: default:
240: VecDot(p, p, &norm_p);
241: break;
242: }
244: /***************************************************************************/
245: /* Check for negative curvature. */
246: /***************************************************************************/
248: if (kappa <= 0.0) {
249: /*************************************************************************/
250: /* In this case, the matrix is indefinite and we have encountered a */
251: /* direction of negative curvature. Because negative curvature occurs */
252: /* during the first step, we must follow a direction. */
253: /*************************************************************************/
255: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
256: PetscInfo1(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
258: if (cg->radius && norm_p > 0.0) {
259: /***********************************************************************/
260: /* Follow direction of negative curvature to the boundary of the */
261: /* trust region. */
262: /***********************************************************************/
264: step = PetscSqrtReal(r2 / norm_p);
265: cg->norm_d = cg->radius;
267: VecAXPY(d, step, p); /* d = d + step p */
269: /***********************************************************************/
270: /* Update objective function. */
271: /***********************************************************************/
273: cg->o_fcn += step * (0.5 * step * kappa - rz);
274: } else if (cg->radius) {
275: /***********************************************************************/
276: /* The norm of the preconditioned direction is zero; use the gradient */
277: /* step. */
278: /***********************************************************************/
280: if (r2 >= rr) {
281: alpha = 1.0;
282: cg->norm_d = PetscSqrtReal(rr);
283: } else {
284: alpha = PetscSqrtReal(r2 / rr);
285: cg->norm_d = cg->radius;
286: }
288: VecAXPY(d, alpha, r); /* d = d + alpha r */
290: /***********************************************************************/
291: /* Compute objective function. */
292: /***********************************************************************/
294: KSP_MatMult(ksp, Qmat, d, z);
295: VecAYPX(z, -0.5, ksp->vec_rhs);
296: VecDot(d, z, &cg->o_fcn);
297: cg->o_fcn = -cg->o_fcn;
298: ++ksp->its;
299: }
300: return(0);
301: }
303: /***************************************************************************/
304: /* Run the conjugate gradient method until either the problem is solved, */
305: /* we encounter the boundary of the trust region, or the conjugate */
306: /* gradient method breaks down. */
307: /***************************************************************************/
309: while (1) {
310: /*************************************************************************/
311: /* Know that kappa is nonzero, because we have not broken down, so we */
312: /* can compute the steplength. */
313: /*************************************************************************/
315: alpha = rz / kappa;
317: /*************************************************************************/
318: /* Compute the steplength and check for intersection with the trust */
319: /* region. */
320: /*************************************************************************/
322: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
323: if (cg->radius && norm_dp1 >= r2) {
324: /***********************************************************************/
325: /* In this case, the matrix is positive definite as far as we know. */
326: /* However, the full step goes beyond the trust region. */
327: /***********************************************************************/
329: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
330: PetscInfo1(ksp, "KSPCGSolve_NASH: constrained step: radius=%g\n", (double)cg->radius);
332: if (norm_p > 0.0) {
333: /*********************************************************************/
334: /* Follow the direction to the boundary of the trust region. */
335: /*********************************************************************/
337: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
338: cg->norm_d = cg->radius;
340: VecAXPY(d, step, p); /* d = d + step p */
342: /*********************************************************************/
343: /* Update objective function. */
344: /*********************************************************************/
346: cg->o_fcn += step * (0.5 * step * kappa - rz);
347: } else {
348: /*********************************************************************/
349: /* The norm of the direction is zero; there is nothing to follow. */
350: /*********************************************************************/
351: }
352: break;
353: }
355: /*************************************************************************/
356: /* Now we can update the direction and residual. */
357: /*************************************************************************/
359: VecAXPY(d, alpha, p); /* d = d + alpha p */
360: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
361: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
363: switch (cg->dtype) {
364: case NASH_PRECONDITIONED_DIRECTION:
365: norm_d = norm_dp1;
366: break;
368: default:
369: VecDot(d, d, &norm_d);
370: break;
371: }
372: cg->norm_d = PetscSqrtReal(norm_d);
374: /*************************************************************************/
375: /* Update objective function. */
376: /*************************************************************************/
378: cg->o_fcn -= 0.5 * alpha * rz;
380: /*************************************************************************/
381: /* Check that the preconditioner appears positive semidefinite. */
382: /*************************************************************************/
384: rzm1 = rz;
385: VecDot(r, z, &rz); /* rz = r^T z */
386: if (rz < 0.0) {
387: /***********************************************************************/
388: /* The preconditioner is indefinite. */
389: /***********************************************************************/
391: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
392: PetscInfo1(ksp, "KSPCGSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz);
393: break;
394: }
396: /*************************************************************************/
397: /* As far as we know, the preconditioner is positive semidefinite. */
398: /* Compute the residual and check for convergence. */
399: /*************************************************************************/
401: switch (ksp->normtype) {
402: case KSP_NORM_PRECONDITIONED:
403: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
404: break;
406: case KSP_NORM_UNPRECONDITIONED:
407: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
408: break;
410: case KSP_NORM_NATURAL:
411: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
412: break;
414: default:
415: norm_r = 0.;
416: break;
417: }
419: KSPLogResidualHistory(ksp, norm_r);
420: KSPMonitor(ksp, ksp->its, norm_r);
421: ksp->rnorm = norm_r;
423: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
424: if (ksp->reason) {
425: /***********************************************************************/
426: /* The method has converged. */
427: /***********************************************************************/
429: PetscInfo2(ksp, "KSPCGSolve_NASH: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
430: break;
431: }
433: /*************************************************************************/
434: /* We have not converged yet. Check for breakdown. */
435: /*************************************************************************/
437: beta = rz / rzm1;
438: if (PetscAbsReal(beta) <= 0.0) {
439: /***********************************************************************/
440: /* Conjugate gradients has broken down. */
441: /***********************************************************************/
443: ksp->reason = KSP_DIVERGED_BREAKDOWN;
444: PetscInfo1(ksp, "KSPCGSolve_NASH: breakdown: beta=%g\n", (double)beta);
445: break;
446: }
448: /*************************************************************************/
449: /* Check iteration limit. */
450: /*************************************************************************/
452: if (ksp->its >= max_cg_its) {
453: ksp->reason = KSP_DIVERGED_ITS;
454: PetscInfo1(ksp, "KSPCGSolve_NASH: iterlim: its=%D\n", ksp->its);
455: break;
456: }
458: /*************************************************************************/
459: /* Update p and the norms. */
460: /*************************************************************************/
462: VecAYPX(p, beta, z); /* p = z + beta p */
464: switch (cg->dtype) {
465: case NASH_PRECONDITIONED_DIRECTION:
466: dMp = beta*(dMp + alpha*norm_p);
467: norm_p = beta*(rzm1 + beta*norm_p);
468: break;
470: default:
471: VecDot(d, p, &dMp);
472: VecDot(p, p, &norm_p);
473: break;
474: }
476: /*************************************************************************/
477: /* Compute the new direction and update the iteration. */
478: /*************************************************************************/
480: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
481: VecDot(p, z, &kappa); /* kappa = p^T Q p */
482: ++ksp->its;
484: /*************************************************************************/
485: /* Check for negative curvature. */
486: /*************************************************************************/
488: if (kappa <= 0.0) {
489: /***********************************************************************/
490: /* In this case, the matrix is indefinite and we have encountered */
491: /* a direction of negative curvature. Stop at the base. */
492: /***********************************************************************/
494: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
495: PetscInfo1(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
496: break;
497: }
498: }
499: return(0);
500: #endif
501: }
503: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
504: {
507: /***************************************************************************/
508: /* Set work vectors needed by conjugate gradient method and allocate */
509: /***************************************************************************/
512: KSPSetWorkVecs(ksp,3);
513: return(0);
514: }
516: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
517: {
521: /***************************************************************************/
522: /* Clear composed functions */
523: /***************************************************************************/
525: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
526: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
527: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
529: /***************************************************************************/
530: /* Destroy KSP object. */
531: /***************************************************************************/
533: KSPDestroyDefault(ksp);
534: return(0);
535: }
537: static PetscErrorCode KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
538: {
539: KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;
542: cg->radius = radius;
543: return(0);
544: }
546: static PetscErrorCode KSPCGGetNormD_NASH(KSP ksp, PetscReal *norm_d)
547: {
548: KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;
551: *norm_d = cg->norm_d;
552: return(0);
553: }
555: static PetscErrorCode KSPCGGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
556: {
557: KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;
560: *o_fcn = cg->o_fcn;
561: return(0);
562: }
564: static PetscErrorCode KSPCGSetFromOptions_NASH(PetscOptionItems *PetscOptionsObject,KSP ksp)
565: {
567: KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;
570: PetscOptionsHead(PetscOptionsObject,"KSPCG NASH options");
572: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
574: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
576: PetscOptionsTail();
577: return(0);
578: }
580: /*MC
581: KSPCGNASH - Code to run conjugate gradient method subject to a constraint
582: on the solution norm. This is used in Trust Region methods for
583: nonlinear equations, SNESNEWTONTR
585: Options Database Keys:
586: . -ksp_cg_radius <r> - Trust Region Radius
588: Notes:
589: This is rarely used directly
591: Level: developer
593: Use preconditioned conjugate gradient to compute
594: an approximate minimizer of the quadratic function
596: q(s) = g^T * s + 0.5 * s^T * H * s
598: subject to the trust region constraint
600: || s || <= delta,
602: where
604: delta is the trust region radius,
605: g is the gradient vector,
606: H is the Hessian approximation, and
607: M is the positive definite preconditioner matrix.
609: KSPConvergedReason may be
610: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
611: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
612: $ other KSP converged/diverged reasons
614: Notes:
615: The preconditioner supplied should be symmetric and positive definite.
617: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn()
618: M*/
620: PETSC_EXTERN PetscErrorCode KSPCreate_CGNASH(KSP ksp)
621: {
623: KSPCG_NASH *cg;
626: PetscNewLog(ksp,&cg);
627: cg->radius = 0.0;
628: cg->dtype = NASH_UNPRECONDITIONED_DIRECTION;
630: ksp->data = (void*) cg;
631: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
632: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
633: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
635: /***************************************************************************/
636: /* Sets the functions that are associated with this data structure */
637: /* (in C++ this is the same as defining virtual functions). */
638: /***************************************************************************/
640: ksp->ops->setup = KSPCGSetUp_NASH;
641: ksp->ops->solve = KSPCGSolve_NASH;
642: ksp->ops->destroy = KSPCGDestroy_NASH;
643: ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
644: ksp->ops->buildsolution = KSPBuildSolutionDefault;
645: ksp->ops->buildresidual = KSPBuildResidualDefault;
646: ksp->ops->view = 0;
648: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_NASH);
649: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",KSPCGGetNormD_NASH);
650: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_NASH);
651: return(0);
652: }