Actual source code: stcg.c
petsc-3.9.4 2018-09-11
2: #include <../src/ksp/ksp/impls/cg/stcg/stcgimpl.h>
4: #define STCG_PRECONDITIONED_DIRECTION 0
5: #define STCG_UNPRECONDITIONED_DIRECTION 1
6: #define STCG_DIRECTION_TYPES 2
8: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
10: static PetscErrorCode KSPCGSolve_STCG(KSP ksp)
11: {
12: #if defined(PETSC_USE_COMPLEX)
13: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "STCG is not available for complex systems");
14: #else
15: KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;
17: Mat Qmat, Mmat;
18: Vec r, z, p, d;
19: 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;
23: PetscInt max_cg_its;
24: PetscBool diagonalscale;
26: /***************************************************************************/
27: /* Check the arguments and parameters. */
28: /***************************************************************************/
31: PCGetDiagonalScale(ksp->pc, &diagonalscale);
32: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
33: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
35: /***************************************************************************/
36: /* Get the workspace vectors and initialize variables */
37: /***************************************************************************/
39: r2 = cg->radius * cg->radius;
40: r = ksp->work[0];
41: z = ksp->work[1];
42: p = ksp->work[2];
43: d = ksp->vec_sol;
44: pc = ksp->pc;
46: PCGetOperators(pc, &Qmat, &Mmat);
48: VecGetSize(d, &max_cg_its);
49: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
50: ksp->its = 0;
52: /***************************************************************************/
53: /* Initialize objective function and direction. */
54: /***************************************************************************/
56: cg->o_fcn = 0.0;
58: VecSet(d, 0.0); /* d = 0 */
59: cg->norm_d = 0.0;
61: /***************************************************************************/
62: /* Begin the conjugate gradient method. Check the right-hand side for */
63: /* numerical problems. The check for not-a-number and infinite values */
64: /* need be performed only once. */
65: /***************************************************************************/
67: VecCopy(ksp->vec_rhs, r); /* r = -grad */
68: VecDot(r, r, &rr); /* rr = r^T r */
69: KSPCheckDot(ksp,rr);
71: /***************************************************************************/
72: /* Check the preconditioner for numerical problems and for positive */
73: /* definiteness. The check for not-a-number and infinite values need be */
74: /* performed only once. */
75: /***************************************************************************/
77: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
78: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
79: if (PetscIsInfOrNanScalar(rz)) {
80: /*************************************************************************/
81: /* The preconditioner contains not-a-number or an infinite value. */
82: /* Return the gradient direction intersected with the trust region. */
83: /*************************************************************************/
85: ksp->reason = KSP_DIVERGED_NANORINF;
86: PetscInfo1(ksp, "KSPCGSolve_STCG: bad preconditioner: rz=%g\n", (double)rz);
88: if (cg->radius != 0) {
89: if (r2 >= rr) {
90: alpha = 1.0;
91: cg->norm_d = PetscSqrtReal(rr);
92: } else {
93: alpha = PetscSqrtReal(r2 / rr);
94: cg->norm_d = cg->radius;
95: }
97: VecAXPY(d, alpha, r); /* d = d + alpha r */
99: /***********************************************************************/
100: /* Compute objective function. */
101: /***********************************************************************/
103: KSP_MatMult(ksp, Qmat, d, z);
104: VecAYPX(z, -0.5, ksp->vec_rhs);
105: VecDot(d, z, &cg->o_fcn);
106: cg->o_fcn = -cg->o_fcn;
107: ++ksp->its;
108: }
109: return(0);
110: }
112: if (rz < 0.0) {
113: /*************************************************************************/
114: /* The preconditioner is indefinite. Because this is the first */
115: /* and we do not have a direction yet, we use the gradient step. Note */
116: /* that we cannot use the preconditioned norm when computing the step */
117: /* because the matrix is indefinite. */
118: /*************************************************************************/
120: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
121: PetscInfo1(ksp, "KSPCGSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz);
123: if (cg->radius != 0.0) {
124: if (r2 >= rr) {
125: alpha = 1.0;
126: cg->norm_d = PetscSqrtReal(rr);
127: } else {
128: alpha = PetscSqrtReal(r2 / rr);
129: cg->norm_d = cg->radius;
130: }
132: VecAXPY(d, alpha, r); /* d = d + alpha r */
134: /***********************************************************************/
135: /* Compute objective function. */
136: /***********************************************************************/
138: KSP_MatMult(ksp, Qmat, d, z);
139: VecAYPX(z, -0.5, ksp->vec_rhs);
140: VecDot(d, z, &cg->o_fcn);
141: cg->o_fcn = -cg->o_fcn;
142: ++ksp->its;
143: }
144: return(0);
145: }
147: /***************************************************************************/
148: /* As far as we know, the preconditioner is positive semidefinite. */
149: /* Compute and log the residual. Check convergence because this */
150: /* initializes things, but do not terminate until at least one conjugate */
151: /* gradient iteration has been performed. */
152: /***************************************************************************/
154: switch (ksp->normtype) {
155: case KSP_NORM_PRECONDITIONED:
156: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
157: break;
159: case KSP_NORM_UNPRECONDITIONED:
160: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
161: break;
163: case KSP_NORM_NATURAL:
164: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
165: break;
167: default:
168: norm_r = 0.0;
169: break;
170: }
172: KSPLogResidualHistory(ksp, norm_r);
173: KSPMonitor(ksp, ksp->its, norm_r);
174: ksp->rnorm = norm_r;
176: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
178: /***************************************************************************/
179: /* Compute the first direction and update the iteration. */
180: /***************************************************************************/
182: VecCopy(z, p); /* p = z */
183: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
184: ++ksp->its;
186: /***************************************************************************/
187: /* Check the matrix for numerical problems. */
188: /***************************************************************************/
190: VecDot(p, z, &kappa); /* kappa = p^T Q p */
191: if (PetscIsInfOrNanScalar(kappa)) {
192: /*************************************************************************/
193: /* The matrix produced not-a-number or an infinite value. In this case, */
194: /* we must stop and use the gradient direction. This condition need */
195: /* only be checked once. */
196: /*************************************************************************/
198: ksp->reason = KSP_DIVERGED_NANORINF;
199: PetscInfo1(ksp, "KSPCGSolve_STCG: bad matrix: kappa=%g\n", (double)kappa);
201: if (cg->radius) {
202: if (r2 >= rr) {
203: alpha = 1.0;
204: cg->norm_d = PetscSqrtReal(rr);
205: } else {
206: alpha = PetscSqrtReal(r2 / rr);
207: cg->norm_d = cg->radius;
208: }
210: VecAXPY(d, alpha, r); /* d = d + alpha r */
212: /***********************************************************************/
213: /* Compute objective function. */
214: /***********************************************************************/
216: KSP_MatMult(ksp, Qmat, d, z);
217: VecAYPX(z, -0.5, ksp->vec_rhs);
218: VecDot(d, z, &cg->o_fcn);
219: cg->o_fcn = -cg->o_fcn;
220: ++ksp->its;
221: }
222: return(0);
223: }
225: /***************************************************************************/
226: /* Initialize variables for calculating the norm of the direction. */
227: /***************************************************************************/
229: dMp = 0.0;
230: norm_d = 0.0;
231: switch (cg->dtype) {
232: case STCG_PRECONDITIONED_DIRECTION:
233: norm_p = rz;
234: break;
236: default:
237: VecDot(p, p, &norm_p);
238: break;
239: }
241: /***************************************************************************/
242: /* Check for negative curvature. */
243: /***************************************************************************/
245: if (kappa <= 0.0) {
246: /*************************************************************************/
247: /* In this case, the matrix is indefinite and we have encountered a */
248: /* direction of negative curvature. Because negative curvature occurs */
249: /* during the first step, we must follow a direction. */
250: /*************************************************************************/
252: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
253: PetscInfo1(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);
255: if (cg->radius != 0.0 && norm_p > 0.0) {
256: /***********************************************************************/
257: /* Follow direction of negative curvature to the boundary of the */
258: /* trust region. */
259: /***********************************************************************/
261: step = PetscSqrtReal(r2 / norm_p);
262: cg->norm_d = cg->radius;
264: VecAXPY(d, step, p); /* d = d + step p */
266: /***********************************************************************/
267: /* Update objective function. */
268: /***********************************************************************/
270: cg->o_fcn += step * (0.5 * step * kappa - rz);
271: } else if (cg->radius != 0.0) {
272: /***********************************************************************/
273: /* The norm of the preconditioned direction is zero; use the gradient */
274: /* step. */
275: /***********************************************************************/
277: if (r2 >= rr) {
278: alpha = 1.0;
279: cg->norm_d = PetscSqrtReal(rr);
280: } else {
281: alpha = PetscSqrtReal(r2 / rr);
282: cg->norm_d = cg->radius;
283: }
285: VecAXPY(d, alpha, r); /* d = d + alpha r */
287: /***********************************************************************/
288: /* Compute objective function. */
289: /***********************************************************************/
291: KSP_MatMult(ksp, Qmat, d, z);
292: VecAYPX(z, -0.5, ksp->vec_rhs);
293: 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 != 0.0 && 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: PetscInfo1(ksp, "KSPCGSolve_STCG: 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 STCG_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: PetscInfo1(ksp, "KSPCGSolve_STCG: 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.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: PetscInfo2(ksp, "KSPCGSolve_STCG: 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 (PetscAbsScalar(beta) <= 0.0) {
437: /***********************************************************************/
438: /* Conjugate gradients has broken down. */
439: /***********************************************************************/
441: ksp->reason = KSP_DIVERGED_BREAKDOWN;
442: PetscInfo1(ksp, "KSPCGSolve_STCG: 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: PetscInfo1(ksp, "KSPCGSolve_STCG: iterlim: its=%D\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 STCG_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. Follow the direction to the */
490: /* boundary of the trust region. */
491: /***********************************************************************/
493: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
494: PetscInfo1(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);
496: if (cg->radius != 0.0 && norm_p > 0.0) {
497: /*********************************************************************/
498: /* Follow direction of negative curvature to boundary. */
499: /*********************************************************************/
501: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
502: cg->norm_d = cg->radius;
504: VecAXPY(d, step, p); /* d = d + step p */
506: /*********************************************************************/
507: /* Update objective function. */
508: /*********************************************************************/
510: cg->o_fcn += step * (0.5 * step * kappa - rz);
511: } else if (cg->radius != 0.0) {
512: /*********************************************************************/
513: /* The norm of the direction is zero; there is nothing to follow. */
514: /*********************************************************************/
515: }
516: break;
517: }
518: }
519: return(0);
520: #endif
521: }
523: static PetscErrorCode KSPCGSetUp_STCG(KSP ksp)
524: {
528: /***************************************************************************/
529: /* Set work vectors needed by conjugate gradient method and allocate */
530: /***************************************************************************/
532: KSPSetWorkVecs(ksp, 3);
533: return(0);
534: }
536: static PetscErrorCode KSPCGDestroy_STCG(KSP ksp)
537: {
541: /***************************************************************************/
542: /* Clear composed functions */
543: /***************************************************************************/
545: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
546: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
547: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
549: /***************************************************************************/
550: /* Destroy KSP object. */
551: /***************************************************************************/
553: KSPDestroyDefault(ksp);
554: return(0);
555: }
557: static PetscErrorCode KSPCGSetRadius_STCG(KSP ksp, PetscReal radius)
558: {
559: KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;
562: cg->radius = radius;
563: return(0);
564: }
566: static PetscErrorCode KSPCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
567: {
568: KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;
571: *norm_d = cg->norm_d;
572: return(0);
573: }
575: static PetscErrorCode KSPCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
576: {
577: KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;
580: *o_fcn = cg->o_fcn;
581: return(0);
582: }
584: static PetscErrorCode KSPCGSetFromOptions_STCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
585: {
587: KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;
590: PetscOptionsHead(PetscOptionsObject,"KSPCG STCG options");
591: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
592: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
593: PetscOptionsTail();
594: return(0);
595: }
597: /*MC
598: KSPCGSTCG - Code to run conjugate gradient method subject to a constraint
599: on the solution norm. This is used in Trust Region methods for
600: nonlinear equations, SNESNEWTONTR
602: Options Database Keys:
603: . -ksp_cg_radius <r> - Trust Region Radius
605: Notes: This is rarely used directly
607: Use preconditioned conjugate gradient to compute
608: an approximate minimizer of the quadratic function
610: q(s) = g^T * s + 0.5 * s^T * H * s
612: subject to the trust region constraint
614: || s || <= delta,
616: where
618: delta is the trust region radius,
619: g is the gradient vector,
620: H is the Hessian approximation, and
621: M is the positive definite preconditioner matrix.
623: KSPConvergedReason may be
624: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
625: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
626: $ other KSP converged/diverged reasons
628: Notes:
629: The preconditioner supplied should be symmetric and positive definite.
631: Level: developer
633: .seealso: KSPCreate(), KSPCGSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn()
634: M*/
636: PETSC_EXTERN PetscErrorCode KSPCreate_CGSTCG(KSP ksp)
637: {
639: KSPCG_STCG *cg;
642: PetscNewLog(ksp,&cg);
644: cg->radius = 0.0;
645: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
647: ksp->data = (void*) cg;
648: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
649: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
650: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
652: /***************************************************************************/
653: /* Sets the functions that are associated with this data structure */
654: /* (in C++ this is the same as defining virtual functions). */
655: /***************************************************************************/
657: ksp->ops->setup = KSPCGSetUp_STCG;
658: ksp->ops->solve = KSPCGSolve_STCG;
659: ksp->ops->destroy = KSPCGDestroy_STCG;
660: ksp->ops->setfromoptions = KSPCGSetFromOptions_STCG;
661: ksp->ops->buildsolution = KSPBuildSolutionDefault;
662: ksp->ops->buildresidual = KSPBuildResidualDefault;
663: ksp->ops->view = 0;
665: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_STCG);
666: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",KSPCGGetNormD_STCG);
667: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_STCG);
668: return(0);
669: }