Actual source code: stcg.c


  2: #include <../src/ksp/ksp/impls/cg/stcg/stcgimpl.h>

  4: #define STCG_PRECONDITIONED_DIRECTION   0
  5: #define STCG_UNPRECONDITIONED_DIRECTION 1
  6: #define STCG_DIRECTION_TYPES            2

  8: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};

 10: static PetscErrorCode KSPCGSolve_STCG(KSP ksp)
 11: {
 12: #if defined(PETSC_USE_COMPLEX)
 13:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "STCG is not available for complex systems");
 14: #else
 15:   KSPCG_STCG     *cg = (KSPCG_STCG*)ksp->data;
 17:   Mat            Qmat, Mmat;
 18:   Vec            r, z, p, d;
 19:   PC             pc;
 20:   PetscReal      norm_r, norm_d, norm_dp1, norm_p, dMp;
 21:   PetscReal      alpha, beta, kappa, rz, rzm1;
 22:   PetscReal      rr, r2, step;
 23:   PetscInt       max_cg_its;
 24:   PetscBool      diagonalscale;

 26:   /***************************************************************************/
 27:   /* Check the arguments and parameters.                                     */
 28:   /***************************************************************************/

 31:   PCGetDiagonalScale(ksp->pc, &diagonalscale);
 32:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
 33:   if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");

 35:   /***************************************************************************/
 36:   /* Get the workspace vectors and initialize variables                      */
 37:   /***************************************************************************/

 39:   r2 = cg->radius * cg->radius;
 40:   r  = ksp->work[0];
 41:   z  = ksp->work[1];
 42:   p  = ksp->work[2];
 43:   d  = ksp->vec_sol;
 44:   pc = ksp->pc;

 46:   PCGetOperators(pc, &Qmat, &Mmat);

 48:   VecGetSize(d, &max_cg_its);
 49:   max_cg_its = PetscMin(max_cg_its, ksp->max_it);
 50:   ksp->its   = 0;

 52:   /***************************************************************************/
 53:   /* Initialize objective function and direction.                            */
 54:   /***************************************************************************/

 56:   cg->o_fcn = 0.0;

 58:   VecSet(d, 0.0);            /* d = 0             */
 59:   cg->norm_d = 0.0;

 61:   /***************************************************************************/
 62:   /* Begin the conjugate gradient method.  Check the right-hand side for     */
 63:   /* numerical problems.  The check for not-a-number and infinite values     */
 64:   /* need be performed only once.                                            */
 65:   /***************************************************************************/

 67:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
 68:   VecDot(r, r, &rr);               /* rr = r^T r        */
 69:   KSPCheckDot(ksp,rr);

 71:   /***************************************************************************/
 72:   /* Check the preconditioner for numerical problems and for positive        */
 73:   /* definiteness.  The check for not-a-number and infinite values need be   */
 74:   /* performed only once.                                                    */
 75:   /***************************************************************************/

 77:   KSP_PCApply(ksp, r, z);          /* z = inv(M) r      */
 78:   VecDot(r, z, &rz);               /* rz = r^T inv(M) r */
 79:   if (PetscIsInfOrNanScalar(rz)) {
 80:     /*************************************************************************/
 81:     /* The preconditioner contains not-a-number or an infinite value.        */
 82:     /* Return the gradient direction intersected with the trust region.      */
 83:     /*************************************************************************/

 85:     ksp->reason = KSP_DIVERGED_NANORINF;
 86:     PetscInfo1(ksp, "KSPCGSolve_STCG: bad preconditioner: rz=%g\n", (double)rz);

 88:     if (cg->radius != 0) {
 89:       if (r2 >= rr) {
 90:         alpha      = 1.0;
 91:         cg->norm_d = PetscSqrtReal(rr);
 92:       } else {
 93:         alpha      = PetscSqrtReal(r2 / rr);
 94:         cg->norm_d = cg->radius;
 95:       }

 97:       VecAXPY(d, alpha, r);        /* d = d + alpha r   */

 99:       /***********************************************************************/
100:       /* Compute objective function.                                         */
101:       /***********************************************************************/

103:       KSP_MatMult(ksp, Qmat, d, z);
104:       VecAYPX(z, -0.5, ksp->vec_rhs);
105:       VecDot(d, z, &cg->o_fcn);
106:       cg->o_fcn = -cg->o_fcn;
107:       ++ksp->its;
108:     }
109:     return(0);
110:   }

112:   if (rz < 0.0) {
113:     /*************************************************************************/
114:     /* The preconditioner is indefinite.  Because this is the first          */
115:     /* and we do not have a direction yet, we use the gradient step.  Note   */
116:     /* that we cannot use the preconditioned norm when computing the step    */
117:     /* because the matrix is indefinite.                                     */
118:     /*************************************************************************/

120:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
121:     PetscInfo1(ksp, "KSPCGSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz);

123:     if (cg->radius != 0.0) {
124:       if (r2 >= rr) {
125:         alpha      = 1.0;
126:         cg->norm_d = PetscSqrtReal(rr);
127:       } else {
128:         alpha      = PetscSqrtReal(r2 / rr);
129:         cg->norm_d = cg->radius;
130:       }

132:       VecAXPY(d, alpha, r);        /* d = d + alpha r   */

134:       /***********************************************************************/
135:       /* Compute objective function.                                         */
136:       /***********************************************************************/

138:       KSP_MatMult(ksp, Qmat, d, z);
139:       VecAYPX(z, -0.5, ksp->vec_rhs);
140:       VecDot(d, z, &cg->o_fcn);
141:       cg->o_fcn = -cg->o_fcn;
142:       ++ksp->its;
143:     }
144:     return(0);
145:   }

147:   /***************************************************************************/
148:   /* As far as we know, the preconditioner is positive semidefinite.         */
149:   /* Compute and log the residual.  Check convergence because this           */
150:   /* initializes things, but do not terminate until at least one conjugate   */
151:   /* gradient iteration has been performed.                                  */
152:   /***************************************************************************/

154:   switch (ksp->normtype) {
155:   case KSP_NORM_PRECONDITIONED:
156:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
157:     break;

159:   case KSP_NORM_UNPRECONDITIONED:
160:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
161:     break;

163:   case KSP_NORM_NATURAL:
164:     norm_r = PetscSqrtReal(rz);                                 /* norm_r = |r|_M    */
165:     break;

167:   default:
168:     norm_r = 0.0;
169:     break;
170:   }

172:   KSPLogResidualHistory(ksp, norm_r);
173:   KSPMonitor(ksp, ksp->its, norm_r);
174:   ksp->rnorm = norm_r;

176:   (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);

178:   /***************************************************************************/
179:   /* Compute the first direction and update the iteration.                   */
180:   /***************************************************************************/

182:   VecCopy(z, p);                   /* p = z             */
183:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
184:   ++ksp->its;

186:   /***************************************************************************/
187:   /* Check the matrix for numerical problems.                                */
188:   /***************************************************************************/

190:   VecDot(p, z, &kappa);            /* kappa = p^T Q p   */
191:   if (PetscIsInfOrNanScalar(kappa)) {
192:     /*************************************************************************/
193:     /* The matrix produced not-a-number or an infinite value.  In this case, */
194:     /* we must stop and use the gradient direction.  This condition need     */
195:     /* only be checked once.                                                 */
196:     /*************************************************************************/

198:     ksp->reason = KSP_DIVERGED_NANORINF;
199:     PetscInfo1(ksp, "KSPCGSolve_STCG: bad matrix: kappa=%g\n", (double)kappa);

201:     if (cg->radius) {
202:       if (r2 >= rr) {
203:         alpha      = 1.0;
204:         cg->norm_d = PetscSqrtReal(rr);
205:       } else {
206:         alpha      = PetscSqrtReal(r2 / rr);
207:         cg->norm_d = cg->radius;
208:       }

210:       VecAXPY(d, alpha, r);        /* d = d + alpha r   */

212:       /***********************************************************************/
213:       /* Compute objective function.                                         */
214:       /***********************************************************************/

216:       KSP_MatMult(ksp, Qmat, d, z);
217:       VecAYPX(z, -0.5, ksp->vec_rhs);
218:       VecDot(d, z, &cg->o_fcn);
219:       cg->o_fcn = -cg->o_fcn;
220:       ++ksp->its;
221:     }
222:     return(0);
223:   }

225:   /***************************************************************************/
226:   /* Initialize variables for calculating the norm of the direction.         */
227:   /***************************************************************************/

229:   dMp    = 0.0;
230:   norm_d = 0.0;
231:   switch (cg->dtype) {
232:   case STCG_PRECONDITIONED_DIRECTION:
233:     norm_p = rz;
234:     break;

236:   default:
237:     VecDot(p, p, &norm_p);
238:     break;
239:   }

241:   /***************************************************************************/
242:   /* Check for negative curvature.                                           */
243:   /***************************************************************************/

245:   if (kappa <= 0.0) {
246:     /*************************************************************************/
247:     /* In this case, the matrix is indefinite and we have encountered a      */
248:     /* direction of negative curvature.  Because negative curvature occurs   */
249:     /* during the first step, we must follow a direction.                    */
250:     /*************************************************************************/

252:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
253:     PetscInfo1(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);

255:     if (cg->radius != 0.0 && norm_p > 0.0) {
256:       /***********************************************************************/
257:       /* Follow direction of negative curvature to the boundary of the       */
258:       /* trust region.                                                       */
259:       /***********************************************************************/

261:       step       = PetscSqrtReal(r2 / norm_p);
262:       cg->norm_d = cg->radius;

264:       VecAXPY(d, step, p); /* d = d + step p    */

266:       /***********************************************************************/
267:       /* Update objective function.                                          */
268:       /***********************************************************************/

270:       cg->o_fcn += step * (0.5 * step * kappa - rz);
271:     } else if (cg->radius != 0.0) {
272:       /***********************************************************************/
273:       /* The norm of the preconditioned direction is zero; use the gradient  */
274:       /* step.                                                               */
275:       /***********************************************************************/

277:       if (r2 >= rr) {
278:         alpha      = 1.0;
279:         cg->norm_d = PetscSqrtReal(rr);
280:       } else {
281:         alpha      = PetscSqrtReal(r2 / rr);
282:         cg->norm_d = cg->radius;
283:       }

285:       VecAXPY(d, alpha, r);        /* d = d + alpha r   */

287:       /***********************************************************************/
288:       /* Compute objective function.                                         */
289:       /***********************************************************************/

291:       KSP_MatMult(ksp, Qmat, d, z);
292:       VecAYPX(z, -0.5, ksp->vec_rhs);
293:       VecDot(d, z, &cg->o_fcn);

295:       cg->o_fcn = -cg->o_fcn;
296:       ++ksp->its;
297:     }
298:     return(0);
299:   }

301:   /***************************************************************************/
302:   /* Run the conjugate gradient method until either the problem is solved,   */
303:   /* we encounter the boundary of the trust region, or the conjugate         */
304:   /* gradient method breaks down.                                            */
305:   /***************************************************************************/

307:   while (1) {
308:     /*************************************************************************/
309:     /* Know that kappa is nonzero, because we have not broken down, so we    */
310:     /* can compute the steplength.                                           */
311:     /*************************************************************************/

313:     alpha = rz / kappa;

315:     /*************************************************************************/
316:     /* Compute the steplength and check for intersection with the trust      */
317:     /* region.                                                               */
318:     /*************************************************************************/

320:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
321:     if (cg->radius != 0.0 && norm_dp1 >= r2) {
322:       /***********************************************************************/
323:       /* In this case, the matrix is positive definite as far as we know.    */
324:       /* However, the full step goes beyond the trust region.                */
325:       /***********************************************************************/

327:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
328:       PetscInfo1(ksp, "KSPCGSolve_STCG: constrained step: radius=%g\n", (double)cg->radius);

330:       if (norm_p > 0.0) {
331:         /*********************************************************************/
332:         /* Follow the direction to the boundary of the trust region.         */
333:         /*********************************************************************/

335:         step       = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
336:         cg->norm_d = cg->radius;

338:         VecAXPY(d, step, p);       /* d = d + step p    */

340:         /*********************************************************************/
341:         /* Update objective function.                                        */
342:         /*********************************************************************/

344:         cg->o_fcn += step * (0.5 * step * kappa - rz);
345:       } else {
346:         /*********************************************************************/
347:         /* The norm of the direction is zero; there is nothing to follow.    */
348:         /*********************************************************************/
349:       }
350:       break;
351:     }

353:     /*************************************************************************/
354:     /* Now we can update the direction and residual.                         */
355:     /*************************************************************************/

357:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
358:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
359:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

361:     switch (cg->dtype) {
362:     case STCG_PRECONDITIONED_DIRECTION:
363:       norm_d = norm_dp1;
364:       break;

366:     default:
367:       VecDot(d, d, &norm_d);
368:       break;
369:     }
370:     cg->norm_d = PetscSqrtReal(norm_d);

372:     /*************************************************************************/
373:     /* Update objective function.                                            */
374:     /*************************************************************************/

376:     cg->o_fcn -= 0.5 * alpha * rz;

378:     /*************************************************************************/
379:     /* Check that the preconditioner appears positive semidefinite.          */
380:     /*************************************************************************/

382:     rzm1 = rz;
383:     VecDot(r, z, &rz);             /* rz = r^T z        */
384:     if (rz < 0.0) {
385:       /***********************************************************************/
386:       /* The preconditioner is indefinite.                                   */
387:       /***********************************************************************/

389:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
390:       PetscInfo1(ksp, "KSPCGSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz);
391:       break;
392:     }

394:     /*************************************************************************/
395:     /* As far as we know, the preconditioner is positive semidefinite.       */
396:     /* Compute the residual and check for convergence.                       */
397:     /*************************************************************************/

399:     switch (ksp->normtype) {
400:     case KSP_NORM_PRECONDITIONED:
401:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
402:       break;

404:     case KSP_NORM_UNPRECONDITIONED:
405:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
406:       break;

408:     case KSP_NORM_NATURAL:
409:       norm_r = PetscSqrtReal(rz);                               /* norm_r = |r|_M    */
410:       break;

412:     default:
413:       norm_r = 0.0;
414:       break;
415:     }

417:     KSPLogResidualHistory(ksp, norm_r);
418:     KSPMonitor(ksp, ksp->its, norm_r);
419:     ksp->rnorm = norm_r;

421:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
422:     if (ksp->reason) {
423:       /***********************************************************************/
424:       /* The method has converged.                                           */
425:       /***********************************************************************/

427:       PetscInfo2(ksp, "KSPCGSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
428:       break;
429:     }

431:     /*************************************************************************/
432:     /* We have not converged yet.  Check for breakdown.                      */
433:     /*************************************************************************/

435:     beta = rz / rzm1;
436:     if (PetscAbsScalar(beta) <= 0.0) {
437:       /***********************************************************************/
438:       /* Conjugate gradients has broken down.                                */
439:       /***********************************************************************/

441:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
442:       PetscInfo1(ksp, "KSPCGSolve_STCG: breakdown: beta=%g\n", (double)beta);
443:       break;
444:     }

446:     /*************************************************************************/
447:     /* Check iteration limit.                                                */
448:     /*************************************************************************/

450:     if (ksp->its >= max_cg_its) {
451:       ksp->reason = KSP_DIVERGED_ITS;
452:       PetscInfo1(ksp, "KSPCGSolve_STCG: iterlim: its=%D\n", ksp->its);
453:       break;
454:     }

456:     /*************************************************************************/
457:     /* Update p and the norms.                                               */
458:     /*************************************************************************/

460:     VecAYPX(p, beta, z);          /* p = z + beta p    */

462:     switch (cg->dtype) {
463:     case STCG_PRECONDITIONED_DIRECTION:
464:       dMp    = beta*(dMp + alpha*norm_p);
465:       norm_p = beta*(rzm1 + beta*norm_p);
466:       break;

468:     default:
469:       VecDot(d, p, &dMp);
470:       VecDot(p, p, &norm_p);
471:       break;
472:     }

474:     /*************************************************************************/
475:     /* Compute the new direction and update the iteration.                   */
476:     /*************************************************************************/

478:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
479:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
480:     ++ksp->its;

482:     /*************************************************************************/
483:     /* Check for negative curvature.                                         */
484:     /*************************************************************************/

486:     if (kappa <= 0.0) {
487:       /***********************************************************************/
488:       /* In this case, the matrix is indefinite and we have encountered      */
489:       /* a direction of negative curvature.  Follow the direction to the     */
490:       /* boundary of the trust region.                                       */
491:       /***********************************************************************/

493:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
494:       PetscInfo1(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);

496:       if (cg->radius != 0.0 && norm_p > 0.0) {
497:         /*********************************************************************/
498:         /* Follow direction of negative curvature to boundary.               */
499:         /*********************************************************************/

501:         step       = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
502:         cg->norm_d = cg->radius;

504:         VecAXPY(d, step, p);       /* d = d + step p    */

506:         /*********************************************************************/
507:         /* Update objective function.                                        */
508:         /*********************************************************************/

510:         cg->o_fcn += step * (0.5 * step * kappa - rz);
511:       } else if (cg->radius != 0.0) {
512:         /*********************************************************************/
513:         /* The norm of the direction is zero; there is nothing to follow.    */
514:         /*********************************************************************/
515:       }
516:       break;
517:     }
518:   }
519:   return(0);
520: #endif
521: }

523: static PetscErrorCode KSPCGSetUp_STCG(KSP ksp)
524: {

528:   /***************************************************************************/
529:   /* Set work vectors needed by conjugate gradient method and allocate       */
530:   /***************************************************************************/

532:   KSPSetWorkVecs(ksp, 3);
533:   return(0);
534: }

536: static PetscErrorCode KSPCGDestroy_STCG(KSP ksp)
537: {

541:   /***************************************************************************/
542:   /* Clear composed functions                                                */
543:   /***************************************************************************/

545:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
546:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
547:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);

549:   /***************************************************************************/
550:   /* Destroy KSP object.                                                     */
551:   /***************************************************************************/

553:   KSPDestroyDefault(ksp);
554:   return(0);
555: }

557: static PetscErrorCode  KSPCGSetRadius_STCG(KSP ksp, PetscReal radius)
558: {
559:   KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;

562:   cg->radius = radius;
563:   return(0);
564: }

566: static PetscErrorCode  KSPCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
567: {
568:   KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;

571:   *norm_d = cg->norm_d;
572:   return(0);
573: }

575: static PetscErrorCode  KSPCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
576: {
577:   KSPCG_STCG *cg = (KSPCG_STCG*)ksp->data;

580:   *o_fcn = cg->o_fcn;
581:   return(0);
582: }

584: static PetscErrorCode KSPCGSetFromOptions_STCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
585: {
587:   KSPCG_STCG     *cg = (KSPCG_STCG*)ksp->data;

590:   PetscOptionsHead(PetscOptionsObject,"KSPCG STCG options");
591:   PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
592:   PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
593:   PetscOptionsTail();
594:   return(0);
595: }

597: /*MC
598:      KSPSTCG -   Code to run conjugate gradient method subject to a constraint
599:          on the solution norm. This is used in Trust Region methods for
600:          nonlinear equations, SNESNEWTONTR

602:    Options Database Keys:
603: .      -ksp_cg_radius <r> - Trust Region Radius

605:    Notes:
606:     This is rarely used directly

608:   Use preconditioned conjugate gradient to compute
609:   an approximate minimizer of the quadratic function

611:             q(s) = g^T * s + 0.5 * s^T * H * s

613:    subject to the trust region constraint

615:             || s || <= delta,

617:    where

619:      delta is the trust region radius,
620:      g is the gradient vector,
621:      H is the Hessian approximation, and
622:      M is the positive definite preconditioner matrix.

624:    KSPConvergedReason may be
625: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
626: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
627: $  other KSP converged/diverged reasons

629:   Notes:
630:   The preconditioner supplied should be symmetric and positive definite.

632:   References:
633:     1. Steihaug, T. (1983): The conjugate gradient method and trust regions in large scale optimization. SIAM J. Numer. Anal. 20, 626--637
634:     2. Toint, Ph.L. (1981): Towards an efficient sparsity exploiting Newton method for minimization. In: Duff, I., ed., Sparse Matrices and Their Uses, pp. 57--88. Academic Press

636:    Level: developer

638: .seealso:  KSPCreate(), KSPCGSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn()
639: M*/

641: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
642: {
644:   KSPCG_STCG     *cg;

647:   PetscNewLog(ksp,&cg);

649:   cg->radius = 0.0;
650:   cg->dtype  = STCG_UNPRECONDITIONED_DIRECTION;

652:   ksp->data  = (void*) cg;
653:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
654:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
655:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
656:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

658:   /***************************************************************************/
659:   /* Sets the functions that are associated with this data structure         */
660:   /* (in C++ this is the same as defining virtual functions).                */
661:   /***************************************************************************/

663:   ksp->ops->setup          = KSPCGSetUp_STCG;
664:   ksp->ops->solve          = KSPCGSolve_STCG;
665:   ksp->ops->destroy        = KSPCGDestroy_STCG;
666:   ksp->ops->setfromoptions = KSPCGSetFromOptions_STCG;
667:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
668:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
669:   ksp->ops->view           = NULL;

671:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_STCG);
672:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",KSPCGGetNormD_STCG);
673:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_STCG);
674:   return(0);
675: }