Actual source code: stcg.c
1: #include <../src/ksp/ksp/impls/cg/stcg/stcgimpl.h>
3: #define STCG_PRECONDITIONED_DIRECTION 0
4: #define STCG_UNPRECONDITIONED_DIRECTION 1
5: #define STCG_DIRECTION_TYPES 2
7: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
9: static PetscErrorCode KSPCGSolve_STCG(KSP ksp)
10: {
11: #if PetscDefined(USE_COMPLEX)
12: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "STCG is not available for complex systems");
13: #else
14: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
15: Mat Qmat, Mmat;
16: Vec r, z, p, d;
17: PC pc;
18: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
19: PetscReal alpha, beta, kappa, rz, rzm1;
20: PetscReal rr, r2, step;
21: PetscInt max_cg_its;
23: /***************************************************************************/
24: /* Check the arguments and parameters. */
25: /***************************************************************************/
27: PetscFunctionBegin;
28: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
30: /***************************************************************************/
31: /* Get the workspace vectors and initialize variables */
32: /***************************************************************************/
34: r2 = cg->radius * cg->radius;
35: r = ksp->work[0];
36: z = ksp->work[1];
37: p = ksp->work[2];
38: d = ksp->vec_sol;
39: pc = ksp->pc;
41: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
43: PetscCall(VecGetSize(d, &max_cg_its));
44: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
45: ksp->its = 0;
47: /***************************************************************************/
48: /* Initialize objective function and direction. */
49: /***************************************************************************/
51: cg->o_fcn = 0.0;
53: PetscCall(VecSet(d, 0.0)); /* d = 0 */
54: cg->norm_d = 0.0;
56: /***************************************************************************/
57: /* Begin the conjugate gradient method. Check the right-hand side for */
58: /* numerical problems. The check for not-a-number and infinite values */
59: /* need be performed only once. */
60: /***************************************************************************/
62: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
63: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
64: KSPCheckDot(ksp, rr);
66: /***************************************************************************/
67: /* Check the preconditioner for numerical problems and for positive */
68: /* definiteness. The check for not-a-number and infinite values need be */
69: /* performed only once. */
70: /***************************************************************************/
72: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
73: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
74: if (PetscIsInfOrNanScalar(rz)) {
75: /*************************************************************************/
76: /* The preconditioner contains not-a-number or an infinite value. */
77: /* Return the gradient direction intersected with the trust region. */
78: /*************************************************************************/
80: ksp->reason = KSP_DIVERGED_NANORINF;
81: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad preconditioner: rz=%g\n", (double)rz));
83: if (cg->radius != 0) {
84: if (r2 >= rr) {
85: alpha = 1.0;
86: cg->norm_d = PetscSqrtReal(rr);
87: } else {
88: alpha = PetscSqrtReal(r2 / rr);
89: cg->norm_d = cg->radius;
90: }
92: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
94: /***********************************************************************/
95: /* Compute objective function. */
96: /***********************************************************************/
98: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
99: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
100: PetscCall(VecDot(d, z, &cg->o_fcn));
101: cg->o_fcn = -cg->o_fcn;
102: ++ksp->its;
103: }
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: if (rz < 0.0) {
108: /*************************************************************************/
109: /* The preconditioner is indefinite. Because this is the first */
110: /* and we do not have a direction yet, we use the gradient step. Note */
111: /* that we cannot use the preconditioned norm when computing the step */
112: /* because the matrix is indefinite. */
113: /*************************************************************************/
115: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
116: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz));
118: if (cg->radius != 0.0) {
119: if (r2 >= rr) {
120: alpha = 1.0;
121: cg->norm_d = PetscSqrtReal(rr);
122: } else {
123: alpha = PetscSqrtReal(r2 / rr);
124: cg->norm_d = cg->radius;
125: }
127: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
129: /***********************************************************************/
130: /* Compute objective function. */
131: /***********************************************************************/
133: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
134: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
135: PetscCall(VecDot(d, z, &cg->o_fcn));
136: cg->o_fcn = -cg->o_fcn;
137: ++ksp->its;
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /***************************************************************************/
143: /* As far as we know, the preconditioner is positive semidefinite. */
144: /* Compute and log the residual. Check convergence because this */
145: /* initializes things, but do not terminate until at least one conjugate */
146: /* gradient iteration has been performed. */
147: /***************************************************************************/
149: switch (ksp->normtype) {
150: case KSP_NORM_PRECONDITIONED:
151: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
152: break;
154: case KSP_NORM_UNPRECONDITIONED:
155: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
156: break;
158: case KSP_NORM_NATURAL:
159: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
160: break;
162: default:
163: norm_r = 0.0;
164: break;
165: }
167: PetscCall(KSPLogResidualHistory(ksp, norm_r));
168: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
169: ksp->rnorm = norm_r;
171: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
173: /***************************************************************************/
174: /* Compute the first direction and update the iteration. */
175: /***************************************************************************/
177: PetscCall(VecCopy(z, p)); /* p = z */
178: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
179: ++ksp->its;
181: /***************************************************************************/
182: /* Check the matrix for numerical problems. */
183: /***************************************************************************/
185: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
186: if (PetscIsInfOrNanScalar(kappa)) {
187: /*************************************************************************/
188: /* The matrix produced not-a-number or an infinite value. In this case, */
189: /* we must stop and use the gradient direction. This condition need */
190: /* only be checked once. */
191: /*************************************************************************/
193: ksp->reason = KSP_DIVERGED_NANORINF;
194: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad matrix: kappa=%g\n", (double)kappa));
196: if (cg->radius) {
197: if (r2 >= rr) {
198: alpha = 1.0;
199: cg->norm_d = PetscSqrtReal(rr);
200: } else {
201: alpha = PetscSqrtReal(r2 / rr);
202: cg->norm_d = cg->radius;
203: }
205: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
207: /***********************************************************************/
208: /* Compute objective function. */
209: /***********************************************************************/
211: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
212: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
213: PetscCall(VecDot(d, z, &cg->o_fcn));
214: cg->o_fcn = -cg->o_fcn;
215: ++ksp->its;
216: }
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /***************************************************************************/
221: /* Initialize variables for calculating the norm of the direction. */
222: /***************************************************************************/
224: dMp = 0.0;
225: norm_d = 0.0;
226: switch (cg->dtype) {
227: case STCG_PRECONDITIONED_DIRECTION:
228: norm_p = rz;
229: break;
231: default:
232: PetscCall(VecDot(p, p, &norm_p));
233: break;
234: }
236: /***************************************************************************/
237: /* Check for negative curvature. */
238: /***************************************************************************/
240: if (kappa <= 0.0) {
241: /*************************************************************************/
242: /* In this case, the matrix is indefinite and we have encountered a */
243: /* direction of negative curvature. Because negative curvature occurs */
244: /* during the first step, we must follow a direction. */
245: /*************************************************************************/
247: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
248: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));
250: if (cg->radius != 0.0 && norm_p > 0.0) {
251: /***********************************************************************/
252: /* Follow direction of negative curvature to the boundary of the */
253: /* trust region. */
254: /***********************************************************************/
256: step = PetscSqrtReal(r2 / norm_p);
257: cg->norm_d = cg->radius;
259: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
261: /***********************************************************************/
262: /* Update objective function. */
263: /***********************************************************************/
265: cg->o_fcn += step * (0.5 * step * kappa - rz);
266: } else if (cg->radius != 0.0) {
267: /***********************************************************************/
268: /* The norm of the preconditioned direction is zero; use the gradient */
269: /* step. */
270: /***********************************************************************/
272: if (r2 >= rr) {
273: alpha = 1.0;
274: cg->norm_d = PetscSqrtReal(rr);
275: } else {
276: alpha = PetscSqrtReal(r2 / rr);
277: cg->norm_d = cg->radius;
278: }
280: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
282: /***********************************************************************/
283: /* Compute objective function. */
284: /***********************************************************************/
286: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
287: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
288: PetscCall(VecDot(d, z, &cg->o_fcn));
290: cg->o_fcn = -cg->o_fcn;
291: ++ksp->its;
292: }
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /***************************************************************************/
297: /* Run the conjugate gradient method until either the problem is solved, */
298: /* we encounter the boundary of the trust region, or the conjugate */
299: /* gradient method breaks down. */
300: /***************************************************************************/
302: while (1) {
303: /*************************************************************************/
304: /* Know that kappa is nonzero, because we have not broken down, so we */
305: /* can compute the steplength. */
306: /*************************************************************************/
308: alpha = rz / kappa;
310: /*************************************************************************/
311: /* Compute the steplength and check for intersection with the trust */
312: /* region. */
313: /*************************************************************************/
315: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
316: if (cg->radius != 0.0 && norm_dp1 >= r2) {
317: /***********************************************************************/
318: /* In this case, the matrix is positive definite as far as we know. */
319: /* However, the full step goes beyond the trust region. */
320: /***********************************************************************/
322: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
323: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: constrained step: radius=%g\n", (double)cg->radius));
325: if (norm_p > 0.0) {
326: /*********************************************************************/
327: /* Follow the direction to the boundary of the trust region. */
328: /* Final residual norm is never computed. */
329: /*********************************************************************/
331: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
332: cg->norm_d = cg->radius;
334: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
336: /*********************************************************************/
337: /* Update objective function. */
338: /*********************************************************************/
340: cg->o_fcn += step * (0.5 * step * kappa - rz);
341: } else {
342: /*********************************************************************/
343: /* The norm of the direction is zero; there is nothing to follow. */
344: /*********************************************************************/
345: }
346: break;
347: }
349: /*************************************************************************/
350: /* Now we can update the direction and residual. */
351: /*************************************************************************/
353: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
354: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
355: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
357: switch (cg->dtype) {
358: case STCG_PRECONDITIONED_DIRECTION:
359: norm_d = norm_dp1;
360: break;
362: default:
363: PetscCall(VecDot(d, d, &norm_d));
364: break;
365: }
366: cg->norm_d = PetscSqrtReal(norm_d);
368: /*************************************************************************/
369: /* Update objective function. */
370: /*************************************************************************/
372: cg->o_fcn -= 0.5 * alpha * rz;
374: /*************************************************************************/
375: /* Check that the preconditioner appears positive semidefinite. */
376: /*************************************************************************/
378: rzm1 = rz;
379: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
380: if (rz < 0.0) {
381: /***********************************************************************/
382: /* The preconditioner is indefinite. */
383: /***********************************************************************/
385: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
386: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz));
387: break;
388: }
390: /*************************************************************************/
391: /* As far as we know, the preconditioner is positive semidefinite. */
392: /* Compute the residual and check for convergence. */
393: /*************************************************************************/
395: switch (ksp->normtype) {
396: case KSP_NORM_PRECONDITIONED:
397: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
398: break;
400: case KSP_NORM_UNPRECONDITIONED:
401: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
402: break;
404: case KSP_NORM_NATURAL:
405: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
406: break;
408: default:
409: norm_r = 0.0;
410: break;
411: }
413: PetscCall(KSPLogResidualHistory(ksp, norm_r));
414: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
415: ksp->rnorm = norm_r;
417: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
418: if (ksp->reason) {
419: /***********************************************************************/
420: /* The method has converged. */
421: /***********************************************************************/
423: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
424: break;
425: }
427: /*************************************************************************/
428: /* We have not converged yet. Check for breakdown. */
429: /*************************************************************************/
431: beta = rz / rzm1;
432: if (PetscAbsScalar(beta) <= 0.0) {
433: /***********************************************************************/
434: /* Conjugate gradients has broken down. */
435: /***********************************************************************/
437: ksp->reason = KSP_DIVERGED_BREAKDOWN;
438: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: breakdown: beta=%g\n", (double)beta));
439: break;
440: }
442: /*************************************************************************/
443: /* Check iteration limit. */
444: /*************************************************************************/
446: if (ksp->its >= max_cg_its) {
447: ksp->reason = KSP_DIVERGED_ITS;
448: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
449: break;
450: }
452: /*************************************************************************/
453: /* Update p and the norms. */
454: /*************************************************************************/
456: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
458: switch (cg->dtype) {
459: case STCG_PRECONDITIONED_DIRECTION:
460: dMp = beta * (dMp + alpha * norm_p);
461: norm_p = beta * (rzm1 + beta * norm_p);
462: break;
464: default:
465: PetscCall(VecDot(d, p, &dMp));
466: PetscCall(VecDot(p, p, &norm_p));
467: break;
468: }
470: /*************************************************************************/
471: /* Compute the new direction and update the iteration. */
472: /*************************************************************************/
474: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
475: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
476: ++ksp->its;
478: /*************************************************************************/
479: /* Check for negative curvature. */
480: /*************************************************************************/
482: if (kappa <= 0.0) {
483: /***********************************************************************/
484: /* In this case, the matrix is indefinite and we have encountered */
485: /* a direction of negative curvature. Follow the direction to the */
486: /* boundary of the trust region. */
487: /***********************************************************************/
489: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
490: PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));
492: if (cg->radius != 0.0 && norm_p > 0.0) {
493: /*********************************************************************/
494: /* Follow direction of negative curvature to boundary. */
495: /* Final residual norm is never computed. */
496: /*********************************************************************/
498: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
499: cg->norm_d = cg->radius;
501: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
503: /*********************************************************************/
504: /* Update objective function. */
505: /*********************************************************************/
507: cg->o_fcn += step * (0.5 * step * kappa - rz);
508: } else if (cg->radius != 0.0) {
509: /*********************************************************************/
510: /* The norm of the direction is zero; there is nothing to follow. */
511: /*********************************************************************/
512: }
513: break;
514: }
515: }
516: PetscFunctionReturn(PETSC_SUCCESS);
517: #endif
518: }
520: static PetscErrorCode KSPCGSetUp_STCG(KSP ksp)
521: {
522: PetscFunctionBegin;
523: /***************************************************************************/
524: /* Set work vectors needed by conjugate gradient method and allocate */
525: /***************************************************************************/
527: PetscCall(KSPSetWorkVecs(ksp, 3));
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: static PetscErrorCode KSPCGDestroy_STCG(KSP ksp)
532: {
533: PetscFunctionBegin;
534: /***************************************************************************/
535: /* Clear composed functions */
536: /***************************************************************************/
538: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
539: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
540: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
542: /***************************************************************************/
543: /* Destroy KSP object. */
544: /***************************************************************************/
546: PetscCall(KSPDestroyDefault(ksp));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: static PetscErrorCode KSPCGSetRadius_STCG(KSP ksp, PetscReal radius)
551: {
552: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
554: PetscFunctionBegin;
555: cg->radius = radius;
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: static PetscErrorCode KSPCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
560: {
561: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
563: PetscFunctionBegin;
564: *norm_d = cg->norm_d;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode KSPCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
569: {
570: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
572: PetscFunctionBegin;
573: *o_fcn = cg->o_fcn;
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode KSPCGSetFromOptions_STCG(KSP ksp, PetscOptionItems PetscOptionsObject)
578: {
579: KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
581: PetscFunctionBegin;
582: PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG STCG options");
583: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
584: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
585: PetscOptionsHeadEnd();
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*MC
590: KSPSTCG - Code to run conjugate gradient method subject to a constraint on the solution norm for use in a trust region method
591: {cite}`steihaug:83`, {cite}`toint1981towards`
593: Options Database Key:
594: . -ksp_cg_radius radius - Trust Region Radius
596: Level: developer
598: Notes:
599: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
601: Use preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
603: $$
604: q(s) = g^T s + \frac{1}{2} s^T H s
605: $$
607: subject to the trust region constraint
609: $$
610: || s || \le \delta,
611: $$
613: where
614: .vb
615: delta is the trust region radius,
616: g is the gradient vector,
617: H is the Hessian approximation, and
618: .ve
620: `KSPConvergedReason` may include
621: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
622: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,
624: The preconditioner supplied should be symmetric and positive definite.
626: .seealso: [](ch_ksp), `KSPCreate()`, `KSPCGSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPNASH`, `KSPGLTR`, `KSPQCG`
627: M*/
629: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
630: {
631: KSPCG_STCG *cg;
633: PetscFunctionBegin;
634: PetscCall(PetscNew(&cg));
636: cg->radius = 0.0;
637: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
639: ksp->data = (void *)cg;
640: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
641: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
642: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
643: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
644: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
646: /***************************************************************************/
647: /* Sets the functions that are associated with this data structure */
648: /* (in C++ this is the same as defining virtual functions). */
649: /***************************************************************************/
651: ksp->ops->setup = KSPCGSetUp_STCG;
652: ksp->ops->solve = KSPCGSolve_STCG;
653: ksp->ops->destroy = KSPCGDestroy_STCG;
654: ksp->ops->setfromoptions = KSPCGSetFromOptions_STCG;
655: ksp->ops->buildsolution = KSPBuildSolutionDefault;
656: ksp->ops->buildresidual = KSPBuildResidualDefault;
657: ksp->ops->view = NULL;
659: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_STCG));
660: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_STCG));
661: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_STCG));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }