Actual source code: stcg.c
1: #define PETSCKSP_DLL
3: #include include/private/kspimpl.h
4: #include src/ksp/ksp/impls/cg/stcg/stcg.h
8: /*@
9: KSPSTCGSetRadius - Sets the radius of the trust region.
11: Collective on KSP
13: Input Parameters:
14: + ksp - the iterative context
15: - radius - the trust region radius (Infinity is the default)
17: Options Database Key:
18: . -ksp_stcg_radius <r>
20: Level: advanced
22: .keywords: KSP, STCG, set, trust region radius
23: @*/
24: PetscErrorCode KSPSTCGSetRadius(KSP ksp,PetscReal radius)
25: {
26: PetscErrorCode ierr, (*f)(KSP, PetscReal);
30: if (radius < 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be positive");
31: PetscObjectQueryFunction((PetscObject)ksp, "KSPSTCGSetRadius_C", (void (**)(void))&f);
32: if (f) {
33: (*f)(ksp, radius);
34: }
35: return(0);
36: }
40: /*@
41: KSPSTCGGetNormD - Got norm of the direction.
43: Collective on KSP
45: Input Parameters:
46: + ksp - the iterative context
47: - norm_d - the norm of the direction
49: Level: advanced
51: .keywords: KSP, STCG, get, norm direction
52: @*/
53: PetscErrorCode KSPSTCGGetNormD(KSP ksp,PetscReal *norm_d)
54: {
55: PetscErrorCode ierr, (*f)(KSP, PetscReal *);
59: PetscObjectQueryFunction((PetscObject)ksp, "KSPSTCGGetNormD_C", (void (**)(void))&f);
60: if (f) {
61: (*f)(ksp, norm_d);
62: }
63: return(0);
64: }
68: /*
69: KSPSolve_STCG - Use preconditioned conjugate gradient to compute
70: an approximate minimizer of the quadratic function
72: q(s) = g^T * s + .5 * s^T * H * s
74: subject to the trust region constraint
76: || s ||_M <= delta,
78: where
80: delta is the trust region radius,
81: g is the gradient vector, and
82: H is the Hessian matrix,
83: M is the positive definite preconditioner matrix.
85: KSPConvergedReason may be
86: $ KSP_CONVERGED_STCG_NEG_CURVE if convergence is reached along a negative curvature direction,
87: $ KSP_CONVERGED_STCG_CONSTRAINED if convergence is reached along a constrained step,
88: $ other KSP converged/diverged reasons
91: Notes:
92: The preconditioner supplied should be symmetric and positive definite.
93: */
95: #ifdef PETSC_USE_COMPLEX
96: PetscErrorCode STCG_VecDot(Vec x, Vec y, PetscReal *a)
97: {
98: PetscScalar ca;
101: VecDot(x,y,&ca);
102: *a = PetscRealPart(ca);
103: return ierr;
104: }
105: #else
106: PetscErrorCode STCG_VecDot(Vec x, Vec y, PetscReal *a)
107: {
108: return VecDot(x,y,a);
109: }
110: #endif
112: PetscErrorCode KSPSolve_STCG(KSP ksp)
113: {
115: MatStructure pflag;
116: Mat Qmat, Mmat;
117: Vec r, z, p, d;
118: PC pc;
119: KSP_STCG *cg;
120: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
121: PetscReal alpha, beta, kappa, rz, rzm1;
122: PetscReal r2;
123: PetscInt i, maxit;
124: PetscTruth diagonalscale;
127: PCDiagonalScale(ksp->pc, &diagonalscale);
128: if (diagonalscale) SETERRQ1(PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ksp->type_name);
130: cg = (KSP_STCG *)ksp->data;
131: r2 = cg->radius * cg->radius;
132: maxit = ksp->max_it;
133: r = ksp->work[0];
134: z = ksp->work[1];
135: p = ksp->work[2];
136: d = ksp->vec_sol;
137: pc = ksp->pc;
139: ksp->its = 0;
140: if (cg->radius < 0.0) {
141: SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
142: }
144: /* Initialize variables */
145: PCGetOperators(pc, &Qmat, &Mmat, &pflag);
147: VecSet(d, 0.0); /* d = 0 */
148: VecCopy(ksp->vec_rhs, r); /* r = -grad */
149: KSP_PCApply(ksp, r, z); /* z = M_{-1} r */
150: cg->norm_d = 0.0;
152: /* Check for numerical problems with the preconditioner */
153: STCG_VecDot(r, z, &rz); /* rz = r^T z */
154: if ((rz != rz) || (rz && (rz / rz != rz / rz))) {
155: ksp->reason = KSP_DIVERGED_NAN;
156: PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", rz);
158: /* In this case, the preconditioner produced not a number or an */
159: /* infinite value. We just take the gradient step. */
160: /* Only needs to be checked once. */
162: if (cg->radius) {
163: STCG_VecDot(r, r, &rz); /* rz = r^T r */
165: alpha = sqrt(r2 / rz);
166: VecAXPY(d, alpha, r); /* d = d + alpha r */
167: cg->norm_d = cg->radius;
168: }
169: return(0);
170: }
172: /* Check that the preconditioner is positive definite */
173: if (rz <= 0.0) {
174: VecNorm(r, NORM_2, &norm_r);
175: if (rz < 0.0 || norm_r > 0.0) {
176: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
177: PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", rz);
179: /* In this case, the preconditioner is indefinite, so we cannot measure */
180: /* the direction in the preconditioned norm. Therefore, we must use */
181: /* an unpreconditioned calculation. The direction in this case is */
182: /* uses the right hand side, which should be the negative gradient */
183: /* intersected with the trust region. */
185: if (cg->radius) {
186: STCG_VecDot(r, r, &rz); /* rz = r^T r */
187:
188: alpha = sqrt(r2 / rz);
189: VecAXPY(d, alpha, r); /* d = d + alpha r */
190: cg->norm_d = cg->radius;
191: }
192: return(0);
193: }
194: }
196: /* As far as we know, the preconditioner is positive definite. Compute */
197: /* the appropriate residual depending on what the user has set. */
198: if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
199: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
200: }
201: else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
202: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
203: }
204: else if (ksp->normtype == KSP_NATURAL_NORM) {
205: norm_r = sqrt(rz); /* norm_r = |r|_B */
206: }
207: else {
208: norm_r = 0;
209: }
211: /* Log the residual and call any registered moitor routines */
212: KSPLogResidualHistory(ksp, norm_r);
213: KSPMonitor(ksp, 0, norm_r);
214: ksp->rnorm = norm_r;
216: /* Test for convergence */
217: (*ksp->converged)(ksp, 0, norm_r, &ksp->reason, ksp->cnvP);
218: if (ksp->reason) {
219: return(0);
220: }
222: /* Compute the initial vectors and variables for trust-region computations */
223: VecCopy(z, p); /* p = z */
225: if (STCG_PRECONDITIONED_DIRECTION == cg->dtype) {
226: dMp = 0.0;
227: norm_p = rz;
228: norm_d = 0.0;
229: }
230: else {
231: dMp = 0.0;
232: STCG_VecDot(p, p, &norm_p);
233: norm_d = 0.0;
234: }
236: /* Compute the direction */
237: KSP_MatMult(ksp, Qmat, p, z); /* Qp = Q * p */
238: STCG_VecDot(p, z, &kappa); /* kappa = p^T Qp */
240: if ((kappa != kappa) || (kappa && (kappa / kappa != kappa / kappa))) {
241: ksp->reason = KSP_DIVERGED_NAN;
242: PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", kappa);
244: /* In this case, the matrix produced not a number or an infinite value. */
245: /* We just take the gradient step. Only needs to be checked once. */
246: if (cg->radius) {
247: STCG_VecDot(r, r, &rz); /* rz = r^T r */
249: alpha = sqrt(r2 / rz);
250: VecAXPY(d, alpha, r); /* d = d + alpha r */
251: cg->norm_d = cg->radius;
252: }
253: return(0);
254: }
256: /* Begin iterating */
257: for (i = 0; i <= maxit; i++) {
258: ++ksp->its;
260: /* Check for negative curvature */
261: if (kappa <= 0.0) {
262: ksp->reason = KSP_CONVERGED_STCG_NEG_CURVE;
263: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", kappa);
265: /* In this case, the matrix is indefinite and we have encountered */
266: /* a direction of negative curvature. Follow the direction to the */
267: /* boundary of the trust region. */
268: if (cg->radius) {
269: alpha = (sqrt(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
270: VecAXPY(d, alpha, p); /* d = d + alpha p */
271: cg->norm_d = cg->radius;
272: }
273: break;
274: }
275: alpha = rz / kappa;
277: /* Now we can update the direction and residual. This is perhaps not */
278: /* the best order to perform the breakdown checks, but the optimization */
279: /* codes need a direction and this is the best we can do. */
281: /* First test if the new direction intersects the trust region. */
282: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
283: if (cg->radius && norm_dp1 >= r2) {
284: ksp->reason = KSP_CONVERGED_STCG_CONSTRAINED;
285: PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", cg->radius);
287: /* In this case, the matrix is positive definite as far as we know. */
288: /* However, the direction does beyond the trust region. Follow the */
289: /* direction to the boundary of the trust region. */
291: alpha = (sqrt(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
292: VecAXPY(d, alpha, p); /* d = d + alpha p */
293: cg->norm_d = cg->radius;
294: break;
295: }
297: /* Update the direction and residual */
298: VecAXPY(d, alpha, p); /* d = d + alpha p */
299: VecAXPY(r, -alpha, z); /* r = r - alpha Qp */
300: KSP_PCApply(ksp, r, z);
302: if (STCG_PRECONDITIONED_DIRECTION == cg->dtype) {
303: norm_d = norm_dp1;
304: }
305: else {
306: STCG_VecDot(d, d, &norm_d);
307: }
308: cg->norm_d = sqrt(norm_d);
310: /* Check that the preconditioner is positive semidefinite */
311: rzm1 = rz;
312: STCG_VecDot(r, z, &rz); /* rz = r^T z */
313: if (rz <= 0.0) {
314: VecNorm(r, NORM_2, &norm_r);
315: if (rz < 0.0 || norm_r > 0.0) {
316: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
317: PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", rz);
319: /* In this case, the preconditioner is indefinite. We could follow */
320: /* the direction to the boundary of the trust region, but it seems */
321: /* best to stop at the current point. */
322: break;
323: }
324: }
326: /* As far as we know, the matrix and preconditioner are positive */
327: /* definite. Compute the appropriate residual depending on what the */
328: /* user has set. */
329: if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
330: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
331: }
332: else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
333: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
334: }
335: else if (ksp->normtype == KSP_NATURAL_NORM) {
336: norm_r = sqrt(rz); /* norm_r = |r|_B */
337: }
338: else {
339: norm_r = 0;
340: }
342: /* Log the residual and call any registered moitor routines */
343: KSPLogResidualHistory(ksp, norm_r);
344: KSPMonitor(ksp, 0, norm_r);
345: ksp->rnorm = norm_r;
346:
347: /* Test for convergence */
348: (*ksp->converged)(ksp, i+1, norm_r, &ksp->reason, ksp->cnvP);
349: if (ksp->reason) { /* convergence */
350: PetscInfo2(ksp,"KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n",norm_r,cg->radius);
351: break;
352: }
354: /* Update p and the norms */
355: beta = rz / rzm1;
356: VecAYPX(p, beta, z); /* p = z + beta p */
358: if (STCG_PRECONDITIONED_DIRECTION == cg->dtype) {
359: dMp = beta*(dMp + alpha*norm_p);
360: norm_p = beta*(rzm1 + beta*norm_p);
361: }
362: else {
363: STCG_VecDot(d, p, &dMp);
364: STCG_VecDot(p, p, &norm_p);
365: }
367: /* Compute new direction */
368: KSP_MatMult(ksp, Qmat, p, z); /* Qp = Q * p */
369: STCG_VecDot(p, z, &kappa); /* kappa = p^T Qp */
370: }
372: if (!ksp->reason) {
373: ksp->reason = KSP_DIVERGED_ITS;
374: }
375: return(0);
376: }
380: PetscErrorCode KSPSetUp_STCG(KSP ksp)
381: {
385: /* This implementation of CG only handles left preconditioning
386: * so generate an error otherwise.
387: */
388: if (ksp->pc_side == PC_RIGHT) {
389: SETERRQ(PETSC_ERR_SUP, "No right preconditioning for KSPSTCG");
390: }
391: else if (ksp->pc_side == PC_SYMMETRIC) {
392: SETERRQ(PETSC_ERR_SUP, "No symmetric preconditioning for KSPSTCG");
393: }
395: /* get work vectors needed by CG */
396: KSPDefaultGetWork(ksp, 3);
397: return(0);
398: }
402: PetscErrorCode KSPDestroy_STCG(KSP ksp)
403: {
404: KSP_STCG *cg = (KSP_STCG *)ksp->data;
408: KSPDefaultFreeWork(ksp);
410: /* Free the context variable */
411: PetscFree(cg);
412: return(0);
413: }
418: PetscErrorCode KSPSTCGSetRadius_STCG(KSP ksp,PetscReal radius)
419: {
420: KSP_STCG *cg = (KSP_STCG *)ksp->data;
423: cg->radius = radius;
424: return(0);
425: }
429: PetscErrorCode KSPSTCGGetNormD_STCG(KSP ksp,PetscReal *norm_d)
430: {
431: KSP_STCG *cg = (KSP_STCG *)ksp->data;
434: *norm_d = cg->norm_d;
435: return(0);
436: }
439: static const char *DType_Table[64] = {
440: "preconditioned", "unpreconditioned"
441: };
445: PetscErrorCode KSPSetFromOptions_STCG(KSP ksp)
446: {
448: KSP_STCG *cg = (KSP_STCG *)ksp->data;
451: PetscOptionsHead("KSP STCG options");
452: PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, PETSC_NULL);
453: PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, PETSC_NULL);
454: PetscOptionsTail();
455: return(0);
456: }
458: /*MC
459: KSPSTCG - Code to run conjugate gradient method subject to a constraint
460: on the solution norm. This is used in Trust Region methods for
461: nonlinear equations, SNESTR
463: Options Database Keys:
464: . -ksp_stcg_radius <r> - Trust Region Radius
466: Notes: This is rarely used directly
468: Level: developer
470: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD()
471: M*/
476: PetscErrorCode KSPCreate_STCG(KSP ksp)
477: {
479: KSP_STCG *cg;
482: PetscNew(KSP_STCG, &cg);
483: PetscLogObjectMemory(ksp, sizeof(KSP_STCG));
484: cg->radius = PETSC_MAX;
485: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
486: ksp->data = (void *)cg;
487: ksp->pc_side = PC_LEFT;
489: /* Sets the functions that are associated with this data structure
490: * (in C++ this is the same as defining virtual functions)
491: */
492: ksp->ops->setup = KSPSetUp_STCG;
493: ksp->ops->solve = KSPSolve_STCG;
494: ksp->ops->destroy = KSPDestroy_STCG;
495: ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
496: ksp->ops->buildsolution = KSPDefaultBuildSolution;
497: ksp->ops->buildresidual = KSPDefaultBuildResidual;
498: ksp->ops->view = 0;
500: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGSetRadius_C",
501: "KSPSTCGSetRadius_STCG",
502: KSPSTCGSetRadius_STCG);
503: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGGetNormD_C",
504: "KSPSTCGGetNormD_STCG",
505: KSPSTCGGetNormD_STCG);
506: return(0);
507: }