Actual source code: stcg.c
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;
16: Mat Qmat, Mmat;
17: Vec r, z, p, d;
18: PC pc;
19: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
20: PetscReal alpha, beta, kappa, rz, rzm1;
21: PetscReal rr, r2, step;
22: PetscInt max_cg_its;
23: PetscBool diagonalscale;
25: /***************************************************************************/
26: /* Check the arguments and parameters. */
27: /***************************************************************************/
29: PetscFunctionBegin;
30: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
31: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
32: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
34: /***************************************************************************/
35: /* Get the workspace vectors and initialize variables */
36: /***************************************************************************/
38: r2 = cg->radius * cg->radius;
39: r = ksp->work[0];
40: z = ksp->work[1];
41: p = ksp->work[2];
42: d = ksp->vec_sol;
43: pc = ksp->pc;
45: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
47: PetscCall(VecGetSize(d, &max_cg_its));
48: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
49: ksp->its = 0;
51: /***************************************************************************/
52: /* Initialize objective function and direction. */
53: /***************************************************************************/
55: cg->o_fcn = 0.0;
57: PetscCall(VecSet(d, 0.0)); /* d = 0 */
58: cg->norm_d = 0.0;
60: /***************************************************************************/
61: /* Begin the conjugate gradient method. Check the right-hand side for */
62: /* numerical problems. The check for not-a-number and infinite values */
63: /* need be performed only once. */
64: /***************************************************************************/
66: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
67: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
68: KSPCheckDot(ksp, rr);
70: /***************************************************************************/
71: /* Check the preconditioner for numerical problems and for positive */
72: /* definiteness. The check for not-a-number and infinite values need be */
73: /* performed only once. */
74: /***************************************************************************/
76: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
77: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
78: if (PetscIsInfOrNanScalar(rz)) {
79: /*************************************************************************/
80: /* The preconditioner contains not-a-number or an infinite value. */
81: /* Return the gradient direction intersected with the trust region. */
82: /*************************************************************************/
84: ksp->reason = KSP_DIVERGED_NANORINF;
85: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad preconditioner: rz=%g\n", (double)rz));
87: if (cg->radius != 0) {
88: if (r2 >= rr) {
89: alpha = 1.0;
90: cg->norm_d = PetscSqrtReal(rr);
91: } else {
92: alpha = PetscSqrtReal(r2 / rr);
93: cg->norm_d = cg->radius;
94: }
96: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
98: /***********************************************************************/
99: /* Compute objective function. */
100: /***********************************************************************/
102: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
103: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
104: PetscCall(VecDot(d, z, &cg->o_fcn));
105: cg->o_fcn = -cg->o_fcn;
106: ++ksp->its;
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: if (rz < 0.0) {
112: /*************************************************************************/
113: /* The preconditioner is indefinite. Because this is the first */
114: /* and we do not have a direction yet, we use the gradient step. Note */
115: /* that we cannot use the preconditioned norm when computing the step */
116: /* because the matrix is indefinite. */
117: /*************************************************************************/
119: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
120: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz));
122: if (cg->radius != 0.0) {
123: if (r2 >= rr) {
124: alpha = 1.0;
125: cg->norm_d = PetscSqrtReal(rr);
126: } else {
127: alpha = PetscSqrtReal(r2 / rr);
128: cg->norm_d = cg->radius;
129: }
131: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
133: /***********************************************************************/
134: /* Compute objective function. */
135: /***********************************************************************/
137: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
138: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
139: PetscCall(VecDot(d, z, &cg->o_fcn));
140: cg->o_fcn = -cg->o_fcn;
141: ++ksp->its;
142: }
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /***************************************************************************/
147: /* As far as we know, the preconditioner is positive semidefinite. */
148: /* Compute and log the residual. Check convergence because this */
149: /* initializes things, but do not terminate until at least one conjugate */
150: /* gradient iteration has been performed. */
151: /***************************************************************************/
153: switch (ksp->normtype) {
154: case KSP_NORM_PRECONDITIONED:
155: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
156: break;
158: case KSP_NORM_UNPRECONDITIONED:
159: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
160: break;
162: case KSP_NORM_NATURAL:
163: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
164: break;
166: default:
167: norm_r = 0.0;
168: break;
169: }
171: PetscCall(KSPLogResidualHistory(ksp, norm_r));
172: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
173: ksp->rnorm = norm_r;
175: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
177: /***************************************************************************/
178: /* Compute the first direction and update the iteration. */
179: /***************************************************************************/
181: PetscCall(VecCopy(z, p)); /* p = z */
182: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
183: ++ksp->its;
185: /***************************************************************************/
186: /* Check the matrix for numerical problems. */
187: /***************************************************************************/
189: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
190: if (PetscIsInfOrNanScalar(kappa)) {
191: /*************************************************************************/
192: /* The matrix produced not-a-number or an infinite value. In this case, */
193: /* we must stop and use the gradient direction. This condition need */
194: /* only be checked once. */
195: /*************************************************************************/
197: ksp->reason = KSP_DIVERGED_NANORINF;
198: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad matrix: kappa=%g\n", (double)kappa));
200: if (cg->radius) {
201: if (r2 >= rr) {
202: alpha = 1.0;
203: cg->norm_d = PetscSqrtReal(rr);
204: } else {
205: alpha = PetscSqrtReal(r2 / rr);
206: cg->norm_d = cg->radius;
207: }
209: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
211: /***********************************************************************/
212: /* Compute objective function. */
213: /***********************************************************************/
215: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
216: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
217: PetscCall(VecDot(d, z, &cg->o_fcn));
218: cg->o_fcn = -cg->o_fcn;
219: ++ksp->its;
220: }
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /***************************************************************************/
225: /* Initialize variables for calculating the norm of the direction. */
226: /***************************************************************************/
228: dMp = 0.0;
229: norm_d = 0.0;
230: switch (cg->dtype) {
231: case STCG_PRECONDITIONED_DIRECTION:
232: norm_p = rz;
233: break;
235: default:
236: PetscCall(VecDot(p, p, &norm_p));
237: break;
238: }
240: /***************************************************************************/
241: /* Check for negative curvature. */
242: /***************************************************************************/
244: if (kappa <= 0.0) {
245: /*************************************************************************/
246: /* In this case, the matrix is indefinite and we have encountered a */
247: /* direction of negative curvature. Because negative curvature occurs */
248: /* during the first step, we must follow a direction. */
249: /*************************************************************************/
251: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
252: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));
254: if (cg->radius != 0.0 && norm_p > 0.0) {
255: /***********************************************************************/
256: /* Follow direction of negative curvature to the boundary of the */
257: /* trust region. */
258: /***********************************************************************/
260: step = PetscSqrtReal(r2 / norm_p);
261: cg->norm_d = cg->radius;
263: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
265: /***********************************************************************/
266: /* Update objective function. */
267: /***********************************************************************/
269: cg->o_fcn += step * (0.5 * step * kappa - rz);
270: } else if (cg->radius != 0.0) {
271: /***********************************************************************/
272: /* The norm of the preconditioned direction is zero; use the gradient */
273: /* step. */
274: /***********************************************************************/
276: if (r2 >= rr) {
277: alpha = 1.0;
278: cg->norm_d = PetscSqrtReal(rr);
279: } else {
280: alpha = PetscSqrtReal(r2 / rr);
281: cg->norm_d = cg->radius;
282: }
284: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
286: /***********************************************************************/
287: /* Compute objective function. */
288: /***********************************************************************/
290: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
291: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
292: PetscCall(VecDot(d, z, &cg->o_fcn));
294: cg->o_fcn = -cg->o_fcn;
295: ++ksp->its;
296: }
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /***************************************************************************/
301: /* Run the conjugate gradient method until either the problem is solved, */
302: /* we encounter the boundary of the trust region, or the conjugate */
303: /* gradient method breaks down. */
304: /***************************************************************************/
306: while (1) {
307: /*************************************************************************/
308: /* Know that kappa is nonzero, because we have not broken down, so we */
309: /* can compute the steplength. */
310: /*************************************************************************/
312: alpha = rz / kappa;
314: /*************************************************************************/
315: /* Compute the steplength and check for intersection with the trust */
316: /* region. */
317: /*************************************************************************/
319: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
320: if (cg->radius != 0.0 && norm_dp1 >= r2) {
321: /***********************************************************************/
322: /* In this case, the matrix is positive definite as far as we know. */
323: /* However, the full step goes beyond the trust region. */
324: /***********************************************************************/
326: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
327: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: constrained step: radius=%g\n", (double)cg->radius));
329: if (norm_p > 0.0) {
330: /*********************************************************************/
331: /* Follow the direction to the boundary of the trust region. */
332: /*********************************************************************/
334: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
335: cg->norm_d = cg->radius;
337: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
339: /*********************************************************************/
340: /* Update objective function. */
341: /*********************************************************************/
343: cg->o_fcn += step * (0.5 * step * kappa - rz);
344: } else {
345: /*********************************************************************/
346: /* The norm of the direction is zero; there is nothing to follow. */
347: /*********************************************************************/
348: }
349: break;
350: }
352: /*************************************************************************/
353: /* Now we can update the direction and residual. */
354: /*************************************************************************/
356: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
357: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
358: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
360: switch (cg->dtype) {
361: case STCG_PRECONDITIONED_DIRECTION:
362: norm_d = norm_dp1;
363: break;
365: default:
366: PetscCall(VecDot(d, d, &norm_d));
367: break;
368: }
369: cg->norm_d = PetscSqrtReal(norm_d);
371: /*************************************************************************/
372: /* Update objective function. */
373: /*************************************************************************/
375: cg->o_fcn -= 0.5 * alpha * rz;
377: /*************************************************************************/
378: /* Check that the preconditioner appears positive semidefinite. */
379: /*************************************************************************/
381: rzm1 = rz;
382: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
383: if (rz < 0.0) {
384: /***********************************************************************/
385: /* The preconditioner is indefinite. */
386: /***********************************************************************/
388: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
389: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz));
390: break;
391: }
393: /*************************************************************************/
394: /* As far as we know, the preconditioner is positive semidefinite. */
395: /* Compute the residual and check for convergence. */
396: /*************************************************************************/
398: switch (ksp->normtype) {
399: case KSP_NORM_PRECONDITIONED:
400: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
401: break;
403: case KSP_NORM_UNPRECONDITIONED:
404: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
405: break;
407: case KSP_NORM_NATURAL:
408: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
409: break;
411: default:
412: norm_r = 0.0;
413: break;
414: }
416: PetscCall(KSPLogResidualHistory(ksp, norm_r));
417: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
418: ksp->rnorm = norm_r;
420: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
421: if (ksp->reason) {
422: /***********************************************************************/
423: /* The method has converged. */
424: /***********************************************************************/
426: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
427: break;
428: }
430: /*************************************************************************/
431: /* We have not converged yet. Check for breakdown. */
432: /*************************************************************************/
434: beta = rz / rzm1;
435: if (PetscAbsScalar(beta) <= 0.0) {
436: /***********************************************************************/
437: /* Conjugate gradients has broken down. */
438: /***********************************************************************/
440: ksp->reason = KSP_DIVERGED_BREAKDOWN;
441: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: breakdown: beta=%g\n", (double)beta));
442: break;
443: }
445: /*************************************************************************/
446: /* Check iteration limit. */
447: /*************************************************************************/
449: if (ksp->its >= max_cg_its) {
450: ksp->reason = KSP_DIVERGED_ITS;
451: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
452: break;
453: }
455: /*************************************************************************/
456: /* Update p and the norms. */
457: /*************************************************************************/
459: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
461: switch (cg->dtype) {
462: case STCG_PRECONDITIONED_DIRECTION:
463: dMp = beta * (dMp + alpha * norm_p);
464: norm_p = beta * (rzm1 + beta * norm_p);
465: break;
467: default:
468: PetscCall(VecDot(d, p, &dMp));
469: PetscCall(VecDot(p, p, &norm_p));
470: break;
471: }
473: /*************************************************************************/
474: /* Compute the new direction and update the iteration. */
475: /*************************************************************************/
477: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
478: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
479: ++ksp->its;
481: /*************************************************************************/
482: /* Check for negative curvature. */
483: /*************************************************************************/
485: if (kappa <= 0.0) {
486: /***********************************************************************/
487: /* In this case, the matrix is indefinite and we have encountered */
488: /* a direction of negative curvature. Follow the direction to the */
489: /* boundary of the trust region. */
490: /***********************************************************************/
492: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
493: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));
495: if (cg->radius != 0.0 && norm_p > 0.0) {
496: /*********************************************************************/
497: /* Follow direction of negative curvature to boundary. */
498: /*********************************************************************/
500: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
501: cg->norm_d = cg->radius;
503: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
505: /*********************************************************************/
506: /* Update objective function. */
507: /*********************************************************************/
509: cg->o_fcn += step * (0.5 * step * kappa - rz);
510: } else if (cg->radius != 0.0) {
511: /*********************************************************************/
512: /* The norm of the direction is zero; there is nothing to follow. */
513: /*********************************************************************/
514: }
515: break;
516: }
517: }
518: PetscFunctionReturn(PETSC_SUCCESS);
519: #endif
520: }
522: static PetscErrorCode KSPCGSetUp_STCG(KSP ksp)
523: {
524: PetscFunctionBegin;
525: /***************************************************************************/
526: /* Set work vectors needed by conjugate gradient method and allocate */
527: /***************************************************************************/
529: PetscCall(KSPSetWorkVecs(ksp, 3));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: static PetscErrorCode KSPCGDestroy_STCG(KSP ksp)
534: {
535: PetscFunctionBegin;
536: /***************************************************************************/
537: /* Clear composed functions */
538: /***************************************************************************/
540: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
541: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
542: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
544: /***************************************************************************/
545: /* Destroy KSP object. */
546: /***************************************************************************/
548: PetscCall(KSPDestroyDefault(ksp));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode KSPCGSetRadius_STCG(KSP ksp, PetscReal radius)
553: {
554: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
556: PetscFunctionBegin;
557: cg->radius = radius;
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode KSPCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
562: {
563: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
565: PetscFunctionBegin;
566: *norm_d = cg->norm_d;
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: static PetscErrorCode KSPCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
571: {
572: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
574: PetscFunctionBegin;
575: *o_fcn = cg->o_fcn;
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode KSPCGSetFromOptions_STCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
580: {
581: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
583: PetscFunctionBegin;
584: PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG STCG options");
585: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
586: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
587: PetscOptionsHeadEnd();
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: /*MC
592: KSPSTCG - Code to run conjugate gradient method subject to a constraint on the solution norm.
594: Options Database Keys:
595: . -ksp_cg_radius <r> - Trust Region Radius
597: Notes:
598: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
600: Use preconditioned conjugate gradient to compute
601: an approximate minimizer of the quadratic function
603: q(s) = g^T * s + 0.5 * s^T * H * s
605: subject to the trust region constraint
607: || s || <= delta,
609: where
611: delta is the trust region radius,
612: g is the gradient vector,
613: H is the Hessian approximation, and
614: M is the positive definite preconditioner matrix.
616: `KSPConvergedReason` may be
617: .vb
618: KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
619: KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step,
620: .ve
621: other `KSP` converged/diverged reasons
623: The preconditioner supplied should be symmetric and positive definite.
625: References:
626: + * - Steihaug, T. (1983): The conjugate gradient method and trust regions in large scale optimization. SIAM J. Numer. Anal. 20, 626--637
627: - * - Toint, Ph.L. (1981): Towards an efficient sparsity exploiting Newton method for minimization. In: Duff, I., ed., Sparse Matrices and Their Uses, pp. 57--88. Academic Press
629: Level: developer
631: .seealso: [](ch_ksp), `KSPCreate()`, `KSPCGSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPNASH`, `KSPGLTR`, `KSPQCG`
632: M*/
634: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
635: {
636: KSPCG_STCG *cg;
638: PetscFunctionBegin;
639: PetscCall(PetscNew(&cg));
641: cg->radius = 0.0;
642: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
644: ksp->data = (void *)cg;
645: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
646: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
647: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
648: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
649: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
651: /***************************************************************************/
652: /* Sets the functions that are associated with this data structure */
653: /* (in C++ this is the same as defining virtual functions). */
654: /***************************************************************************/
656: ksp->ops->setup = KSPCGSetUp_STCG;
657: ksp->ops->solve = KSPCGSolve_STCG;
658: ksp->ops->destroy = KSPCGDestroy_STCG;
659: ksp->ops->setfromoptions = KSPCGSetFromOptions_STCG;
660: ksp->ops->buildsolution = KSPBuildSolutionDefault;
661: ksp->ops->buildresidual = KSPBuildResidualDefault;
662: ksp->ops->view = NULL;
664: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_STCG));
665: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_STCG));
666: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_STCG));
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }