Actual source code: stcg.c
petsc-3.3-p7 2013-05-11
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(((PetscObject)ksp)->comm,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: #ifdef PETSC_USE_COMPLEX
96: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP, "STCG is not available for complex systems");
97: #else
98: KSP_STCG *cg = (KSP_STCG *)ksp->data;
101: MatStructure pflag;
102: Mat Qmat, Mmat;
103: Vec r, z, p, d;
104: PC pc;
106: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
107: PetscReal alpha, beta, kappa, rz, rzm1;
108: PetscReal rr, r2, step;
110: PetscInt max_cg_its;
112: PetscBool diagonalscale;
115: /***************************************************************************/
116: /* Check the arguments and parameters. */
117: /***************************************************************************/
119: PCGetDiagonalScale(ksp->pc, &diagonalscale);
120: if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
121: if (cg->radius < 0.0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
123: /***************************************************************************/
124: /* Get the workspace vectors and initialize variables */
125: /***************************************************************************/
127: r2 = cg->radius * cg->radius;
128: r = ksp->work[0];
129: z = ksp->work[1];
130: p = ksp->work[2];
131: d = ksp->vec_sol;
132: pc = ksp->pc;
134: PCGetOperators(pc, &Qmat, &Mmat, &pflag);
136: VecGetSize(d, &max_cg_its);
137: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
138: ksp->its = 0;
140: /***************************************************************************/
141: /* Initialize objective function and direction. */
142: /***************************************************************************/
144: cg->o_fcn = 0.0;
146: VecSet(d, 0.0); /* d = 0 */
147: cg->norm_d = 0.0;
149: /***************************************************************************/
150: /* Begin the conjugate gradient method. Check the right-hand side for */
151: /* numerical problems. The check for not-a-number and infinite values */
152: /* need be performed only once. */
153: /***************************************************************************/
155: VecCopy(ksp->vec_rhs, r); /* r = -grad */
156: VecDot(r, r, &rr); /* rr = r^T r */
157: if (PetscIsInfOrNanScalar(rr)) {
158: /*************************************************************************/
159: /* The right-hand side contains not-a-number or an infinite value. */
160: /* The gradient step does not work; return a zero value for the step. */
161: /*************************************************************************/
163: ksp->reason = KSP_DIVERGED_NAN;
164: PetscInfo1(ksp, "KSPSolve_STCG: bad right-hand side: rr=%g\n", rr);
165: return(0);
166: }
168: /***************************************************************************/
169: /* Check the preconditioner for numerical problems and for positive */
170: /* definiteness. The check for not-a-number and infinite values need be */
171: /* performed only once. */
172: /***************************************************************************/
174: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
175: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
176: if (PetscIsInfOrNanScalar(rz)) {
177: /*************************************************************************/
178: /* The preconditioner contains not-a-number or an infinite value. */
179: /* Return the gradient direction intersected with the trust region. */
180: /*************************************************************************/
182: ksp->reason = KSP_DIVERGED_NAN;
183: PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", rz);
185: if (cg->radius != 0) {
186: if (r2 >= rr) {
187: alpha = 1.0;
188: cg->norm_d = PetscSqrtReal(rr);
189: }
190: else {
191: alpha = PetscSqrtReal(r2 / rr);
192: cg->norm_d = cg->radius;
193: }
195: VecAXPY(d, alpha, r); /* d = d + alpha r */
197: /***********************************************************************/
198: /* Compute objective function. */
199: /***********************************************************************/
201: KSP_MatMult(ksp, Qmat, d, z);
202: VecAYPX(z, -0.5, ksp->vec_rhs);
203: VecDot(d, z, &cg->o_fcn);
204: cg->o_fcn = -cg->o_fcn;
205: ++ksp->its;
206: }
207: return(0);
208: }
210: if (rz < 0.0) {
211: /*************************************************************************/
212: /* The preconditioner is indefinite. Because this is the first */
213: /* and we do not have a direction yet, we use the gradient step. Note */
214: /* that we cannot use the preconditioned norm when computing the step */
215: /* because the matrix is indefinite. */
216: /*************************************************************************/
218: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
219: PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", rz);
221: if (cg->radius != 0.0) {
222: if (r2 >= rr) {
223: alpha = 1.0;
224: cg->norm_d = PetscSqrtReal(rr);
225: }
226: else {
227: alpha = PetscSqrtReal(r2 / rr);
228: cg->norm_d = cg->radius;
229: }
231: VecAXPY(d, alpha, r); /* d = d + alpha r */
233: /***********************************************************************/
234: /* Compute objective function. */
235: /***********************************************************************/
237: KSP_MatMult(ksp, Qmat, d, z);
238: VecAYPX(z, -0.5, ksp->vec_rhs);
239: VecDot(d, z, &cg->o_fcn);
240: cg->o_fcn = -cg->o_fcn;
241: ++ksp->its;
242: }
243: return(0);
244: }
246: /***************************************************************************/
247: /* As far as we know, the preconditioner is positive semidefinite. */
248: /* Compute and log the residual. Check convergence because this */
249: /* initializes things, but do not terminate until at least one conjugate */
250: /* gradient iteration has been performed. */
251: /***************************************************************************/
253: switch(ksp->normtype) {
254: case KSP_NORM_PRECONDITIONED:
255: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
256: break;
258: case KSP_NORM_UNPRECONDITIONED:
259: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
260: break;
262: case KSP_NORM_NATURAL:
263: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
264: break;
266: default:
267: norm_r = 0.0;
268: break;
269: }
271: KSPLogResidualHistory(ksp, norm_r);
272: KSPMonitor(ksp, ksp->its, norm_r);
273: ksp->rnorm = norm_r;
275: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
277: /***************************************************************************/
278: /* Compute the first direction and update the iteration. */
279: /***************************************************************************/
281: VecCopy(z, p); /* p = z */
282: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
283: ++ksp->its;
285: /***************************************************************************/
286: /* Check the matrix for numerical problems. */
287: /***************************************************************************/
289: VecDot(p, z, &kappa); /* kappa = p^T Q p */
290: if (PetscIsInfOrNanScalar(kappa)) {
291: /*************************************************************************/
292: /* The matrix produced not-a-number or an infinite value. In this case, */
293: /* we must stop and use the gradient direction. This condition need */
294: /* only be checked once. */
295: /*************************************************************************/
297: ksp->reason = KSP_DIVERGED_NAN;
298: PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", kappa);
300: if (cg->radius) {
301: if (r2 >= rr) {
302: alpha = 1.0;
303: cg->norm_d = PetscSqrtReal(rr);
304: }
305: else {
306: alpha = PetscSqrtReal(r2 / rr);
307: cg->norm_d = cg->radius;
308: }
310: VecAXPY(d, alpha, r); /* d = d + alpha r */
312: /***********************************************************************/
313: /* Compute objective function. */
314: /***********************************************************************/
316: KSP_MatMult(ksp, Qmat, d, z);
317: VecAYPX(z, -0.5, ksp->vec_rhs);
318: VecDot(d, z, &cg->o_fcn);
319: cg->o_fcn = -cg->o_fcn;
320: ++ksp->its;
321: }
322: return(0);
323: }
325: /***************************************************************************/
326: /* Initialize variables for calculating the norm of the direction. */
327: /***************************************************************************/
329: dMp = 0.0;
330: norm_d = 0.0;
331: switch(cg->dtype) {
332: case STCG_PRECONDITIONED_DIRECTION:
333: norm_p = rz;
334: break;
336: default:
337: VecDot(p, p, &norm_p);
338: break;
339: }
341: /***************************************************************************/
342: /* Check for negative curvature. */
343: /***************************************************************************/
345: if (kappa <= 0.0) {
346: /*************************************************************************/
347: /* In this case, the matrix is indefinite and we have encountered a */
348: /* direction of negative curvature. Because negative curvature occurs */
349: /* during the first step, we must follow a direction. */
350: /*************************************************************************/
352: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
353: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", kappa);
355: if (cg->radius != 0.0 && norm_p > 0.0) {
356: /***********************************************************************/
357: /* Follow direction of negative curvature to the boundary of the */
358: /* trust region. */
359: /***********************************************************************/
361: step = PetscSqrtReal(r2 / norm_p);
362: cg->norm_d = cg->radius;
364: VecAXPY(d, step, p); /* d = d + step p */
366: /***********************************************************************/
367: /* Update objective function. */
368: /***********************************************************************/
370: cg->o_fcn += step * (0.5 * step * kappa - rz);
371: }
372: else if (cg->radius != 0.0) {
373: /***********************************************************************/
374: /* The norm of the preconditioned direction is zero; use the gradient */
375: /* step. */
376: /***********************************************************************/
378: if (r2 >= rr) {
379: alpha = 1.0;
380: cg->norm_d = PetscSqrtReal(rr);
381: }
382: else {
383: alpha = PetscSqrtReal(r2 / rr);
384: cg->norm_d = cg->radius;
385: }
387: VecAXPY(d, alpha, r); /* d = d + alpha r */
389: /***********************************************************************/
390: /* Compute objective function. */
391: /***********************************************************************/
393: KSP_MatMult(ksp, Qmat, d, z);
394: VecAYPX(z, -0.5, ksp->vec_rhs);
395: VecDot(d, z, &cg->o_fcn);
396: cg->o_fcn = -cg->o_fcn;
397: ++ksp->its;
398: }
399: return(0);
400: }
402: /***************************************************************************/
403: /* Run the conjugate gradient method until either the problem is solved, */
404: /* we encounter the boundary of the trust region, or the conjugate */
405: /* gradient method breaks down. */
406: /***************************************************************************/
408: while(1) {
409: /*************************************************************************/
410: /* Know that kappa is nonzero, because we have not broken down, so we */
411: /* can compute the steplength. */
412: /*************************************************************************/
414: alpha = rz / kappa;
416: /*************************************************************************/
417: /* Compute the steplength and check for intersection with the trust */
418: /* region. */
419: /*************************************************************************/
421: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
422: if (cg->radius != 0.0 && norm_dp1 >= r2) {
423: /***********************************************************************/
424: /* In this case, the matrix is positive definite as far as we know. */
425: /* However, the full step goes beyond the trust region. */
426: /***********************************************************************/
428: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
429: PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", cg->radius);
431: if (norm_p > 0.0) {
432: /*********************************************************************/
433: /* Follow the direction to the boundary of the trust region. */
434: /*********************************************************************/
436: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
437: cg->norm_d = cg->radius;
439: VecAXPY(d, step, p); /* d = d + step p */
441: /*********************************************************************/
442: /* Update objective function. */
443: /*********************************************************************/
445: cg->o_fcn += step * (0.5 * step * kappa - rz);
446: }
447: else {
448: /*********************************************************************/
449: /* The norm of the direction is zero; there is nothing to follow. */
450: /*********************************************************************/
451: }
452: break;
453: }
455: /*************************************************************************/
456: /* Now we can update the direction and residual. */
457: /*************************************************************************/
459: VecAXPY(d, alpha, p); /* d = d + alpha p */
460: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
461: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
463: switch(cg->dtype) {
464: case STCG_PRECONDITIONED_DIRECTION:
465: norm_d = norm_dp1;
466: break;
468: default:
469: VecDot(d, d, &norm_d);
470: break;
471: }
472: cg->norm_d = PetscSqrtReal(norm_d);
474: /*************************************************************************/
475: /* Update objective function. */
476: /*************************************************************************/
478: cg->o_fcn -= 0.5 * alpha * rz;
480: /*************************************************************************/
481: /* Check that the preconditioner appears positive semidefinite. */
482: /*************************************************************************/
484: rzm1 = rz;
485: VecDot(r, z, &rz); /* rz = r^T z */
486: if (rz < 0.0) {
487: /***********************************************************************/
488: /* The preconditioner is indefinite. */
489: /***********************************************************************/
491: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
492: PetscInfo1(ksp, "KSPSolve_STCG: cg indefinite preconditioner: rz=%g\n", rz);
493: break;
494: }
496: /*************************************************************************/
497: /* As far as we know, the preconditioner is positive semidefinite. */
498: /* Compute the residual and check for convergence. */
499: /*************************************************************************/
501: switch(ksp->normtype) {
502: case KSP_NORM_PRECONDITIONED:
503: VecNorm(z, NORM_2, &norm_r);/* norm_r = |z| */
504: break;
506: case KSP_NORM_UNPRECONDITIONED:
507: VecNorm(r, NORM_2, &norm_r);/* norm_r = |r| */
508: break;
510: case KSP_NORM_NATURAL:
511: norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
512: break;
514: default:
515: norm_r = 0.0;
516: break;
517: }
519: KSPLogResidualHistory(ksp, norm_r);
520: KSPMonitor(ksp, ksp->its, norm_r);
521: ksp->rnorm = norm_r;
522:
523: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
524: if (ksp->reason) {
525: /***********************************************************************/
526: /* The method has converged. */
527: /***********************************************************************/
529: PetscInfo2(ksp, "KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
530: break;
531: }
533: /*************************************************************************/
534: /* We have not converged yet. Check for breakdown. */
535: /*************************************************************************/
537: beta = rz / rzm1;
538: if (fabs(beta) <= 0.0) {
539: /***********************************************************************/
540: /* Conjugate gradients has broken down. */
541: /***********************************************************************/
543: ksp->reason = KSP_DIVERGED_BREAKDOWN;
544: PetscInfo1(ksp, "KSPSolve_STCG: breakdown: beta=%g\n", beta);
545: break;
546: }
548: /*************************************************************************/
549: /* Check iteration limit. */
550: /*************************************************************************/
552: if (ksp->its >= max_cg_its) {
553: ksp->reason = KSP_DIVERGED_ITS;
554: PetscInfo1(ksp, "KSPSolve_STCG: iterlim: its=%d\n", ksp->its);
555: break;
556: }
558: /*************************************************************************/
559: /* Update p and the norms. */
560: /*************************************************************************/
562: VecAYPX(p, beta, z); /* p = z + beta p */
564: switch(cg->dtype) {
565: case STCG_PRECONDITIONED_DIRECTION:
566: dMp = beta*(dMp + alpha*norm_p);
567: norm_p = beta*(rzm1 + beta*norm_p);
568: break;
570: default:
571: VecDot(d, p, &dMp);
572: VecDot(p, p, &norm_p);
573: break;
574: }
576: /*************************************************************************/
577: /* Compute the new direction and update the iteration. */
578: /*************************************************************************/
580: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
581: VecDot(p, z, &kappa); /* kappa = p^T Q p */
582: ++ksp->its;
584: /*************************************************************************/
585: /* Check for negative curvature. */
586: /*************************************************************************/
588: if (kappa <= 0.0) {
589: /***********************************************************************/
590: /* In this case, the matrix is indefinite and we have encountered */
591: /* a direction of negative curvature. Follow the direction to the */
592: /* boundary of the trust region. */
593: /***********************************************************************/
595: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
596: PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", kappa);
598: if (cg->radius != 0.0 && norm_p > 0.0) {
599: /*********************************************************************/
600: /* Follow direction of negative curvature to boundary. */
601: /*********************************************************************/
603: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
604: cg->norm_d = cg->radius;
606: VecAXPY(d, step, p); /* d = d + step p */
608: /*********************************************************************/
609: /* Update objective function. */
610: /*********************************************************************/
612: cg->o_fcn += step * (0.5 * step * kappa - rz);
613: }
614: else if (cg->radius != 0.0) {
615: /*********************************************************************/
616: /* The norm of the direction is zero; there is nothing to follow. */
617: /*********************************************************************/
618: }
619: break;
620: }
621: }
623: return(0);
624: #endif
625: }
629: PetscErrorCode KSPSetUp_STCG(KSP ksp)
630: {
634: /***************************************************************************/
635: /* Set work vectors needed by conjugate gradient method and allocate */
636: /***************************************************************************/
638: KSPDefaultGetWork(ksp, 3);
639: return(0);
640: }
644: PetscErrorCode KSPDestroy_STCG(KSP ksp)
645: {
650: /***************************************************************************/
651: /* Clear composed functions */
652: /***************************************************************************/
654: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGSetRadius_C","",PETSC_NULL);
655: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGGetNormD_C","",PETSC_NULL);
656: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGGetObjFcn_C","",PETSC_NULL);
658: /***************************************************************************/
659: /* Destroy KSP object. */
660: /***************************************************************************/
662: KSPDefaultDestroy(ksp);
663: return(0);
664: }
666: EXTERN_C_BEGIN
669: PetscErrorCode KSPSTCGSetRadius_STCG(KSP ksp, PetscReal radius)
670: {
671: KSP_STCG *cg = (KSP_STCG *)ksp->data;
674: cg->radius = radius;
675: return(0);
676: }
680: PetscErrorCode KSPSTCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
681: {
682: KSP_STCG *cg = (KSP_STCG *)ksp->data;
685: *norm_d = cg->norm_d;
686: return(0);
687: }
691: PetscErrorCode KSPSTCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn){
692: KSP_STCG *cg = (KSP_STCG *)ksp->data;
695: *o_fcn = cg->o_fcn;
696: return(0);
697: }
698: EXTERN_C_END
702: PetscErrorCode KSPSetFromOptions_STCG(KSP ksp)
703: {
705: KSP_STCG *cg = (KSP_STCG *)ksp->data;
708: PetscOptionsHead("KSP STCG options");
709: PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, PETSC_NULL);
710: PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, PETSC_NULL);
711: PetscOptionsTail();
712: return(0);
713: }
715: /*MC
716: KSPSTCG - Code to run conjugate gradient method subject to a constraint
717: on the solution norm. This is used in Trust Region methods for
718: nonlinear equations, SNESTR
720: Options Database Keys:
721: . -ksp_stcg_radius <r> - Trust Region Radius
723: Notes: This is rarely used directly
725: Use preconditioned conjugate gradient to compute
726: an approximate minimizer of the quadratic function
728: q(s) = g^T * s + 0.5 * s^T * H * s
730: subject to the trust region constraint
732: || s || <= delta,
734: where
736: delta is the trust region radius,
737: g is the gradient vector,
738: H is the Hessian approximation, and
739: M is the positive definite preconditioner matrix.
741: KSPConvergedReason may be
742: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
743: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
744: $ other KSP converged/diverged reasons
746: Notes:
747: The preconditioner supplied should be symmetric and positive definite.
749: Level: developer
751: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD(), KSPSTCGGetObjFcn()
752: M*/
754: EXTERN_C_BEGIN
757: PetscErrorCode KSPCreate_STCG(KSP ksp)
758: {
760: KSP_STCG *cg;
764: PetscNewLog(ksp,KSP_STCG, &cg);
766: cg->radius = 0.0;
767: cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
769: ksp->data = (void *) cg;
770: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
771: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
772: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
774: /***************************************************************************/
775: /* Sets the functions that are associated with this data structure */
776: /* (in C++ this is the same as defining virtual functions). */
777: /***************************************************************************/
779: ksp->ops->setup = KSPSetUp_STCG;
780: ksp->ops->solve = KSPSolve_STCG;
781: ksp->ops->destroy = KSPDestroy_STCG;
782: ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
783: ksp->ops->buildsolution = KSPDefaultBuildSolution;
784: ksp->ops->buildresidual = KSPDefaultBuildResidual;
785: ksp->ops->view = 0;
787: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
788: "KSPSTCGSetRadius_C",
789: "KSPSTCGSetRadius_STCG",
790: KSPSTCGSetRadius_STCG);
791: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
792: "KSPSTCGGetNormD_C",
793: "KSPSTCGGetNormD_STCG",
794: KSPSTCGGetNormD_STCG);
795: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
796: "KSPSTCGGetObjFcn_C",
797: "KSPSTCGGetObjFcn_STCG",
798: KSPSTCGGetObjFcn_STCG);
799: return(0);
800: }
801: EXTERN_C_END