Actual source code: stcg.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: #include <petsc/private/kspimpl.h>             /*I "petscksp.h" I*/
  3: #include <../src/ksp/ksp/impls/cg/stcg/stcgimpl.h>

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

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

 13: /*@
 14:     KSPSTCGSetRadius - Sets the radius of the trust region.

 16:     Logically Collective on KSP

 18:     Input Parameters:
 19: +   ksp    - the iterative context
 20: -   radius - the trust region radius (Infinity is the default)

 22:     Options Database Key:
 23: .   -ksp_stcg_radius <r>

 25:     Level: advanced

 27: .keywords: KSP, STCG, set, trust region radius
 28: @*/
 29: PetscErrorCode  KSPSTCGSetRadius(KSP ksp, PetscReal radius)
 30: {

 35:   if (radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Radius negative");
 37:   PetscTryMethod(ksp,"KSPSTCGSetRadius_C",(KSP,PetscReal),(ksp,radius));
 38:   return(0);
 39: }

 43: /*@
 44:     KSPSTCGGetNormD - Got norm of the direction.

 46:     Collective on KSP

 48:     Input Parameters:
 49: +   ksp    - the iterative context
 50: -   norm_d - the norm of the direction

 52:     Level: advanced

 54: .keywords: KSP, STCG, get, norm direction
 55: @*/
 56: PetscErrorCode  KSPSTCGGetNormD(KSP ksp, PetscReal *norm_d)
 57: {

 62:   PetscUseMethod(ksp,"KSPSTCGGetNormD_C",(KSP,PetscReal*),(ksp,norm_d));
 63:   return(0);
 64: }

 68: /*@
 69:     KSPSTCGGetObjFcn - Get objective function value.

 71:     Collective on KSP

 73:     Input Parameters:
 74: +   ksp   - the iterative context
 75: -   o_fcn - the objective function value

 77:     Level: advanced

 79: .keywords: KSP, STCG, get, objective function
 80: @*/
 81: PetscErrorCode  KSPSTCGGetObjFcn(KSP ksp, PetscReal *o_fcn)
 82: {

 87:   PetscUseMethod(ksp,"KSPSTCGGetObjFcn_C",(KSP,PetscReal*),(ksp,o_fcn));
 88:   return(0);
 89: }

 93: PetscErrorCode KSPSolve_STCG(KSP ksp)
 94: {
 95: #if defined(PETSC_USE_COMPLEX)
 96:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "STCG is not available for complex systems");
 97: #else
 98:   KSP_STCG       *cg = (KSP_STCG*)ksp->data;
100:   Mat            Qmat, Mmat;
101:   Vec            r, z, p, d;
102:   PC             pc;
103:   PetscReal      norm_r, norm_d, norm_dp1, norm_p, dMp;
104:   PetscReal      alpha, beta, kappa, rz, rzm1;
105:   PetscReal      rr, r2, step;
106:   PetscInt       max_cg_its;
107:   PetscBool      diagonalscale;

109:   /***************************************************************************/
110:   /* Check the arguments and parameters.                                     */
111:   /***************************************************************************/

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

118:   /***************************************************************************/
119:   /* Get the workspace vectors and initialize variables                      */
120:   /***************************************************************************/

122:   r2 = cg->radius * cg->radius;
123:   r  = ksp->work[0];
124:   z  = ksp->work[1];
125:   p  = ksp->work[2];
126:   d  = ksp->vec_sol;
127:   pc = ksp->pc;

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

131:   VecGetSize(d, &max_cg_its);
132:   max_cg_its = PetscMin(max_cg_its, ksp->max_it);
133:   ksp->its   = 0;

135:   /***************************************************************************/
136:   /* Initialize objective function and direction.                            */
137:   /***************************************************************************/

139:   cg->o_fcn = 0.0;

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

144:   /***************************************************************************/
145:   /* Begin the conjugate gradient method.  Check the right-hand side for     */
146:   /* numerical problems.  The check for not-a-number and infinite values     */
147:   /* need be performed only once.                                            */
148:   /***************************************************************************/

150:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
151:   VecDot(r, r, &rr);               /* rr = r^T r        */
152:   KSPCheckDot(ksp,rr);

154:   /***************************************************************************/
155:   /* Check the preconditioner for numerical problems and for positive        */
156:   /* definiteness.  The check for not-a-number and infinite values need be   */
157:   /* performed only once.                                                    */
158:   /***************************************************************************/

160:   KSP_PCApply(ksp, r, z);          /* z = inv(M) r      */
161:   VecDot(r, z, &rz);               /* rz = r^T inv(M) r */
162:   if (PetscIsInfOrNanScalar(rz)) {
163:     /*************************************************************************/
164:     /* The preconditioner contains not-a-number or an infinite value.        */
165:     /* Return the gradient direction intersected with the trust region.      */
166:     /*************************************************************************/

168:     ksp->reason = KSP_DIVERGED_NANORINF;
169:     PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", (double)rz);

171:     if (cg->radius != 0) {
172:       if (r2 >= rr) {
173:         alpha      = 1.0;
174:         cg->norm_d = PetscSqrtReal(rr);
175:       } else {
176:         alpha      = PetscSqrtReal(r2 / rr);
177:         cg->norm_d = cg->radius;
178:       }

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

182:       /***********************************************************************/
183:       /* Compute objective function.                                         */
184:       /***********************************************************************/

186:       KSP_MatMult(ksp, Qmat, d, z);
187:       VecAYPX(z, -0.5, ksp->vec_rhs);
188:       VecDot(d, z, &cg->o_fcn);
189:       cg->o_fcn = -cg->o_fcn;
190:       ++ksp->its;
191:     }
192:     return(0);
193:   }

195:   if (rz < 0.0) {
196:     /*************************************************************************/
197:     /* The preconditioner is indefinite.  Because this is the first          */
198:     /* and we do not have a direction yet, we use the gradient step.  Note   */
199:     /* that we cannot use the preconditioned norm when computing the step    */
200:     /* because the matrix is indefinite.                                     */
201:     /*************************************************************************/

203:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
204:     PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz);

206:     if (cg->radius != 0.0) {
207:       if (r2 >= rr) {
208:         alpha      = 1.0;
209:         cg->norm_d = PetscSqrtReal(rr);
210:       } else {
211:         alpha      = PetscSqrtReal(r2 / rr);
212:         cg->norm_d = cg->radius;
213:       }

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

217:       /***********************************************************************/
218:       /* Compute objective function.                                         */
219:       /***********************************************************************/

221:       KSP_MatMult(ksp, Qmat, d, z);
222:       VecAYPX(z, -0.5, ksp->vec_rhs);
223:       VecDot(d, z, &cg->o_fcn);
224:       cg->o_fcn = -cg->o_fcn;
225:       ++ksp->its;
226:     }
227:     return(0);
228:   }

