Actual source code: stcg.c

petsc-3.4.4 2014-03-13
  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: }