Actual source code: stcg.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/stcg/stcgimpl.h>
5: #define STCG_PRECONDITIONED_DIRECTION 0
6: #define STCG_UNPRECONDITIONED_DIRECTION 1
7: #define STCG_DIRECTION_TYPES 2
9: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
13: /*@
14: KSPSTCGSetRadius - 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_stcg_radius <r>
25: Level: advanced
27: .keywords: KSP, STCG, set, trust region radius
28: @*/
29: PetscErrorCode KSPSTCGSetRadius(KSP ksp, PetscReal radius)
30: {
35: if (radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Radius negative");
37: PetscTryMethod(ksp,"KSPSTCGSetRadius_C",(KSP,PetscReal),(ksp,radius));
38: return(0);
39: }
43: /*@
44: KSPSTCGGetNormD - 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, STCG, get, norm direction
55: @*/
56: PetscErrorCode KSPSTCGGetNormD(KSP ksp, PetscReal *norm_d)
57: {
62: PetscUseMethod(ksp,"KSPSTCGGetNormD_C",(KSP,PetscReal*),(ksp,norm_d));
63: return(0);
64: }
68: /*@
69: KSPSTCGGetObjFcn - 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, STCG, get, objective function
80: @*/
81: PetscErrorCode KSPSTCGGetObjFcn(KSP ksp, PetscReal *o_fcn)
82: {
87: PetscUseMethod(ksp,"KSPSTCGGetObjFcn_C",(KSP,PetscReal*),(ksp,o_fcn));
88: return(0);
89: }
93: PetscErrorCode KSPSolve_STCG(KSP ksp)
94: {
95: #if defined(PETSC_USE_COMPLEX)
96: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "STCG is not available for complex systems");
97: #else
98: KSP_STCG *cg = (KSP_STCG*)ksp->data;
100: Mat Qmat, Mmat;
101: Vec r, z, p, d;
102: PC pc;
103: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
104: PetscReal alpha, beta, kappa, rz, rzm1;
105: PetscReal rr, r2, step;
106: PetscInt max_cg_its;
107: PetscBool diagonalscale;
109: /***************************************************************************/
110: /* Check the arguments and parameters. */
111: /***************************************************************************/
114: PCGetDiagonalScale(ksp->pc, &diagonalscale);
115: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
116: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
118: /***************************************************************************/
119: /* Get the workspace vectors and initialize variables */
120: /***************************************************************************/
122: r2 = cg->radius * cg->radius;
123: r = ksp->work[0];
124: z = ksp->work[1];
125: p = ksp->work[2];
126: d = ksp->vec_sol;
127: pc = ksp->pc;
129: PCGetOperators(pc, &Qmat, &Mmat);
131: VecGetSize(d, &max_cg_its);
132: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
133: ksp->its = 0;
135: /***************************************************************************/
136: /* Initialize objective function and direction. */
137: /***************************************************************************/
139: cg->o_fcn = 0.0;
141: VecSet(d, 0.0); /* d = 0 */
142: cg->norm_d = 0.0;
144: /***************************************************************************/
145: /* Begin the conjugate gradient method. Check the right-hand side for */
146: /* numerical problems. The check for not-a-number and infinite values */
147: /* need be performed only once. */
148: /***************************************************************************/
150: VecCopy(ksp->vec_rhs, r); /* r = -grad */
151: VecDot(r, r, &rr); /* rr = r^T r */
152: KSPCheckDot(ksp,rr);
154: /***************************************************************************/
155: /* Check the preconditioner for numerical problems and for positive */
156: /* definiteness. The check for not-a-number and infinite values need be */
157: /* performed only once. */
158: /***************************************************************************/
160: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
161: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
162: if (PetscIsInfOrNanScalar(rz)) {
163: /*************************************************************************/
164: /* The preconditioner contains not-a-number or an infinite value. */
165: /* Return the gradient direction intersected with the trust region. */
166: /*************************************************************************/
168: ksp->reason = KSP_DIVERGED_NANORINF;
169: PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", (double)rz);
171: if (cg->radius != 0) {
172: if (r2 >= rr) {
173: alpha = 1.0;
174: cg->norm_d = PetscSqrtReal(rr);
175: } else {
176: alpha = PetscSqrtReal(r2 / rr);
177: cg->norm_d = cg->radius;
178: }
180: VecAXPY(d, alpha, r); /* d = d + alpha r */
182: /***********************************************************************/
183: /* Compute objective function. */
184: /***********************************************************************/
186: KSP_MatMult(ksp, Qmat, d, z);
187: VecAYPX(z, -0.5, ksp->vec_rhs);
188: VecDot(d, z, &cg->o_fcn);
189: cg->o_fcn = -cg->o_fcn;
190: ++ksp->its;
191: }
192: return(0);
193: }
195: if (rz < 0.0) {
196: /*************************************************************************/
197: /* The preconditioner is indefinite. Because this is the first */
198: /* and we do not have a direction yet, we use the gradient step. Note */
199: /* that we cannot use the preconditioned norm when computing the step */
200: /* because the matrix is indefinite. */
201: /*************************************************************************/
203: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
204: PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz);
206: if (cg->radius != 0.0) {
207: if (r2 >= rr) {
208: alpha = 1.0;
209: cg->norm_d = PetscSqrtReal(rr);
210: } else {
211: alpha = PetscSqrtReal(r2 / rr);
212: cg->norm_d = cg->radius;
213: }
215: VecAXPY(d, alpha, r); /* d = d + alpha r */
217: /***********************************************************************/
218: /* Compute objective function. */
219: /***********************************************************************/
221: KSP_MatMult(ksp, Qmat, d, z);
222: VecAYPX(z, -0.5, ksp->vec_rhs);
223: VecDot(d, z, &cg->o_fcn);
224: cg->o_fcn = -cg->o_fcn;
225: ++ksp->its;
226: }
227: return(0);
228: }
230: /***************************************************************************/
231: /* As far as we know, the preconditioner is positive semidefinite. */
232: /* Compute and log the residual. Check convergence because this */
233: /* initializes things, but do not terminate until at least one conjugate */
234: /* gradient iteration has been performed. */
235: /***************************************************************************/
237: switch (ksp->normtype) {
238: case KSP_NORM_PRECONDITIONED:
239: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
240: break;
242: case KSP_NORM_UNPRECONDITIONED:
243: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
244: break;
246: case KSP_NORM_NATURAL:
247: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
248: break;
250: default:
251: norm_r = 0.0;
252: break;
253: }
255: KSPLogResidualHistory(ksp, norm_r);
256: KSPMonitor(ksp, ksp->its, norm_r);
257: ksp->rnorm = norm_r;
259: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
261: /***************************************************************************/
262: /* Compute the first direction and update the iteration. */
263: /***************************************************************************/
265: VecCopy(z, p); /* p = z */
266: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
267: ++ksp->its;
269: /***************************************************************************/
270: /* Check the matrix for numerical problems. */
271: /***************************************************************************/
273: VecDot(p, z, &kappa); /* kappa = p^T Q p */
274: if (PetscIsInfOrNanScalar(kappa)) {
275: /*************************************************************************/
276: /* The matrix produced not-a-number or an infinite value. In this case, */
277: /* we must stop and use the gradient direction. This condition need */
278: /* only be checked once. */
279: /*************************************************************************/
281: ksp->reason = KSP_DIVERGED_NANORINF;
282: PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", (double)kappa);
284: if (cg->radius) {
285: if (r2 >= rr) {
286: alpha = 1.0;
287: cg->norm_d = PetscSqrtReal(rr);
288: } else {
289: alpha = PetscSqrtReal(r2 / rr);
290: cg->norm_d = cg->radius;
291: }
293: VecAXPY(d, alpha, r); /* d = d + alpha r */
295: /***********************************************************************/
296: /* Compute objective function. */
297: /***********************************************************************/
299: KSP_MatMult(ksp, Qmat, d, z);
300: VecAYPX(z, -0.5, ksp->vec_rhs);
301: VecDot(d, z, &cg->o_fcn);
302: cg->o_fcn = -cg->o_fcn;
303: ++ksp->its;
304: }
305: return(0);
306: }
308: /***************************************************************************/
309: /* Initialize variables for calculating the norm of the direction. */
310: /***************************************************************************/
312: dMp = 0.0;
313: norm_d = 0.0;
314: switch (cg->dtype) {
315: case STCG_PRECONDITIONED_DIRECTION:
316: norm_p = rz;
317: break;
319: default:
320: VecDot(p, p, &norm_p);
321: break;
322: }
324: /***************************************************************************/
325: /* Check for negative curvature. */
326: /***************************************************************************/
328: if (kappa <= 0.0) {
329: /*************************************************************************/
330: /* In this case, the matrix is indefinite and we have encountered a */
331: /* direction of negative curvature. Because negative curvature occurs */
332: /* during the first step, we must follow a direction. */
333: /*************************************************************************/
335: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
336: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);
338: if (cg->radius != 0.0 && norm_p > 0.0) {
339: /***********************************************************************/
340: /* Follow direction of negative curvature to the boundary of the */
341: /* trust region. */
342: /***********************************************************************/
344: step = PetscSqrtReal(r2 / norm_p);
345: cg->norm_d = cg->radius;
347: VecAXPY(d, step, p); /* d = d + step p */
349: /***********************************************************************/
350: /* Update objective function. */
351: /***********************************************************************/
353: cg->o_fcn += step * (0.5 * step * kappa - rz);
354: } else if (cg->radius != 0.0) {
355: /***********************************************************************/
356: /* The norm of the preconditioned direction is zero; use the gradient */
357: /* step. */
358: /***********************************************************************/
360: if (r2 >= rr) {
361: alpha = 1.0;
362: cg->norm_d = PetscSqrtReal(rr);
363: } else {
364: alpha = PetscSqrtReal(r2 / rr);
365: cg->norm_d = cg->radius;
366: }
368: VecAXPY(d, alpha, r); /* d = d + alpha r */
370: /***********************************************************************/
371: /* Compute objective function. */
372: /***********************************************************************/
374: KSP_MatMult(ksp, Qmat, d, z);
375: VecAYPX(z, -0.5, ksp->vec_rhs);
376: VecDot(d, z, &cg->o_fcn);
378: cg->o_fcn = -cg->o_fcn;
379: ++ksp->its;
380: }
381: return(0);
382: }
384: /***************************************************************************/
385: /* Run the conjugate gradient method until either the problem is solved, */
386: /* we encounter the boundary of the trust region, or the conjugate */
387: /* gradient method breaks down. */
388: /***************************************************************************/
390: while (1) {
391: /*************************************************************************/
392: /* Know that kappa is nonzero, because we have not broken down, so we */
393: /* can compute the steplength. */
394: /*************************************************************************/
396: alpha = rz / kappa;
398: /*************************************************************************/
399: /* Compute the steplength and check for intersection with the trust */
400: /* region. */
401: /*************************************************************************/
403: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
404: if (cg->radius != 0.0 && norm_dp1 >= r2) {
405: /***********************************************************************/
406: /* In this case, the matrix is positive definite as far as we know. */
407: /* However, the full step goes beyond the trust region. */
408: /***********************************************************************/
410: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
411: PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", (double)cg->radius);
413: if (norm_p > 0.0) {
414: /*********************************************************************/
415: /* Follow the direction to the boundary of the trust region. */
416: /*********************************************************************/
418: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
419: cg->norm_d = cg->radius;
421: VecAXPY(d, step, p); /* d = d + step p */
423: /*********************************************************************/
424: /* Update objective function. */
425: /*********************************************************************/
427: cg->o_fcn += step * (0.5 * step * kappa - rz);
428: } else {
429: /*********************************************************************/
430: /* The norm of the direction is zero; there is nothing to follow. */
431: /*********************************************************************/
432: }
433: break;
434: }
436: /*************************************************************************/
437: /* Now we can update the direction and residual. */
438: /*************************************************************************/
440: VecAXPY(d, alpha, p); /* d = d + alpha p */
441: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
442: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
444: switch (cg->dtype) {
445: case STCG_PRECONDITIONED_DIRECTION:
446: norm_d = norm_dp1;
447: break;
449: default:
450: VecDot(d, d, &norm_d);
451: break;
452: }
453: cg->norm_d = PetscSqrtReal(norm_d);
455: /*************************************************************************/
456: /* Update objective function. */
457: /*************************************************************************/
459: cg->o_fcn -= 0.5 * alpha * rz;
461: /*************************************************************************/
462: /* Check that the preconditioner appears positive semidefinite. */
463: /*************************************************************************/
465: rzm1 = rz;
466: VecDot(r, z, &rz); /* rz = r^T z */
467: if (rz < 0.0) {
468: /***********************************************************************/
469: /* The preconditioner is indefinite. */
470: /***********************************************************************/
472: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
473: PetscInfo1(ksp, "KSPSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz);
474: break;
475: }
477: /*************************************************************************/
478: /* As far as we know, the preconditioner is positive semidefinite. */
479: /* Compute the residual and check for convergence. */
480: /*************************************************************************/
482: switch (ksp->normtype) {
483: case KSP_NORM_PRECONDITIONED:
484: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
485: break;
487: case KSP_NORM_UNPRECONDITIONED:
488: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
489: break;
491: case KSP_NORM_NATURAL:
492: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
493: break;
495: default:
496: norm_r = 0.0;
497: break;
498: }
500: KSPLogResidualHistory(ksp, norm_r);
501: KSPMonitor(ksp, ksp->its, norm_r);
502: ksp->rnorm = norm_r;
504: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
505: if (ksp->reason) {
506: /***********************************************************************/
507: /* The method has converged. */
508: /***********************************************************************/
510: PetscInfo2(ksp, "KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
511: break;
512: }
514: /*************************************************************************/
515: /* We have not converged yet. Check for breakdown. */
516: /*************************************************************************/
518: beta = rz / rzm1;
519: if (PetscAbsScalar(beta) <= 0.0) {
520: /***********************************************************************/
521: /* Conjugate gradients has broken down. */
522: /***********************************************************************/
524: ksp->reason = KSP_DIVERGED_BREAKDOWN;
525: PetscInfo1(ksp, "KSPSolve_STCG: breakdown: beta=%g\n", (double)beta);
526: break;
527: }
529: /*************************************************************************/
530: /* Check iteration limit. */
531: /*************************************************************************/
533: if (ksp->its >= max_cg_its) {
534: ksp->reason = KSP_DIVERGED_ITS;
535: PetscInfo1(ksp, "KSPSolve_STCG: iterlim: its=%D\n", ksp->its);
536: break;
537: }
539: /*************************************************************************/
540: /* Update p and the norms. */
541: /*************************************************************************/
543: VecAYPX(p, beta, z); /* p = z + beta p */
545: switch (cg->dtype) {
546: case STCG_PRECONDITIONED_DIRECTION:
547: dMp = beta*(dMp + alpha*norm_p);
548: norm_p = beta*(rzm1 + beta*norm_p);
549: break;
551: default:
552: VecDot(d, p, &dMp);
553: VecDot(p, p, &norm_p);
554: break;
555: }
557: /*************************************************************************/
558: /* Compute the new direction and update the iteration. */
559: /*************************************************************************/
561: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
562: VecDot(p, z, &kappa); /* kappa = p^T Q p */
563: ++ksp->its;
565: /*************************************************************************/
566: /* Check for negative curvature. */
567: /*************************************************************************/
569: if (kappa <= 0.0) {
570: /***********************************************************************/
571: /* In this case, the matrix is indefinite and we have encountered */
572: /* a direction of negative curvature. Follow the direction to the */
573: /* boundary of the trust region. */
574: /***********************************************************************/
576: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
577: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);
579: if (cg->radius != 0.0 && norm_p > 0.0) {
580: /*********************************************************************/
581: /* Follow direction of negative curvature to boundary. */
582: /*********************************************************************/
584: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
585: cg->norm_d = cg->radius;
587: VecAXPY(d, step, p); /* d = d + step p */
589: /*********************************************************************/
590: /* Update objective function. */
591: /*********************************************************************/
593: cg->o_fcn += step * (0.5 * step * kappa - rz);
594: } else if (cg->radius != 0.0) {
595: /*********************************************************************/
596: /* The norm of the direction is zero; there is nothing to follow. */
597: /*********************************************************************/
598: }
599: break;
600: }
601: }
602: return(0);
603: #endif
604: }
608: PetscErrorCode KSPSetUp_STCG(KSP ksp)
609: {
612: /***************************************************************************/
613: /* Set work vectors needed by conjugate gradient method and allocate */
614: /***************************************************************************/
617: KSPSetWorkVecs(ksp, 3);
618: return(0);
619: }
623: PetscErrorCode KSPDestroy_STCG(KSP ksp)
624: {
627: /***************************************************************************/
628: /* Clear composed functions */
629: /***************************************************************************/
632: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",NULL);
633: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",NULL);
634: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",NULL);
636: /***************************************************************************/
637: /* Destroy KSP object. */
638: /***************************************************************************/
640: KSPDestroyDefault(ksp);
641: return(0);
642: }
646: static PetscErrorCode KSPSTCGSetRadius_STCG(KSP ksp, PetscReal radius)
647: {
648: KSP_STCG *cg = (KSP_STCG*)ksp->data;
651: cg->radius = radius;
652: return(0);
653: }
657: static PetscErrorCode KSPSTCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
658: {
659: KSP_STCG *cg = (KSP_STCG*)ksp->data;
662: *norm_d = cg->norm_d;
663: return(0);
664: }
668: static PetscErrorCode KSPSTCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
669: {
670: KSP_STCG *cg = (KSP_STCG*)ksp->data;
673: *o_fcn = cg->o_fcn;
674: return(0);
675: }
679: PetscErrorCode KSPSetFromOptions_STCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
680: {
682: KSP_STCG *cg = (KSP_STCG*)ksp->data;
685: PetscOptionsHead(PetscOptionsObject,"KSP STCG options");
686: PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, NULL);
687: PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
688: PetscOptionsTail();
689: return(0);
690: }
692: /*MC
693: KSPSTCG - Code to run conjugate gradient method subject to a constraint
694: on the solution norm. This is used in Trust Region methods for
695: nonlinear equations, SNESNEWTONTR
697: Options Database Keys:
698: . -ksp_stcg_radius <r> - Trust Region Radius
700: Notes: This is rarely used directly
702: Use preconditioned conjugate gradient to compute
703: an approximate minimizer of the quadratic function
705: q(s) = g^T * s + 0.5 * s^T * H * s
707: subject to the trust region constraint
709: || s || <= delta,
711: where
713: delta is the trust region radius,
714: g is the gradient vector,
715: H is the Hessian approximation, and
716: M is the positive definite preconditioner matrix.
718: KSPConvergedReason may be
719: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
720: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
721: $ other KSP converged/diverged reasons
723: Notes:
724: The preconditioner supplied should be symmetric and positive definite.
726: Level: developer
728: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD(), KSPSTCGGetObjFcn()
729: M*/
733: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
734: {
736: KSP_STCG *cg;
739: PetscNewLog(ksp,&cg);
741: cg->radius = 0.0;
742: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
744: ksp->data = (void*) cg;
745: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
746: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
747: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
749: /***************************************************************************/
750: /* Sets the functions that are associated with this data structure */
751: /* (in C++ this is the same as defining virtual functions). */
752: /***************************************************************************/
754: ksp->ops->setup = KSPSetUp_STCG;
755: ksp->ops->solve = KSPSolve_STCG;
756: ksp->ops->destroy = KSPDestroy_STCG;
757: ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
758: ksp->ops->buildsolution = KSPBuildSolutionDefault;
759: ksp->ops->buildresidual = KSPBuildResidualDefault;
760: ksp->ops->view = 0;
762: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",KSPSTCGSetRadius_STCG);
763: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",KSPSTCGGetNormD_STCG);
764: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",KSPSTCGGetObjFcn_STCG);
765: return(0);
766: }