Actual source code: stcg.c
petsc-3.4.5 2014-06-29
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: PetscUseMethod(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: MatStructure pflag;
101: Mat Qmat, Mmat;
102: Vec r, z, p, d;
103: PC pc;
104: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
105: PetscReal alpha, beta, kappa, rz, rzm1;
106: PetscReal rr, r2, step;
107: PetscInt max_cg_its;
108: PetscBool diagonalscale;
110: /***************************************************************************/
111: /* Check the arguments and parameters. */
112: /***************************************************************************/
115: PCGetDiagonalScale(ksp->pc, &diagonalscale);
116: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
117: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
119: /***************************************************************************/
120: /* Get the workspace vectors and initialize variables */
121: /***************************************************************************/
123: r2 = cg->radius * cg->radius;
124: r = ksp->work[0];
125: z = ksp->work[1];
126: p = ksp->work[2];
127: d = ksp->vec_sol;
128: pc = ksp->pc;
130: PCGetOperators(pc, &Qmat, &Mmat, &pflag);
132: VecGetSize(d, &max_cg_its);
133: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
134: ksp->its = 0;
136: /***************************************************************************/
137: /* Initialize objective function and direction. */
138: /***************************************************************************/
140: cg->o_fcn = 0.0;
142: VecSet(d, 0.0); /* d = 0 */
143: cg->norm_d = 0.0;
145: /***************************************************************************/
146: /* Begin the conjugate gradient method. Check the right-hand side for */
147: /* numerical problems. The check for not-a-number and infinite values */
148: /* need be performed only once. */
149: /***************************************************************************/
151: VecCopy(ksp->vec_rhs, r); /* r = -grad */
152: VecDot(r, r, &rr); /* rr = r^T r */
153: if (PetscIsInfOrNanScalar(rr)) {
154: /*************************************************************************/
155: /* The right-hand side contains not-a-number or an infinite value. */
156: /* The gradient step does not work; return a zero value for the step. */
157: /*************************************************************************/
159: ksp->reason = KSP_DIVERGED_NANORINF;
160: PetscInfo1(ksp, "KSPSolve_STCG: bad right-hand side: rr=%g\n", rr);
161: return(0);
162: }
164: /***************************************************************************/
165: /* Check the preconditioner for numerical problems and for positive */
166: /* definiteness. The check for not-a-number and infinite values need be */
167: /* performed only once. */
168: /***************************************************************************/
170: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
171: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
172: if (PetscIsInfOrNanScalar(rz)) {
173: /*************************************************************************/
174: /* The preconditioner contains not-a-number or an infinite value. */
175: /* Return the gradient direction intersected with the trust region. */
176: /*************************************************************************/
178: ksp->reason = KSP_DIVERGED_NANORINF;
179: PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", rz);
181: if (cg->radius != 0) {
182: if (r2 >= rr) {
183: alpha = 1.0;
184: cg->norm_d = PetscSqrtReal(rr);
185: } else {
186: alpha = PetscSqrtReal(r2 / rr);
187: cg->norm_d = cg->radius;
188: }
190: VecAXPY(d, alpha, r); /* d = d + alpha r */
192: /***********************************************************************/
193: /* Compute objective function. */
194: /***********************************************************************/
196: KSP_MatMult(ksp, Qmat, d, z);
197: VecAYPX(z, -0.5, ksp->vec_rhs);
198: VecDot(d, z, &cg->o_fcn);
199: cg->o_fcn = -cg->o_fcn;
200: ++ksp->its;
201: }
202: return(0);
203: }
205: if (rz < 0.0) {
206: /*************************************************************************/
207: /* The preconditioner is indefinite. Because this is the first */
208: /* and we do not have a direction yet, we use the gradient step. Note */
209: /* that we cannot use the preconditioned norm when computing the step */
210: /* because the matrix is indefinite. */
211: /*************************************************************************/
213: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
214: PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", rz);
216: if (cg->radius != 0.0) {
217: if (r2 >= rr) {
218: alpha = 1.0;
219: cg->norm_d = PetscSqrtReal(rr);
220: } else {
221: alpha = PetscSqrtReal(r2 / rr);
222: cg->norm_d = cg->radius;
223: }
225: VecAXPY(d, alpha, r); /* d = d + alpha r */
227: /***********************************************************************/
228: /* Compute objective function. */
229: /***********************************************************************/
231: KSP_MatMult(ksp, Qmat, d, z);
232: VecAYPX(z, -0.5, ksp->vec_rhs);
233: VecDot(d, z, &cg->o_fcn);
234: cg->o_fcn = -cg->o_fcn;
235: ++ksp->its;
236: }
237: return(0);
238: }
240: /***************************************************************************/
241: /* As far as we know, the preconditioner is positive semidefinite. */
242: /* Compute and log the residual. Check convergence because this */
243: /* initializes things, but do not terminate until at least one conjugate */
244: /* gradient iteration has been performed. */
245: /***************************************************************************/
247: switch (ksp->normtype) {
248: case KSP_NORM_PRECONDITIONED:
249: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
250: break;
252: case KSP_NORM_UNPRECONDITIONED:
253: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
254: break;
256: case KSP_NORM_NATURAL:
257: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
258: break;
260: default:
261: norm_r = 0.0;
262: break;
263: }
265: KSPLogResidualHistory(ksp, norm_r);
266: KSPMonitor(ksp, ksp->its, norm_r);
267: ksp->rnorm = norm_r;
269: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
271: /***************************************************************************/
272: /* Compute the first direction and update the iteration. */
273: /***************************************************************************/
275: VecCopy(z, p); /* p = z */
276: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
277: ++ksp->its;
279: /***************************************************************************/
280: /* Check the matrix for numerical problems. */
281: /***************************************************************************/
283: VecDot(p, z, &kappa); /* kappa = p^T Q p */
284: if (PetscIsInfOrNanScalar(kappa)) {
285: /*************************************************************************/
286: /* The matrix produced not-a-number or an infinite value. In this case, */
287: /* we must stop and use the gradient direction. This condition need */
288: /* only be checked once. */
289: /*************************************************************************/
291: ksp->reason = KSP_DIVERGED_NANORINF;
292: PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", kappa);
294: if (cg->radius) {
295: if (r2 >= rr) {
296: alpha = 1.0;
297: cg->norm_d = PetscSqrtReal(rr);
298: } else {
299: alpha = PetscSqrtReal(r2 / rr);
300: cg->norm_d = cg->radius;
301: }
303: VecAXPY(d, alpha, r); /* d = d + alpha r */
305: /***********************************************************************/
306: /* Compute objective function. */
307: /***********************************************************************/
309: KSP_MatMult(ksp, Qmat, d, z);
310: VecAYPX(z, -0.5, ksp->vec_rhs);
311: VecDot(d, z, &cg->o_fcn);
312: cg->o_fcn = -cg->o_fcn;
313: ++ksp->its;
314: }
315: return(0);
316: }
318: /***************************************************************************/
319: /* Initialize variables for calculating the norm of the direction. */
320: /***************************************************************************/
322: dMp = 0.0;
323: norm_d = 0.0;
324: switch (cg->dtype) {
325: case STCG_PRECONDITIONED_DIRECTION:
326: norm_p = rz;
327: break;
329: default:
330: VecDot(p, p, &norm_p);
331: break;
332: }
334: /***************************************************************************/
335: /* Check for negative curvature. */
336: /***************************************************************************/
338: if (kappa <= 0.0) {
339: /*************************************************************************/
340: /* In this case, the matrix is indefinite and we have encountered a */
341: /* direction of negative curvature. Because negative curvature occurs */
342: /* during the first step, we must follow a direction. */
343: /*************************************************************************/
345: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
346: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", kappa);
348: if (cg->radius != 0.0 && norm_p > 0.0) {
349: /***********************************************************************/
350: /* Follow direction of negative curvature to the boundary of the */
351: /* trust region. */
352: /***********************************************************************/
354: step = PetscSqrtReal(r2 / norm_p);
355: cg->norm_d = cg->radius;
357: VecAXPY(d, step, p); /* d = d + step p */
359: /***********************************************************************/
360: /* Update objective function. */
361: /***********************************************************************/
363: cg->o_fcn += step * (0.5 * step * kappa - rz);
364: } else if (cg->radius != 0.0) {
365: /***********************************************************************/
366: /* The norm of the preconditioned direction is zero; use the gradient */
367: /* step. */
368: /***********************************************************************/
370: if (r2 >= rr) {
371: alpha = 1.0;
372: cg->norm_d = PetscSqrtReal(rr);
373: } else {
374: alpha = PetscSqrtReal(r2 / rr);
375: cg->norm_d = cg->radius;
376: }
378: VecAXPY(d, alpha, r); /* d = d + alpha r */
380: /***********************************************************************/
381: /* Compute objective function. */
382: /***********************************************************************/
384: KSP_MatMult(ksp, Qmat, d, z);
385: VecAYPX(z, -0.5, ksp->vec_rhs);
386: VecDot(d, z, &cg->o_fcn);
388: cg->o_fcn = -cg->o_fcn;
389: ++ksp->its;
390: }
391: return(0);
392: }
394: /***************************************************************************/
395: /* Run the conjugate gradient method until either the problem is solved, */
396: /* we encounter the boundary of the trust region, or the conjugate */
397: /* gradient method breaks down. */
398: /***************************************************************************/
400: while (1) {
401: /*************************************************************************/
402: /* Know that kappa is nonzero, because we have not broken down, so we */
403: /* can compute the steplength. */
404: /*************************************************************************/
406: alpha = rz / kappa;
408: /*************************************************************************/
409: /* Compute the steplength and check for intersection with the trust */
410: /* region. */
411: /*************************************************************************/
413: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
414: if (cg->radius != 0.0 && norm_dp1 >= r2) {
415: /***********************************************************************/
416: /* In this case, the matrix is positive definite as far as we know. */
417: /* However, the full step goes beyond the trust region. */
418: /***********************************************************************/
420: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
421: PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", cg->radius);
423: if (norm_p > 0.0) {
424: /*********************************************************************/
425: /* Follow the direction to the boundary of the trust region. */
426: /*********************************************************************/
428: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
429: cg->norm_d = cg->radius;
431: VecAXPY(d, step, p); /* d = d + step p */
433: /*********************************************************************/
434: /* Update objective function. */
435: /*********************************************************************/
437: cg->o_fcn += step * (0.5 * step * kappa - rz);
438: } else {
439: /*********************************************************************/
440: /* The norm of the direction is zero; there is nothing to follow. */
441: /*********************************************************************/
442: }
443: break;
444: }
446: /*************************************************************************/
447: /* Now we can update the direction and residual. */
448: /*************************************************************************/
450: VecAXPY(d, alpha, p); /* d = d + alpha p */
451: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
452: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
454: switch (cg->dtype) {
455: case STCG_PRECONDITIONED_DIRECTION:
456: norm_d = norm_dp1;
457: break;
459: default:
460: VecDot(d, d, &norm_d);
461: break;
462: }
463: cg->norm_d = PetscSqrtReal(norm_d);
465: /*************************************************************************/
466: /* Update objective function. */
467: /*************************************************************************/
469: cg->o_fcn -= 0.5 * alpha * rz;
471: /*************************************************************************/
472: /* Check that the preconditioner appears positive semidefinite. */
473: /*************************************************************************/
475: rzm1 = rz;
476: VecDot(r, z, &rz); /* rz = r^T z */
477: if (rz < 0.0) {
478: /***********************************************************************/
479: /* The preconditioner is indefinite. */
480: /***********************************************************************/
482: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
483: PetscInfo1(ksp, "KSPSolve_STCG: cg indefinite preconditioner: rz=%g\n", rz);
484: break;
485: }
487: /*************************************************************************/
488: /* As far as we know, the preconditioner is positive semidefinite. */
489: /* Compute the residual and check for convergence. */
490: /*************************************************************************/
492: switch (ksp->normtype) {
493: case KSP_NORM_PRECONDITIONED:
494: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
495: break;
497: case KSP_NORM_UNPRECONDITIONED:
498: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
499: break;
501: case KSP_NORM_NATURAL:
502: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
503: break;
505: default:
506: norm_r = 0.0;
507: break;
508: }
510: KSPLogResidualHistory(ksp, norm_r);
511: KSPMonitor(ksp, ksp->its, norm_r);
512: ksp->rnorm = norm_r;
514: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
515: if (ksp->reason) {
516: /***********************************************************************/
517: /* The method has converged. */
518: /***********************************************************************/
520: PetscInfo2(ksp, "KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
521: break;
522: }
524: /*************************************************************************/
525: /* We have not converged yet. Check for breakdown. */
526: /*************************************************************************/
528: beta = rz / rzm1;
529: if (fabs(beta) <= 0.0) {
530: /***********************************************************************/
531: /* Conjugate gradients has broken down. */
532: /***********************************************************************/
534: ksp->reason = KSP_DIVERGED_BREAKDOWN;
535: PetscInfo1(ksp, "KSPSolve_STCG: breakdown: beta=%g\n", beta);
536: break;
537: }
539: /*************************************************************************/
540: /* Check iteration limit. */
541: /*************************************************************************/
543: if (ksp->its >= max_cg_its) {
544: ksp->reason = KSP_DIVERGED_ITS;
545: PetscInfo1(ksp, "KSPSolve_STCG: iterlim: its=%d\n", ksp->its);
546: break;
547: }
549: /*************************************************************************/
550: /* Update p and the norms. */
551: /*************************************************************************/
553: VecAYPX(p, beta, z); /* p = z + beta p */
555: switch (cg->dtype) {
556: case STCG_PRECONDITIONED_DIRECTION:
557: dMp = beta*(dMp + alpha*norm_p);
558: norm_p = beta*(rzm1 + beta*norm_p);
559: break;
561: default:
562: VecDot(d, p, &dMp);
563: VecDot(p, p, &norm_p);
564: break;
565: }
567: /*************************************************************************/
568: /* Compute the new direction and update the iteration. */
569: /*************************************************************************/
571: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
572: VecDot(p, z, &kappa); /* kappa = p^T Q p */
573: ++ksp->its;
575: /*************************************************************************/
576: /* Check for negative curvature. */
577: /*************************************************************************/
579: if (kappa <= 0.0) {
580: /***********************************************************************/
581: /* In this case, the matrix is indefinite and we have encountered */
582: /* a direction of negative curvature. Follow the direction to the */
583: /* boundary of the trust region. */
584: /***********************************************************************/
586: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
587: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", kappa);
589: if (cg->radius != 0.0 && norm_p > 0.0) {
590: /*********************************************************************/
591: /* Follow direction of negative curvature to boundary. */
592: /*********************************************************************/
594: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
595: cg->norm_d = cg->radius;
597: VecAXPY(d, step, p); /* d = d + step p */
599: /*********************************************************************/
600: /* Update objective function. */
601: /*********************************************************************/
603: cg->o_fcn += step * (0.5 * step * kappa - rz);
604: } else if (cg->radius != 0.0) {
605: /*********************************************************************/
606: /* The norm of the direction is zero; there is nothing to follow. */
607: /*********************************************************************/
608: }
609: break;
610: }
611: }
612: return(0);
613: #endif
614: }
618: PetscErrorCode KSPSetUp_STCG(KSP ksp)
619: {
622: /***************************************************************************/
623: /* Set work vectors needed by conjugate gradient method and allocate */
624: /***************************************************************************/
627: KSPSetWorkVecs(ksp, 3);
628: return(0);
629: }
633: PetscErrorCode KSPDestroy_STCG(KSP ksp)
634: {
637: /***************************************************************************/
638: /* Clear composed functions */
639: /***************************************************************************/
642: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",NULL);
643: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",NULL);
644: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",NULL);
646: /***************************************************************************/
647: /* Destroy KSP object. */
648: /***************************************************************************/
650: KSPDestroyDefault(ksp);
651: return(0);
652: }
656: static PetscErrorCode KSPSTCGSetRadius_STCG(KSP ksp, PetscReal radius)
657: {
658: KSP_STCG *cg = (KSP_STCG*)ksp->data;
661: cg->radius = radius;
662: return(0);
663: }
667: static PetscErrorCode KSPSTCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
668: {
669: KSP_STCG *cg = (KSP_STCG*)ksp->data;
672: *norm_d = cg->norm_d;
673: return(0);
674: }
678: static PetscErrorCode KSPSTCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
679: {
680: KSP_STCG *cg = (KSP_STCG*)ksp->data;
683: *o_fcn = cg->o_fcn;
684: return(0);
685: }
689: PetscErrorCode KSPSetFromOptions_STCG(KSP ksp)
690: {
692: KSP_STCG *cg = (KSP_STCG*)ksp->data;
695: PetscOptionsHead("KSP STCG options");
696: PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, NULL);
697: PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
698: PetscOptionsTail();
699: return(0);
700: }
702: /*MC
703: KSPSTCG - Code to run conjugate gradient method subject to a constraint
704: on the solution norm. This is used in Trust Region methods for
705: nonlinear equations, SNESNEWTONTR
707: Options Database Keys:
708: . -ksp_stcg_radius <r> - Trust Region Radius
710: Notes: This is rarely used directly
712: Use preconditioned conjugate gradient to compute
713: an approximate minimizer of the quadratic function
715: q(s) = g^T * s + 0.5 * s^T * H * s
717: subject to the trust region constraint
719: || s || <= delta,
721: where
723: delta is the trust region radius,
724: g is the gradient vector,
725: H is the Hessian approximation, and
726: M is the positive definite preconditioner matrix.
728: KSPConvergedReason may be
729: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
730: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
731: $ other KSP converged/diverged reasons
733: Notes:
734: The preconditioner supplied should be symmetric and positive definite.
736: Level: developer
738: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD(), KSPSTCGGetObjFcn()
739: M*/
743: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
744: {
746: KSP_STCG *cg;
749: PetscNewLog(ksp,KSP_STCG, &cg);
751: cg->radius = 0.0;
752: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
754: ksp->data = (void*) cg;
755: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
756: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
757: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
759: /***************************************************************************/
760: /* Sets the functions that are associated with this data structure */
761: /* (in C++ this is the same as defining virtual functions). */
762: /***************************************************************************/
764: ksp->ops->setup = KSPSetUp_STCG;
765: ksp->ops->solve = KSPSolve_STCG;
766: ksp->ops->destroy = KSPDestroy_STCG;
767: ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
768: ksp->ops->buildsolution = KSPBuildSolutionDefault;
769: ksp->ops->buildresidual = KSPBuildResidualDefault;
770: ksp->ops->view = 0;
772: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",KSPSTCGSetRadius_STCG);
773: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",KSPSTCGGetNormD_STCG);
774: PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",KSPSTCGGetObjFcn_STCG);
775: return(0);
776: }