230:   /***************************************************************************/
231:   /* As far as we know, the preconditioner is positive semidefinite.         */
232:   /* Compute and log the residual.  Check convergence because this           */
233:   /* initializes things, but do not terminate until at least one conjugate   */
234:   /* gradient iteration has been performed.                                  */
235:   /***************************************************************************/

237:   switch (ksp->normtype) {
238:   case KSP_NORM_PRECONDITIONED:
239:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
240:     break;

242:   case KSP_NORM_UNPRECONDITIONED:
243:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
244:     break;

246:   case KSP_NORM_NATURAL:
247:     norm_r = PetscSqrtReal(rz);                                 /* norm_r = |r|_M    */
248:     break;

250:   default:
251:     norm_r = 0.0;
252:     break;
253:   }

255:   KSPLogResidualHistory(ksp, norm_r);
256:   KSPMonitor(ksp, ksp->its, norm_r);
257:   ksp->rnorm = norm_r;

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

261:   /***************************************************************************/
262:   /* Compute the first direction and update the iteration.                   */
263:   /***************************************************************************/

265:   VecCopy(z, p);                   /* p = z             */
266:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
267:   ++ksp->its;

269:   /***************************************************************************/
270:   /* Check the matrix for numerical problems.                                */
271:   /***************************************************************************/

