Actual source code: nash.c
petsc-3.7.5 2017-01-01
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: Mat Qmat, Mmat;
102: Vec r, z, p, d;
103: PC pc;
105: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
106: PetscReal alpha, beta, kappa, rz, rzm1;
107: PetscReal rr, r2, step;
109: PetscInt max_cg_its;
111: PetscBool diagonalscale;
114: /***************************************************************************/
115: /* Check the arguments and parameters. */
116: /***************************************************************************/
118: PCGetDiagonalScale(ksp->pc, &diagonalscale);
119: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
120: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
122: /***************************************************************************/
123: /* Get the workspace vectors and initialize variables */
124: /***************************************************************************/
126: r2 = cg->radius * cg->radius;
127: r = ksp->work[0];
128: z = ksp->work[1];
129: p = ksp->work[2];
130: d = ksp->vec_sol;
131: pc = ksp->pc;
133: PCGetOperators(pc, &Qmat, &Mmat);
135: VecGetSize(d, &max_cg_its);
136: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
137: ksp->its = 0;
139: /***************************************************************************/
140: /* Initialize objective function and direction. */
141: /***************************************************************************/
143: cg->o_fcn = 0.0;
145: VecSet(d, 0.0); /* d = 0 */
146: cg->norm_d = 0.0;
148: /***************************************************************************/
149: /* Begin the conjugate gradient method. Check the right-hand side for */
150: /* numerical problems. The check for not-a-number and infinite values */
151: /* need be performed only once. */
152: /***************************************************************************/
154: VecCopy(ksp->vec_rhs, r); /* r = -grad */
155: VecDot(r, r, &rr); /* rr = r^T r */
156: KSPCheckDot(ksp,rr);
158: /***************************************************************************/
159: /* Check the preconditioner for numerical problems and for positive */
160: /* definiteness. The check for not-a-number and infinite values need be */
161: /* performed only once. */
162: /***************************************************************************/
164: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
165: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
166: if (PetscIsInfOrNanScalar(rz)) {
167: /*************************************************************************/
168: /* The preconditioner contains not-a-number or an infinite value. */
169: /* Return the gradient direction intersected with the trust region. */
170: /*************************************************************************/
172: ksp->reason = KSP_DIVERGED_NANORINF;
173: PetscInfo1(ksp, "KSPSolve_NASH: bad preconditioner: rz=%g\n", (double)rz);
175: if (cg->radius) {
176: if (r2 >= rr) {
177: alpha = 1.0;
178: cg->norm_d = PetscSqrtReal(rr);
179: } else {
180: alpha = PetscSqrtReal(r2 / rr);
181: cg->norm_d = cg->radius;
182: }
184: VecAXPY(d, alpha, r); /* d = d + alpha r */
186: /***********************************************************************/
187: /* Compute objective function. */
188: /***********************************************************************/
190: KSP_MatMult(ksp, Qmat, d, z);
191: VecAYPX(z, -0.5, ksp->vec_rhs);
192: VecDot(d, z, &cg->o_fcn);
193: cg->o_fcn = -cg->o_fcn;
194: ++ksp->its;
195: }
196: return(0);
197: }
199: if (rz < 0.0) {
200: /*************************************************************************/
201: /* The preconditioner is indefinite. Because this is the first */
202: /* and we do not have a direction yet, we use the gradient step. Note */
203: /* that we cannot use the preconditioned norm when computing the step */
204: /* because the matrix is indefinite. */
205: /*************************************************************************/
207: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
208: PetscInfo1(ksp, "KSPSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz);
210: if (cg->radius) {
211: if (r2 >= rr) {
212: alpha = 1.0;
213: cg->norm_d = PetscSqrtReal(rr);
214: } else {
215: alpha = PetscSqrtReal(r2 / rr);
216: cg->norm_d = cg->radius;
217: }
219: VecAXPY(d, alpha, r); /* d = d + alpha r */
221: /***********************************************************************/
222: /* Compute objective function. */
223: /***********************************************************************/
225: KSP_MatMult(ksp, Qmat, d, z);
226: VecAYPX(z, -0.5, ksp->vec_rhs);
227: VecDot(d, z, &cg->o_fcn);
228: cg->o_fcn = -cg->o_fcn;
229: ++ksp->its;
230: }
231: return(0);
232: }
234: /***************************************************************************/
235: /* As far as we know, the preconditioner is positive semidefinite. */
236: /* Compute and log the residual. Check convergence because this */
237: /* initializes things, but do not terminate until at least one conjugate */
238: /* gradient iteration has been performed. */
239: /***************************************************************************/
241: switch (ksp->normtype) {
242: case KSP_NORM_PRECONDITIONED:
243: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
244: break;
246: case KSP_NORM_UNPRECONDITIONED:
247: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
248: break;
250: case KSP_NORM_NATURAL:
251: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
252: break;
254: default:
255: norm_r = 0.0;
256: break;
257: }
259: KSPLogResidualHistory(ksp, norm_r);
260: KSPMonitor(ksp, ksp->its, norm_r);
261: ksp->rnorm = norm_r;
263: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
265: /***************************************************************************/
266: /* Compute the first direction and update the iteration. */
267: /***************************************************************************/
269: VecCopy(z, p); /* p = z */
270: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
271: ++ksp->its;
273: /***************************************************************************/
274: /* Check the matrix for numerical problems. */
275: /***************************************************************************/
277: VecDot(p, z, &kappa); /* kappa = p^T Q p */
278: if (PetscIsInfOrNanScalar(kappa)) {
279: /*************************************************************************/
280: /* The matrix produced not-a-number or an infinite value. In this case, */
281: /* we must stop and use the gradient direction. This condition need */
282: /* only be checked once. */
283: /*************************************************************************/
285: ksp->reason = KSP_DIVERGED_NANORINF;
286: PetscInfo1(ksp, "KSPSolve_NASH: bad matrix: kappa=%g\n", (double)kappa);
288: if (cg->radius) {
289: if (r2 >= rr) {
290: alpha = 1.0;
291: cg->norm_d = PetscSqrtReal(rr);
292: } else {
293: alpha = PetscSqrtReal(r2 / rr);
294: cg->norm_d = cg->radius;
295: }
297: VecAXPY(d, alpha, r); /* d = d + alpha r */
299: /***********************************************************************/
300: /* Compute objective function. */
301: /***********************************************************************/
303: KSP_MatMult(ksp, Qmat, d, z);
304: VecAYPX(z, -0.5, ksp->vec_rhs);
305: VecDot(d, z, &cg->o_fcn);
306: cg->o_fcn = -cg->o_fcn;
307: ++ksp->its;
308: }
309: return(0);
310: }
312: /***************************************************************************/
313: /* Initialize variables for calculating the norm of the direction. */
314: /***************************************************************************/
316: dMp = 0.0;
317: norm_d = 0.0;
318: switch (cg->dtype) {
319: case NASH_PRECONDITIONED_DIRECTION:
320: norm_p = rz;
321: break;
323: default:
324: VecDot(p, p, &norm_p);
325: break;
326: }
328: /***************************************************************************/
329: /* Check for negative curvature. */
330: /***************************************************************************/
332: if (kappa <= 0.0) {
333: /*************************************************************************/
334: /* In this case, the matrix is indefinite and we have encountered a */
335: /* direction of negative curvature. Because negative curvature occurs */
336: /* during the first step, we must follow a direction. */
337: /*************************************************************************/
339: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
340: PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
342: if (cg->radius && norm_p > 0.0) {
343: /***********************************************************************/
344: /* Follow direction of negative curvature to the boundary of the */
345: /* trust region. */
346: /***********************************************************************/
348: step = PetscSqrtReal(r2 / norm_p);
349: cg->norm_d = cg->radius;
351: VecAXPY(d, step, p); /* d = d + step p */
353: /***********************************************************************/
354: /* Update objective function. */
355: /***********************************************************************/
357: cg->o_fcn += step * (0.5 * step * kappa - rz);
358: } else if (cg->radius) {
359: /***********************************************************************/
360: /* The norm of the preconditioned direction is zero; use the gradient */
361: /* step. */
362: /***********************************************************************/
364: if (r2 >= rr) {
365: alpha = 1.0;
366: cg->norm_d = PetscSqrtReal(rr);
367: } else {
368: alpha = PetscSqrtReal(r2 / rr);
369: cg->norm_d = cg->radius;
370: }
372: VecAXPY(d, alpha, r); /* d = d + alpha r */
374: /***********************************************************************/
375: /* Compute objective function. */
376: /***********************************************************************/
378: KSP_MatMult(ksp, Qmat, d, z);
379: VecAYPX(z, -0.5, ksp->vec_rhs);
380: VecDot(d, z, &cg->o_fcn);
381: cg->o_fcn = -cg->o_fcn;
382: ++ksp->its;
383: }
384: return(0);
385: }
387: /***************************************************************************/
388: /* Run the conjugate gradient method until either the problem is solved, */
389: /* we encounter the boundary of the trust region, or the conjugate */
390: /* gradient method breaks down. */
391: /***************************************************************************/
393: while (1) {
394: /*************************************************************************/
395: /* Know that kappa is nonzero, because we have not broken down, so we */
396: /* can compute the steplength. */
397: /*************************************************************************/
399: alpha = rz / kappa;
401: /*************************************************************************/
402: /* Compute the steplength and check for intersection with the trust */
403: /* region. */
404: /*************************************************************************/
406: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
407: if (cg->radius && norm_dp1 >= r2) {
408: /***********************************************************************/
409: /* In this case, the matrix is positive definite as far as we know. */
410: /* However, the full step goes beyond the trust region. */
411: /***********************************************************************/
413: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
414: PetscInfo1(ksp, "KSPSolve_NASH: constrained step: radius=%g\n", (double)cg->radius);
416: if (norm_p > 0.0) {
417: /*********************************************************************/
418: /* Follow the direction to the boundary of the trust region. */
419: /*********************************************************************/
421: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
422: cg->norm_d = cg->radius;
424: VecAXPY(d, step, p); /* d = d + step p */
426: /*********************************************************************/
427: /* Update objective function. */
428: /*********************************************************************/
430: cg->o_fcn += step * (0.5 * step * kappa - rz);
431: } else {
432: /*********************************************************************/
433: /* The norm of the direction is zero; there is nothing to follow. */
434: /*********************************************************************/
435: }
436: break;
437: }
439: /*************************************************************************/
440: /* Now we can update the direction and residual. */
441: /*************************************************************************/
443: VecAXPY(d, alpha, p); /* d = d + alpha p */
444: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
445: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
447: switch (cg->dtype) {
448: case NASH_PRECONDITIONED_DIRECTION:
449: norm_d = norm_dp1;
450: break;
452: default:
453: VecDot(d, d, &norm_d);
454: break;
455: }
456: cg->norm_d = PetscSqrtReal(norm_d);
458: /*************************************************************************/
459: /* Update objective function. */
460: /*************************************************************************/
462: cg->o_fcn -= 0.5 * alpha * rz;
464: /*************************************************************************/
465: /* Check that the preconditioner appears positive semidefinite. */
466: /*************************************************************************/
468: rzm1 = rz;
469: VecDot(r, z, &rz); /* rz = r^T z */
470: if (rz < 0.0) {
471: /***********************************************************************/
472: /* The preconditioner is indefinite. */
473: /***********************************************************************/
475: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
476: PetscInfo1(ksp, "KSPSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz);
477: break;
478: }
480: /*************************************************************************/
481: /* As far as we know, the preconditioner is positive semidefinite. */
482: /* Compute the residual and check for convergence. */
483: /*************************************************************************/
485: switch (ksp->normtype) {
486: case KSP_NORM_PRECONDITIONED:
487: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
488: break;
490: case KSP_NORM_UNPRECONDITIONED:
491: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
492: break;
494: case KSP_NORM_NATURAL:
495: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
496: break;
498: default:
499: norm_r = 0.;
500: break;
501: }
503: KSPLogResidualHistory(ksp, norm_r);
504: KSPMonitor(ksp, ksp->its, norm_r);
505: ksp->rnorm = norm_r;
507: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
508: if (ksp->reason) {
509: /***********************************************************************/
510: /* The method has converged. */
511: /***********************************************************************/
513: PetscInfo2(ksp, "KSPSolve_NASH: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
514: break;
515: }
517: /*************************************************************************/
518: /* We have not converged yet. Check for breakdown. */
519: /*************************************************************************/
521: beta = rz / rzm1;
522: if (PetscAbsReal(beta) <= 0.0) {
523: /***********************************************************************/
524: /* Conjugate gradients has broken down. */
525: /***********************************************************************/
527: ksp->reason = KSP_DIVERGED_BREAKDOWN;
528: PetscInfo1(ksp, "KSPSolve_NASH: breakdown: beta=%g\n", (double)beta);
529: break;
530: }
532: /*************************************************************************/
533: /* Check iteration limit. */
534: /*************************************************************************/
536: if (ksp->its >= max_cg_its) {
537: ksp->reason = KSP_DIVERGED_ITS;
538: PetscInfo1(ksp, "KSPSolve_NASH: iterlim: its=%D\n", ksp->its);
539: break;
540: }
542: /*************************************************************************/
543: /* Update p and the norms. */
544: /*************************************************************************/
546: VecAYPX(p, beta, z); /* p = z + beta p */
548: switch (cg->dtype) {
549: case NASH_PRECONDITIONED_DIRECTION:
550: dMp = beta*(dMp + alpha*norm_p);
551: norm_p = beta*(rzm1 + beta*norm_p);
552: break;
554: default:
555: VecDot(d, p, &dMp);
556: VecDot(p, p, &norm_p);
557: break;
558: }
560: /*************************************************************************/
561: /* Compute the new direction and update the iteration. */
562: /*************************************************************************/
564: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
565: VecDot(p, z, &kappa); /* kappa = p^T Q p */
566: ++ksp->its;
568: /*************************************************************************/
569: /* Check for negative curvature. */
570: /*************************************************************************/
572: if (kappa <= 0.0) {
573: /***********************************************************************/
574: /* In this case, the matrix is indefinite and we have encountered */
575: /* a direction of negative curvature. Stop at the base. */
576: /***********************************************************************/
578: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
579: PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
580: break;
581: }
582: }
583: return(0);
584: #endif
585: }
589: PetscErrorCode KSPSetUp_NASH(KSP ksp)
590: {
594: /***************************************************************************/
595: /* Set work vectors needed by conjugate gradient method and allocate */
596: /***************************************************************************/
598: KSPSetWorkVecs(ksp, 3);
599: return(0);
600: }
604: PetscErrorCode KSPDestroy_NASH(KSP ksp)
605: {
609: /***************************************************************************/
610: /* Clear composed functions */
611: /***************************************************************************/
613: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",NULL);
614: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",NULL);
615: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",NULL);
617: /***************************************************************************/
618: /* Destroy KSP object. */
619: /***************************************************************************/
621: KSPDestroyDefault(ksp);
622: return(0);
623: }
627: static PetscErrorCode KSPNASHSetRadius_NASH(KSP ksp, PetscReal radius)
628: {
629: KSP_NASH *cg = (KSP_NASH*)ksp->data;
632: cg->radius = radius;
633: return(0);
634: }
638: static PetscErrorCode KSPNASHGetNormD_NASH(KSP ksp, PetscReal *norm_d)
639: {
640: KSP_NASH *cg = (KSP_NASH*)ksp->data;
643: *norm_d = cg->norm_d;
644: return(0);
645: }
649: static PetscErrorCode KSPNASHGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
650: {
651: KSP_NASH *cg = (KSP_NASH*)ksp->data;
654: *o_fcn = cg->o_fcn;
655: return(0);
656: }
660: PetscErrorCode KSPSetFromOptions_NASH(PetscOptionItems *PetscOptionsObject,KSP ksp)
661: {
663: KSP_NASH *cg = (KSP_NASH*)ksp->data;
666: PetscOptionsHead(PetscOptionsObject,"KSP NASH options");
668: PetscOptionsReal("-ksp_nash_radius", "Trust Region Radius", "KSPNASHSetRadius", cg->radius, &cg->radius, NULL);
670: PetscOptionsEList("-ksp_nash_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
672: PetscOptionsTail();
673: return(0);
674: }
676: /*MC
677: KSPNASH - Code to run conjugate gradient method subject to a constraint
678: on the solution norm. This is used in Trust Region methods for
679: nonlinear equations, SNESNEWTONTR
681: Options Database Keys:
682: . -ksp_nash_radius <r> - Trust Region Radius
684: Notes: This is rarely used directly
686: Level: developer
688: Use preconditioned conjugate gradient to compute
689: an approximate minimizer of the quadratic function
691: q(s) = g^T * s + 0.5 * s^T * H * s
693: subject to the trust region constraint
695: || s || <= delta,
697: where
699: delta is the trust region radius,
700: g is the gradient vector,
701: H is the Hessian approximation, and
702: M is the positive definite preconditioner matrix.
704: KSPConvergedReason may be
705: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
706: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
707: $ other KSP converged/diverged reasons
709: Notes:
710: The preconditioner supplied should be symmetric and positive definite.
712: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPNASHSetRadius(), KSPNASHGetNormD(), KSPNASHGetObjFcn()
713: M*/
717: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
718: {
720: KSP_NASH *cg;
723: PetscNewLog(ksp,&cg);
724: cg->radius = 0.0;
725: cg->dtype = NASH_UNPRECONDITIONED_DIRECTION;
727: ksp->data = (void*) cg;
728: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
729: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
730: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
732: /***************************************************************************/
733: /* Sets the functions that are associated with this data structure */
734: /* (in C++ this is the same as defining virtual functions). */
735: /***************************************************************************/
737: ksp->ops->setup = KSPSetUp_NASH;
738: ksp->ops->solve = KSPSolve_NASH;
739: ksp->ops->destroy = KSPDestroy_NASH;
740: ksp->ops->setfromoptions = KSPSetFromOptions_NASH;
741: ksp->ops->buildsolution = KSPBuildSolutionDefault;
742: ksp->ops->buildresidual = KSPBuildResidualDefault;
743: ksp->ops->view = 0;
745: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",KSPNASHSetRadius_NASH);
746: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",KSPNASHGetNormD_NASH);
747: PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",KSPNASHGetObjFcn_NASH);
748: return(0);
749: }