Actual source code: stcg.c
petsc-3.5.4 2015-05-23
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: if (PetscIsInfOrNanScalar(rr)) {
153: /*************************************************************************/
154: /* The right-hand side contains not-a-number or an infinite value. */
155: /* The gradient step does not work; return a zero value for the step. */
156: /*************************************************************************/
158: ksp->reason = KSP_DIVERGED_NANORINF;
159: PetscInfo1(ksp, "KSPSolve_STCG: bad right-hand side: rr=%g\n", (double)rr);
160: return(0);
161: }
163: /***************************************************************************/
164: /* Check the preconditioner for numerical problems and for positive */
165: /* definiteness. The check for not-a-number and infinite values need be */
166: /* performed only once. */
167: /***************************************************************************/
169: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
170: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
171: if (PetscIsInfOrNanScalar(rz)) {
172: /*************************************************************************/
173: /* The preconditioner contains not-a-number or an infinite value. */
174: /* Return the gradient direction intersected with the trust region. */
175: /*************************************************************************/
177: ksp->reason = KSP_DIVERGED_NANORINF;
178: PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", (double)rz);
180: if (cg->radius != 0) {
181: if (r2 >= rr) {
182: alpha = 1.0;
183: cg->norm_d = PetscSqrtReal(rr);
184: } else {
185: alpha = PetscSqrtReal(r2 / rr);
186: cg->norm_d = cg->radius;
187: }
189: VecAXPY(d, alpha, r); /* d = d + alpha r */
191: /***********************************************************************/
192: /* Compute objective function. */
193: /***********************************************************************/
195: KSP_MatMult(ksp, Qmat, d, z);
196: VecAYPX(z, -0.5, ksp->vec_rhs);
197: VecDot(d, z, &cg->o_fcn);
198: cg->o_fcn = -cg->o_fcn;
199: ++ksp->its;
200: }
201: return(0);
202: }
204: if (rz < 0.0) {
205: /*************************************************************************/
206: /* The preconditioner is indefinite. Because this is the first */
207: /* and we do not have a direction yet, we use the gradient step. Note */
208: /* that we cannot use the preconditioned norm when computing the step */
209: /* because the matrix is indefinite. */
210: /*************************************************************************/
212: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
213: PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz);
215: if (cg->radius != 0.0) {
216: if (r2 >= rr) {
217: alpha = 1.0;
218: cg->norm_d = PetscSqrtReal(rr);
219: } else {
220: alpha = PetscSqrtReal(r2 / rr);
221: cg->norm_d = cg->radius;
222: }
224: VecAXPY(d, alpha, r); /* d = d + alpha r */
226: /***********************************************************************/
227: /* Compute objective function. */
228: /***********************************************************************/
230: KSP_MatMult(ksp, Qmat, d, z);
231: VecAYPX(z, -0.5, ksp->vec_rhs);
232: VecDot(d, z, &cg->o_fcn);
233: cg->o_fcn = -cg->o_fcn;
234: ++ksp->its;
235: }
236: return(0);
237: }
239: /***************************************************************************/
240: /* As far as we know, the preconditioner is positive semidefinite. */
241: /* Compute and log the residual. Check convergence because this */
242: /* initializes things, but do not terminate until at least one conjugate */
243: /* gradient iteration has been performed. */
244: /***************************************************************************/
246: switch (ksp->normtype) {
247: case KSP_NORM_PRECONDITIONED:
248: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
249: break;
251: case KSP_NORM_UNPRECONDITIONED:
252: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
253: break;
255: case KSP_NORM_NATURAL:
256: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
257: break;
259: default:
260: norm_r = 0.0;
261: break;
262: }
264: KSPLogResidualHistory(ksp, norm_r);
265: KSPMonitor(ksp, ksp->its, norm_r);
266: ksp->rnorm = norm_r;
268: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
270: /***************************************************************************/
271: /* Compute the first direction and update the iteration. */
272: /***************************************************************************/
274: VecCopy(z, p); /* p = z */
275: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
276: ++ksp->its;
278: /***************************************************************************/
279: /* Check the matrix for numerical problems. */
280: /***************************************************************************/
282: VecDot(p, z, &kappa); /* kappa = p^T Q p */
283: if (PetscIsInfOrNanScalar(kappa)) {
284: /*************************************************************************/
285: /* The matrix produced not-a-number or an infinite value. In this case, */
286: /* we must stop and use the gradient direction. This condition need */
287: /* only be checked once. */
288: /*************************************************************************/
290: ksp->reason = KSP_DIVERGED_NANORINF;
291: PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", (double)kappa);
293: if (cg->radius) {
294: if (r2 >= rr) {
295: alpha = 1.0;
296: cg->norm_d = PetscSqrtReal(rr);
297: } else {
298: alpha = PetscSqrtReal(r2 / rr);
299: cg->norm_d = cg->radius;
300: }
302: VecAXPY(d, alpha, r); /* d = d + alpha r */
304: /***********************************************************************/
305: /* Compute objective function. */
306: /***********************************************************************/
308: KSP_MatMult(ksp, Qmat, d, z);
309: VecAYPX(z, -0.5, ksp->vec_rhs);
310: VecDot(d, z, &cg->o_fcn);
311: cg->o_fcn = -cg->o_fcn;
312: ++ksp->its;
313: }
314: return(0);
315: }
317: /***************************************************************************/
318: /* Initialize variables for calculating the norm of the direction. */
319: /***************************************************************************/
321: dMp = 0.0;
322: norm_d = 0.0;
323: switch (cg->dtype) {
324: case STCG_PRECONDITIONED_DIRECTION:
325: norm_p = rz;
326: break;
328: default:
329: VecDot(p, p, &norm_p);
330: break;
331: }
333: /***************************************************************************/
334: /* Check for negative curvature. */
335: /***************************************************************************/
337: if (kappa <= 0.0) {
338: /*************************************************************************/
339: /* In this case, the matrix is indefinite and we have encountered a */
340: /* direction of negative curvature. Because negative curvature occurs */
341: /* during the first step, we must follow a direction. */
342: /*************************************************************************/
344: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
345: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);
347: if (cg->radius != 0.0 && norm_p > 0.0) {
348: /***********************************************************************/
349: /* Follow direction of negative curvature to the boundary of the */
350: /* trust region. */
351: /***********************************************************************/
353: step = PetscSqrtReal(r2 / norm_p);
354: cg->norm_d = cg->radius;
356: VecAXPY(d, step, p); /* d = d + step p */
358: /***********************************************************************/
359: /* Update objective function. */
360: /***********************************************************************/
362: cg->o_fcn += step * (0.5 * step * kappa - rz);
363: } else if (cg->radius != 0.0) {
364: /***********************************************************************/
365: /* The norm of the preconditioned direction is zero; use the gradient */
366: /* step. */
367: /***********************************************************************/
369: if (r2 >= rr) {
370: alpha = 1.0;
371: cg->norm_d = PetscSqrtReal(rr);
372: } else {
373: alpha = PetscSqrtReal(r2 / rr);
374: cg->norm_d = cg->radius;
375: }
377: VecAXPY(d, alpha, r); /* d = d + alpha r */
379: /***********************************************************************/
380: /* Compute objective function. */
381: /***********************************************************************/
383: KSP_MatMult(ksp, Qmat, d, z);
384: VecAYPX(z, -0.5, ksp->vec_rhs);
385: VecDot(d, z, &cg->o_fcn);
387: cg->o_fcn = -cg->o_fcn;
388: ++ksp->its;
389: }
390: return(0);
391: }
393: /***************************************************************************/
394: /* Run the conjugate gradient method until either the problem is solved, */
395: /* we encounter the boundary of the trust region, or the conjugate */
396: /* gradient method breaks down. */
397: /***************************************************************************/
399: while (1) {
400: /*************************************************************************/
401: /* Know that kappa is nonzero, because we have not broken down, so we */
402: /* can compute the steplength. */
403: /*************************************************************************/
405: alpha = rz / kappa;
407: /*************************************************************************/
408: /* Compute the steplength and check for intersection with the trust */
409: /* region. */
410: /*************************************************************************/
412: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
413: if (cg->radius != 0.0 && norm_dp1 >= r2) {
414: /***********************************************************************/
415: /* In this case, the matrix is positive definite as far as we know. */
416: /* However, the full step goes beyond the trust region. */
417: /***********************************************************************/
419: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
420: PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", (double)cg->radius);
422: if (norm_p > 0.0) {
423: /*********************************************************************/
424: /* Follow the direction to the boundary of the trust region. */
425: /*********************************************************************/
427: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
428: cg->norm_d = cg->radius;
430: VecAXPY(d, step, p); /* d = d + step p */
432: /*********************************************************************/
433: /* Update objective function. */
434: /*********************************************************************/
436: cg->o_fcn += step * (0.5 * step * kappa - rz);
437: } else {
438: /*********************************************************************/
439: /* The norm of the direction is zero; there is nothing to follow. */
440: /*********************************************************************/
441: }
442: break;
443: }
445: /*************************************************************************/
446: /* Now we can update the direction and residual. */
447: /*************************************************************************/
449: VecAXPY(d, alpha, p); /* d = d + alpha p */
450: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
451: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
453: switch (cg->dtype) {
454: case STCG_PRECONDITIONED_DIRECTION:
455: norm_d = norm_dp1;
456: break;
458: default:
459: VecDot(d, d, &norm_d);
460: break;
461: }
462: cg->norm_d = PetscSqrtReal(norm_d);
464: /*************************************************************************/
465: /* Update objective function. */
466: /*************************************************************************/
468: cg->o_fcn -= 0.5 * alpha * rz;
470: /*************************************************************************/
471: /* Check that the preconditioner appears positive semidefinite. */
472: /*************************************************************************/
474: rzm1 = rz;
475: VecDot(r, z, &rz); /* rz = r^T z */
476: if (rz < 0.0) {
477: /***********************************************************************/
478: /* The preconditioner is indefinite. */
479: /***********************************************************************/
481: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
482: PetscInfo1(ksp, "KSPSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz);
483: break;
484: }
486: /*************************************************************************/
487: /* As far as we know, the preconditioner is positive semidefinite. */
488: /* Compute the residual and check for convergence. */
489: /*************************************************************************/
491: switch (ksp->normtype) {
492: case KSP_NORM_PRECONDITIONED:
493: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
494: break;
496: case KSP_NORM_UNPRECONDITIONED:
497: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
498: break;
500: case KSP_NORM_NATURAL:
501: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
502: break;
504: default:
505: norm_r = 0.0;
506: break;
507: }
509: KSPLogResidualHistory(ksp, norm_r);
510: KSPMonitor(ksp, ksp->its, norm_r);
511: ksp->rnorm = norm_r;
513: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
514: if (ksp->reason) {
515: /***********************************************************************/
516: /* The method has converged. */
517: /***********************************************************************/
519: PetscInfo2(ksp, "KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
520: break;
521: }
523: /*************************************************************************/
524: /* We have not converged yet. Check for breakdown. */
525: /*************************************************************************/
527: beta = rz / rzm1;
528: if (PetscAbsScalar(beta) <= 0.0) {
529: /***********************************************************************/
530: /* Conjugate gradients has broken down. */
531: /***********************************************************************/
533: ksp->reason = KSP_DIVERGED_BREAKDOWN;
534: PetscInfo1(ksp, "KSPSolve_STCG: breakdown: beta=%g\n", (double)beta);
535: break;
536: }
538: /*************************************************************************/
539: /* Check iteration limit. */
540: /*************************************************************************/
542: if (ksp->its >= max_cg_its) {
543: ksp->reason = KSP_DIVERGED_ITS;
544: PetscInfo1(ksp, "KSPSolve_STCG: iterlim: its=%D\n", ksp->its);
545: break;
546: }
548: /*************************************************************************/
549: /* Update p and the norms. */
550: /*************************************************************************/
552: VecAYPX(p, beta, z); /* p = z + beta p */
554: switch (cg->dtype) {
555: case STCG_PRECONDITIONED_DIRECTION:
556: dMp = beta*(dMp + alpha*norm_p);
557: norm_p = beta*(rzm1 + beta*norm_p);
558: break;
560: default:
561: VecDot(d, p, &dMp);
562: VecDot(p, p, &norm_p);
563: break;
564: }
566: /*************************************************************************/
567: /* Compute the new direction and update the iteration. */
568: /*************************************************************************/
570: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
571: VecDot(p, z, &kappa); /* kappa = p^T Q p */
572: ++ksp->its;
574: /*************************************************************************/
575: /* Check for negative curvature. */
576: /*************************************************************************/
578: if (kappa <= 0.0) {
579: /***********************************************************************/
580: /* In this case, the matrix is indefinite and we have encountered */
581: /* a direction of negative curvature. Follow the direction to the */
582: /* boundary of the trust region. */
583: /***********************************************************************/
585: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
586: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);
588: if (cg->radius != 0.0 && norm_p > 0.0) {
589: /*********************************************************************/
590: /* Follow direction of negative curvature to boundary. */
591: /*********************************************************************/
593: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
594: cg->norm_d = cg->radius;
596: VecAXPY(d, step, p); /* d = d + step p */
598: /*********************************************************************/
599: /* Update objective function. */
600: /*********************************************************************/
602: cg->o_fcn += step * (0.5 * step * kappa - rz);
603: } else if (cg->radius != 0.0) {
604: /*********************************************************************/
605: /* The norm of the direction is zero; there is nothing to follow. */
606: /*********************************************************************/
607: }
608: break;
609: }
610: }
611: return(0);
612: #endif
613: }
617: PetscErrorCode KSPSetUp_STCG(KSP ksp)
618: {
621: /***************************************************************************/
622: /* Set work vectors needed by conjugate gradient method and allocate */
623: /***************************************************************************/
626: KSPSetWorkVecs(ksp, 3);
627: return(0);
628: }
632: PetscErrorCode KSPDestroy_STCG(KSP ksp)
633: {
636: /***************************************************************************/
637: /* Clear composed functions */
638: /***************************************************************************/
641: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",NULL);
642: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",NULL);
643: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",NULL);
645: /***************************************************************************/
646: /* Destroy KSP object. */
647: /***************************************************************************/
649: KSPDestroyDefault(ksp);
650: return(0);
651: }
655: static PetscErrorCode KSPSTCGSetRadius_STCG(KSP ksp, PetscReal radius)
656: {
657: KSP_STCG *cg = (KSP_STCG*)ksp->data;
660: cg->radius = radius;
661: return(0);
662: }
666: static PetscErrorCode KSPSTCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
667: {
668: KSP_STCG *cg = (KSP_STCG*)ksp->data;
671: *norm_d = cg->norm_d;
672: return(0);
673: }
677: static PetscErrorCode KSPSTCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
678: {
679: KSP_STCG *cg = (KSP_STCG*)ksp->data;
682: *o_fcn = cg->o_fcn;
683: return(0);
684: }
688: PetscErrorCode KSPSetFromOptions_STCG(KSP ksp)
689: {
691: KSP_STCG *cg = (KSP_STCG*)ksp->data;
694: PetscOptionsHead("KSP STCG options");
695: PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, NULL);
696: PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
697: PetscOptionsTail();
698: return(0);
699: }
701: /*MC
702: KSPSTCG - Code to run conjugate gradient method subject to a constraint
703: on the solution norm. This is used in Trust Region methods for
704: nonlinear equations, SNESNEWTONTR
706: Options Database Keys:
707: . -ksp_stcg_radius <r> - Trust Region Radius
709: Notes: This is rarely used directly
711: Use preconditioned conjugate gradient to compute
712: an approximate minimizer of the quadratic function
714: q(s) = g^T * s + 0.5 * s^T * H * s
716: subject to the trust region constraint
718: || s || <= delta,
720: where
722: delta is the trust region radius,
723: g is the gradient vector,
724: H is the Hessian approximation, and
725: M is the positive definite preconditioner matrix.
727: KSPConvergedReason may be
728: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
729: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
730: $ other KSP converged/diverged reasons
732: Notes:
733: The preconditioner supplied should be symmetric and positive definite.
735: Level: developer
737: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD(), KSPSTCGGetObjFcn()
738: M*/
742: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
743: {
745: KSP_STCG *cg;
748: PetscNewLog(ksp,&cg);
750: cg->radius = 0.0;
751: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
753: ksp->data = (void*) cg;
754: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
755: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
756: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
758: /***************************************************************************/
759: /* Sets the functions that are associated with this data structure */
760: /* (in C++ this is the same as defining virtual functions). */
761: /***************************************************************************/
763: ksp->ops->setup = KSPSetUp_STCG;
764: ksp->ops->solve = KSPSolve_STCG;
765: ksp->ops->destroy = KSPDestroy_STCG;
766: ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
767: ksp->ops->buildsolution = KSPBuildSolutionDefault;
768: ksp->ops->buildresidual = KSPBuildResidualDefault;
769: ksp->ops->view = 0;
771: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",KSPSTCGSetRadius_STCG);
772: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",KSPSTCGGetNormD_STCG);
773: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",KSPSTCGGetObjFcn_STCG);
774: return(0);
775: }