273:   VecDot(p, z, &kappa);            /* kappa = p^T Q p   */
274:   if (PetscIsInfOrNanScalar(kappa)) {
275:     /*************************************************************************/
276:     /* The matrix produced not-a-number or an infinite value.  In this case, */
277:     /* we must stop and use the gradient direction.  This condition need     */
278:     /* only be checked once.                                                 */
279:     /*************************************************************************/

281:     ksp->reason = KSP_DIVERGED_NANORINF;
282:     PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", (double)kappa);

284:     if (cg->radius) {
285:       if (r2 >= rr) {
286:         alpha      = 1.0;
287:         cg->norm_d = PetscSqrtReal(rr);
288:       } else {
289:         alpha      = PetscSqrtReal(r2 / rr);
290:         cg->norm_d = cg->radius;
291:       }

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

295:       /***********************************************************************/
296:       /* Compute objective function.                                         */
297:       /***********************************************************************/

299:       KSP_MatMult(ksp, Qmat, d, z);
300:       VecAYPX(z, -0.5, ksp->vec_rhs);
301:       VecDot(d, z, &cg->o_fcn);
302:       cg->o_fcn = -cg->o_fcn;
303:       ++ksp->its;
304:     }
305:     return(0);
306:   }

308:   /***************************************************************************/
309:   /* Initialize variables for calculating the norm of the direction.         */
310:   /***************************************************************************/

312:   dMp    = 0.0;
313:   norm_d = 0.0;
314:   switch (cg->dtype) {
315:   case STCG_PRECONDITIONED_DIRECTION:
316:     norm_p = rz;
317:     break;

319:   default:
320:     VecDot(p, p, &norm_p);
321:     break;
322:   }

324:   /***************************************************************************/
325:   /* Check for negative curvature.                                           */
326:   /***************************************************************************/

328:   if (kappa <= 0.0) {
329:     /*************************************************************************/
330:     /* In this case, the matrix is indefinite and we have encountered a      */
331:     /* direction of negative curvature.  Because negative curvature occurs   */
332:     /* during the first step, we must follow a direction.                    */
333:     /*************************************************************************/

335:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
336:     PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);

338:     if (cg->radius != 0.0 && norm_p > 0.0) {
339:       /***********************************************************************/
340:       /* Follow direction of negative curvature to the boundary of the       */
341:       /* trust region.                                                       */
342:       /***********************************************************************/

344:       step       = PetscSqrtReal(r2 / norm_p);
345:       cg->norm_d = cg->radius;

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

349:       /***********************************************************************/
350:       /* Update objective function.                                          */
351:       /***********************************************************************/

353:       cg->o_fcn += step * (0.5 * step * kappa - rz);
354:     } else if (cg->radius != 0.0) {
355:       /***********************************************************************/
356:       /* The norm of the preconditioned direction is zero; use the gradient  */
357:       /* step.                                                               */
358:       /***********************************************************************/

360:       if (r2 >= rr) {
361:         alpha      = 1.0;
362:         cg->norm_d = PetscSqrtReal(rr);
363:       } else {
364:         alpha      = PetscSqrtReal(r2 / rr);
365:         cg->norm_d = cg->radius;
366:       }

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

370:       /***********************************************************************/
371:       /* Compute objective function.                                         */
372:       /***********************************************************************/

374:       KSP_MatMult(ksp, Qmat, d, z);
375:       VecAYPX(z, -0.5, ksp->vec_rhs);
376:       VecDot(d, z, &cg->o_fcn);

378:       cg->o_fcn = -cg->o_fcn;
379:       ++ksp->its;
380:     }
381:     return(0);
382:   }

384:   /***************************************************************************/
385:   /* Run the conjugate gradient method until either the problem is solved,   */
386:   /* we encounter the boundary of the trust region, or the conjugate         */
387:   /* gradient method breaks down.                                            */
388:   /***************************************************************************/

390:   while (1) {
391:     /*************************************************************************/
392:     /* Know that kappa is nonzero, because we have not broken down, so we    */
393:     /* can compute the steplength.                                           */
394:     /*************************************************************************/

396:     alpha = rz / kappa;

398:     /*************************************************************************/
399:     /* Compute the steplength and check for intersection with the trust      */
400:     /* region.                                                               */
401:     /*************************************************************************/

403:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
404:     if (cg->radius != 0.0 && norm_dp1 >= r2) {
405:       /***********************************************************************/
406:       /* In this case, the matrix is positive definite as far as we know.    */
407:       /* However, the full step goes beyond the trust region.                */
408:       /***********************************************************************/

410:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
411:       PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", (double)cg->radius);

413:       if (norm_p > 0.0) {
414:         /*********************************************************************/
415:         /* Follow the direction to the boundary of the trust region.         */
416:         /*********************************************************************/

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

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

423:         /*********************************************************************/
424:         /* Update objective function.                                        */
425:         /*********************************************************************/

427:         cg->o_fcn += step * (0.5 * step * kappa - rz);
428:       } else {
429:         /*********************************************************************/
430:         /* The norm of the direction is zero; there is nothing to follow.    */
431:         /*********************************************************************/
432:       }
433:       break;
434:     }

436:     /*************************************************************************/
437:     /* Now we can update the direction and residual.                         */
438:     /*************************************************************************/

440:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
441:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
442:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

444:     switch (cg->dtype) {
445:     case STCG_PRECONDITIONED_DIRECTION:
446:       norm_d = norm_dp1;
447:       break;

449:     default:
450:       VecDot(d, d, &norm_d);
451:       break;
452:     }
453:     cg->norm_d = PetscSqrtReal(norm_d);

455:     /*************************************************************************/
456:     /* Update objective function.                                            */
457:     /*************************************************************************/

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

461:     /*************************************************************************/
462:     /* Check that the preconditioner appears positive semidefinite.          */
463:     /*************************************************************************/

465:     rzm1 = rz;
466:     VecDot(r, z, &rz);             /* rz = r^T z        */
467:     if (rz < 0.0) {
468:       /***********************************************************************/
469:       /* The preconditioner is indefinite.                                   */
470:       /***********************************************************************/

472:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
473:       PetscInfo1(ksp, "KSPSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz);
474:       break;
475:     }

477:     /*************************************************************************/
478:     /* As far as we know, the preconditioner is positive semidefinite.       */
479:     /* Compute the residual and check for convergence.                       */
480:     /*************************************************************************/

482:     switch (ksp->normtype) {
483:     case KSP_NORM_PRECONDITIONED:
484:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
485:       break;

487:     case KSP_NORM_UNPRECONDITIONED:
488:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
489:       break;

491:     case KSP_NORM_NATURAL:
492:       norm_r = PetscSqrtReal(rz);                               /* norm_r = |r|_M    */
493:       break;

495:     default:
496:       norm_r = 0.0;
497:       break;
498:     }

500:     KSPLogResidualHistory(ksp, norm_r);
501:     KSPMonitor(ksp, ksp->its, norm_r);
502:     ksp->rnorm = norm_r;

504:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
505:     if (ksp->reason) {
506:       /***********************************************************************/
507:       /* The method has converged.                                           */
508:       /***********************************************************************/

510:       PetscInfo2(ksp, "KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
511:       break;
512:     }

514:     /*************************************************************************/
515:     /* We have not converged yet.  Check for breakdown.                      */
516:     /*************************************************************************/

518:     beta = rz / rzm1;
519:     if (PetscAbsScalar(beta) <= 0.0) {
520:       /***********************************************************************/
521:       /* Conjugate gradients has broken down.                                */
522:       /***********************************************************************/

524:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
525:       PetscInfo1(ksp, "KSPSolve_STCG: breakdown: beta=%g\n", (double)beta);
526:       break;
527:     }

529:     /*************************************************************************/
530:     /* Check iteration limit.                                                */
531:     /*************************************************************************/

533:     if (ksp->its >= max_cg_its) {
534:       ksp->reason = KSP_DIVERGED_ITS;
535:       PetscInfo1(ksp, "KSPSolve_STCG: iterlim: its=%D\n", ksp->its);
536:       break;
537:     }

539:     /*************************************************************************/
540:     /* Update p and the norms.                                               */
541:     /*************************************************************************/

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

545:     switch (cg->dtype) {
546:     case STCG_PRECONDITIONED_DIRECTION:
547:       dMp    = beta*(dMp + alpha*norm_p);
548:       norm_p = beta*(rzm1 + beta*norm_p);
549:       break;

551:     default:
552:       VecDot(d, p, &dMp);
553:       VecDot(p, p, &norm_p);
554:       break;
555:     }

557:     /*************************************************************************/
558:     /* Compute the new direction and update the iteration.                   */
559:     /*************************************************************************/

561:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
562:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
563:     ++ksp->its;

565:     /*************************************************************************/
566:     /* Check for negative curvature.                                         */
567:     /*************************************************************************/

569:     if (kappa <= 0.0) {
570:       /***********************************************************************/
571:       /* In this case, the matrix is indefinite and we have encountered      */
572:       /* a direction of negative curvature.  Follow the direction to the     */
573:       /* boundary of the trust region.                                       */
574:       /***********************************************************************/

576:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
577:       PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", (double)kappa);

579:       if (cg->radius != 0.0 && norm_p > 0.0) {
580:         /*********************************************************************/
581:         /* Follow direction of negative curvature to boundary.               */
582:         /*********************************************************************/

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

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

589:         /*********************************************************************/
590:         /* Update objective function.                                        */
591:         /*********************************************************************/

593:         cg->o_fcn += step * (0.5 * step * kappa - rz);
594:       } else if (cg->radius != 0.0) {
595:         /*********************************************************************/
596:         /* The norm of the direction is zero; there is nothing to follow.    */
597:         /*********************************************************************/
598:       }
599:       break;
600:     }
601:   }
602:   return(0);
603: #endif
604: }

608: PetscErrorCode KSPSetUp_STCG(KSP ksp)
609: {

612:   /***************************************************************************/
613:   /* Set work vectors needed by conjugate gradient method and allocate       */
614:   /***************************************************************************/

617:   KSPSetWorkVecs(ksp, 3);
618:   return(0);
619: }

623: PetscErrorCode KSPDestroy_STCG(KSP ksp)
624: {

627:   /***************************************************************************/
628:   /* Clear composed functions                                                */
629:   /***************************************************************************/

632:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",NULL);
633:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",NULL);
634:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",NULL);

636:   /***************************************************************************/
637:   /* Destroy KSP object.                                                     */
638:   /***************************************************************************/

640:   KSPDestroyDefault(ksp);
641:   return(0);
642: }

646: static PetscErrorCode  KSPSTCGSetRadius_STCG(KSP ksp, PetscReal radius)
647: {
648:   KSP_STCG *cg = (KSP_STCG*)ksp->data;

651:   cg->radius = radius;
652:   return(0);
653: }

657: static PetscErrorCode  KSPSTCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
658: {
659:   KSP_STCG *cg = (KSP_STCG*)ksp->data;

662:   *norm_d = cg->norm_d;
663:   return(0);
664: }

668: static PetscErrorCode  KSPSTCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
669: {
670:   KSP_STCG *cg = (KSP_STCG*)ksp->data;

673:   *o_fcn = cg->o_fcn;
674:   return(0);
675: }

679: PetscErrorCode KSPSetFromOptions_STCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
680: {
682:   KSP_STCG       *cg = (KSP_STCG*)ksp->data;

685:   PetscOptionsHead(PetscOptionsObject,"KSP STCG options");
686:   PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, NULL);
687:   PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
688:   PetscOptionsTail();
689:   return(0);
690: }

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

697:    Options Database Keys:
698: .      -ksp_stcg_radius <r> - Trust Region Radius

700:    Notes: This is rarely used directly

702:   Use preconditioned conjugate gradient to compute
703:   an approximate minimizer of the quadratic function

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

707:    subject to the trust region constraint

709:             || s || <= delta,

711:    where

713:      delta is the trust region radius,
714:      g is the gradient vector,
715:      H is the Hessian approximation, and
716:      M is the positive definite preconditioner matrix.

718:    KSPConvergedReason may be
719: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
720: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
721: $  other KSP converged/diverged reasons

723:   Notes:
724:   The preconditioner supplied should be symmetric and positive definite.

726:    Level: developer

728: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD(), KSPSTCGGetObjFcn()
729: M*/

733: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
734: {
736:   KSP_STCG       *cg;

739:   PetscNewLog(ksp,&cg);

741:   cg->radius = 0.0;
742:   cg->dtype  = STCG_UNPRECONDITIONED_DIRECTION;

744:   ksp->data = (void*) cg;
745:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
746:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
747:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

749:   /***************************************************************************/
750:   /* Sets the functions that are associated with this data structure         */
751:   /* (in C++ this is the same as defining virtual functions).                */
752:   /***************************************************************************/

754:   ksp->ops->setup          = KSPSetUp_STCG;
755:   ksp->ops->solve          = KSPSolve_STCG;
756:   ksp->ops->destroy        = KSPDestroy_STCG;
757:   ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
758:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
759:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
760:   ksp->ops->view           = 0;

762:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",KSPSTCGSetRadius_STCG);
763:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",KSPSTCGGetNormD_STCG);
764:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",KSPSTCGGetObjFcn_STCG);
765:   return(0);
766: }