Actual source code: stcg.c

petsc-3.5.0 2014-06-30
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:   if (PetscIsInfOrNanScalar(rr)) {
153:     /*************************************************************************/
154:     /* The right-hand side contains not-a-number or an infinite value.       */
155:     /* The gradient step does not work; return a zero value for the step.    */
156:     /*************************************************************************/

158:     ksp->reason = KSP_DIVERGED_NANORINF;
159:     PetscInfo1(ksp, "KSPSolve_STCG: bad right-hand side: rr=%g\n", (double)rr);
160:     return(0);
161:   }

163:   /***************************************************************************/
164:   /* Check the preconditioner for numerical problems and for positive        */
165:   /* definiteness.  The check for not-a-number and infinite values need be   */
166:   /* performed only once.                                                    */
167:   /***************************************************************************/

169:   KSP_PCApply(ksp, r, z);          /* z = inv(M) r      */
170:   VecDot(r, z, &rz);               /* rz = r^T inv(M) r */
171:   if (PetscIsInfOrNanScalar(rz)) {
172:     /*************************************************************************/
173:     /* The preconditioner contains not-a-number or an infinite value.        */
174:     /* Return the gradient direction intersected with the trust region.      */
175:     /*************************************************************************/

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

180:     if (cg->radius != 0) {
181:       if (r2 >= rr) {
182:         alpha      = 1.0;
183:         cg->norm_d = PetscSqrtReal(rr);
184:       } else {
185:         alpha      = PetscSqrtReal(r2 / rr);
186:         cg->norm_d = cg->radius;
187:       }

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

191:       /***********************************************************************/
192:       /* Compute objective function.                                         */
193:       /***********************************************************************/

195:       KSP_MatMult(ksp, Qmat, d, z);
196:       VecAYPX(z, -0.5, ksp->vec_rhs);
197:       VecDot(d, z, &cg->o_fcn);
198:       cg->o_fcn = -cg->o_fcn;
199:       ++ksp->its;
200:     }
201:     return(0);
202:   }

204:   if (rz < 0.0) {
205:     /*************************************************************************/
206:     /* The preconditioner is indefinite.  Because this is the first          */
207:     /* and we do not have a direction yet, we use the gradient step.  Note   */
208:     /* that we cannot use the preconditioned norm when computing the step    */
209:     /* because the matrix is indefinite.                                     */
210:     /*************************************************************************/

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

215:     if (cg->radius != 0.0) {
216:       if (r2 >= rr) {
217:         alpha      = 1.0;
218:         cg->norm_d = PetscSqrtReal(rr);
219:       } else {
220:         alpha      = PetscSqrtReal(r2 / rr);
221:         cg->norm_d = cg->radius;
222:       }

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

226:       /***********************************************************************/
227:       /* Compute objective function.                                         */
228:       /***********************************************************************/

230:       KSP_MatMult(ksp, Qmat, d, z);
231:       VecAYPX(z, -0.5, ksp->vec_rhs);
232:       VecDot(d, z, &cg->o_fcn);
233:       cg->o_fcn = -cg->o_fcn;
234:       ++ksp->its;
235:     }
236:     return(0);
237:   }

239:   /***************************************************************************/
240:   /* As far as we know, the preconditioner is positive semidefinite.         */
241:   /* Compute and log the residual.  Check convergence because this           */
242:   /* initializes things, but do not terminate until at least one conjugate   */
243:   /* gradient iteration has been performed.                                  */
244:   /***************************************************************************/

246:   switch (ksp->normtype) {
247:   case KSP_NORM_PRECONDITIONED:
248:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
249:     break;

251:   case KSP_NORM_UNPRECONDITIONED:
252:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
253:     break;

255:   case KSP_NORM_NATURAL:
256:     norm_r = PetscSqrtReal(rz);                                 /* norm_r = |r|_M    */
257:     break;

259:   default:
260:     norm_r = 0.0;
261:     break;
262:   }

264:   KSPLogResidualHistory(ksp, norm_r);
265:   KSPMonitor(ksp, ksp->its, norm_r);
266:   ksp->rnorm = norm_r;

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

270:   /***************************************************************************/
271:   /* Compute the first direction and update the iteration.                   */
272:   /***************************************************************************/

274:   VecCopy(z, p);                   /* p = z             */
275:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
276:   ++ksp->its;

278:   /***************************************************************************/
279:   /* Check the matrix for numerical problems.                                */
280:   /***************************************************************************/

282:   VecDot(p, z, &kappa);            /* kappa = p^T Q p   */
283:   if (PetscIsInfOrNanScalar(kappa)) {
284:     /*************************************************************************/
285:     /* The matrix produced not-a-number or an infinite value.  In this case, */
286:     /* we must stop and use the gradient direction.  This condition need     */
287:     /* only be checked once.                                                 */
288:     /*************************************************************************/

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

293:     if (cg->radius) {
294:       if (r2 >= rr) {
295:         alpha      = 1.0;
296:         cg->norm_d = PetscSqrtReal(rr);
297:       } else {
298:         alpha      = PetscSqrtReal(r2 / rr);
299:         cg->norm_d = cg->radius;
300:       }

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

304:       /***********************************************************************/
305:       /* Compute objective function.                                         */
306:       /***********************************************************************/

308:       KSP_MatMult(ksp, Qmat, d, z);
309:       VecAYPX(z, -0.5, ksp->vec_rhs);
310:       VecDot(d, z, &cg->o_fcn);
311:       cg->o_fcn = -cg->o_fcn;
312:       ++ksp->its;
313:     }
314:     return(0);
315:   }

317:   /***************************************************************************/
318:   /* Initialize variables for calculating the norm of the direction.         */
319:   /***************************************************************************/

321:   dMp    = 0.0;
322:   norm_d = 0.0;
323:   switch (cg->dtype) {
324:   case STCG_PRECONDITIONED_DIRECTION:
325:     norm_p = rz;
326:     break;

328:   default:
329:     VecDot(p, p, &norm_p);
330:     break;
331:   }

333:   /***************************************************************************/
334:   /* Check for negative curvature.                                           */
335:   /***************************************************************************/

337:   if (kappa <= 0.0) {
338:     /*************************************************************************/
339:     /* In this case, the matrix is indefinite and we have encountered a      */
340:     /* direction of negative curvature.  Because negative curvature occurs   */
341:     /* during the first step, we must follow a direction.                    */
342:     /*************************************************************************/

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

347:     if (cg->radius != 0.0 && norm_p > 0.0) {
348:       /***********************************************************************/
349:       /* Follow direction of negative curvature to the boundary of the       */
350:       /* trust region.                                                       */
351:       /***********************************************************************/

353:       step       = PetscSqrtReal(r2 / norm_p);
354:       cg->norm_d = cg->radius;

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

358:       /***********************************************************************/
359:       /* Update objective function.                                          */
360:       /***********************************************************************/

362:       cg->o_fcn += step * (0.5 * step * kappa - rz);
363:     } else if (cg->radius != 0.0) {
364:       /***********************************************************************/
365:       /* The norm of the preconditioned direction is zero; use the gradient  */
366:       /* step.                                                               */
367:       /***********************************************************************/

369:       if (r2 >= rr) {
370:         alpha      = 1.0;
371:         cg->norm_d = PetscSqrtReal(rr);
372:       } else {
373:         alpha      = PetscSqrtReal(r2 / rr);
374:         cg->norm_d = cg->radius;
375:       }

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

379:       /***********************************************************************/
380:       /* Compute objective function.                                         */
381:       /***********************************************************************/

383:       KSP_MatMult(ksp, Qmat, d, z);
384:       VecAYPX(z, -0.5, ksp->vec_rhs);
385:       VecDot(d, z, &cg->o_fcn);

387:       cg->o_fcn = -cg->o_fcn;
388:       ++ksp->its;
389:     }
390:     return(0);
391:   }

393:   /***************************************************************************/
394:   /* Run the conjugate gradient method until either the problem is solved,   */
395:   /* we encounter the boundary of the trust region, or the conjugate         */
396:   /* gradient method breaks down.                                            */
397:   /***************************************************************************/

399:   while (1) {
400:     /*************************************************************************/
401:     /* Know that kappa is nonzero, because we have not broken down, so we    */
402:     /* can compute the steplength.                                           */
403:     /*************************************************************************/

405:     alpha = rz / kappa;

407:     /*************************************************************************/
408:     /* Compute the steplength and check for intersection with the trust      */
409:     /* region.                                                               */
410:     /*************************************************************************/

412:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
413:     if (cg->radius != 0.0 && norm_dp1 >= r2) {
414:       /***********************************************************************/
415:       /* In this case, the matrix is positive definite as far as we know.    */
416:       /* However, the full step goes beyond the trust region.                */
417:       /***********************************************************************/

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

422:       if (norm_p > 0.0) {
423:         /*********************************************************************/
424:         /* Follow the direction to the boundary of the trust region.         */
425:         /*********************************************************************/

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

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

432:         /*********************************************************************/
433:         /* Update objective function.                                        */
434:         /*********************************************************************/

436:         cg->o_fcn += step * (0.5 * step * kappa - rz);
437:       } else {
438:         /*********************************************************************/
439:         /* The norm of the direction is zero; there is nothing to follow.    */
440:         /*********************************************************************/
441:       }
442:       break;
443:     }

445:     /*************************************************************************/
446:     /* Now we can update the direction and residual.                         */
447:     /*************************************************************************/

449:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
450:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
451:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

453:     switch (cg->dtype) {
454:     case STCG_PRECONDITIONED_DIRECTION:
455:       norm_d = norm_dp1;
456:       break;

458:     default:
459:       VecDot(d, d, &norm_d);
460:       break;
461:     }
462:     cg->norm_d = PetscSqrtReal(norm_d);

464:     /*************************************************************************/
465:     /* Update objective function.                                            */
466:     /*************************************************************************/

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

470:     /*************************************************************************/
471:     /* Check that the preconditioner appears positive semidefinite.          */
472:     /*************************************************************************/

474:     rzm1 = rz;
475:     VecDot(r, z, &rz);             /* rz = r^T z        */
476:     if (rz < 0.0) {
477:       /***********************************************************************/
478:       /* The preconditioner is indefinite.                                   */
479:       /***********************************************************************/

481:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
482:       PetscInfo1(ksp, "KSPSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz);
483:       break;
484:     }

486:     /*************************************************************************/
487:     /* As far as we know, the preconditioner is positive semidefinite.       */
488:     /* Compute the residual and check for convergence.                       */
489:     /*************************************************************************/

491:     switch (ksp->normtype) {
492:     case KSP_NORM_PRECONDITIONED:
493:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
494:       break;

496:     case KSP_NORM_UNPRECONDITIONED:
497:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
498:       break;

500:     case KSP_NORM_NATURAL:
501:       norm_r = PetscSqrtReal(rz);                               /* norm_r = |r|_M    */
502:       break;

504:     default:
505:       norm_r = 0.0;
506:       break;
507:     }

509:     KSPLogResidualHistory(ksp, norm_r);
510:     KSPMonitor(ksp, ksp->its, norm_r);
511:     ksp->rnorm = norm_r;

513:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
514:     if (ksp->reason) {
515:       /***********************************************************************/
516:       /* The method has converged.                                           */
517:       /***********************************************************************/

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

523:     /*************************************************************************/
524:     /* We have not converged yet.  Check for breakdown.                      */
525:     /*************************************************************************/

527:     beta = rz / rzm1;
528:     if (PetscAbsScalar(beta) <= 0.0) {
529:       /***********************************************************************/
530:       /* Conjugate gradients has broken down.                                */
531:       /***********************************************************************/

533:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
534:       PetscInfo1(ksp, "KSPSolve_STCG: breakdown: beta=%g\n", (double)beta);
535:       break;
536:     }

538:     /*************************************************************************/
539:     /* Check iteration limit.                                                */
540:     /*************************************************************************/

542:     if (ksp->its >= max_cg_its) {
543:       ksp->reason = KSP_DIVERGED_ITS;
544:       PetscInfo1(ksp, "KSPSolve_STCG: iterlim: its=%D\n", ksp->its);
545:       break;
546:     }

548:     /*************************************************************************/
549:     /* Update p and the norms.                                               */
550:     /*************************************************************************/

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

554:     switch (cg->dtype) {
555:     case STCG_PRECONDITIONED_DIRECTION:
556:       dMp    = beta*(dMp + alpha*norm_p);
557:       norm_p = beta*(rzm1 + beta*norm_p);
558:       break;

560:     default:
561:       VecDot(d, p, &dMp);
562:       VecDot(p, p, &norm_p);
563:       break;
564:     }

566:     /*************************************************************************/
567:     /* Compute the new direction and update the iteration.                   */
568:     /*************************************************************************/

570:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
571:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
572:     ++ksp->its;

574:     /*************************************************************************/
575:     /* Check for negative curvature.                                         */
576:     /*************************************************************************/

578:     if (kappa <= 0.0) {
579:       /***********************************************************************/
580:       /* In this case, the matrix is indefinite and we have encountered      */
581:       /* a direction of negative curvature.  Follow the direction to the     */
582:       /* boundary of the trust region.                                       */
583:       /***********************************************************************/

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

588:       if (cg->radius != 0.0 && norm_p > 0.0) {
589:         /*********************************************************************/
590:         /* Follow direction of negative curvature to boundary.               */
591:         /*********************************************************************/

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

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

598:         /*********************************************************************/
599:         /* Update objective function.                                        */
600:         /*********************************************************************/

602:         cg->o_fcn += step * (0.5 * step * kappa - rz);
603:       } else if (cg->radius != 0.0) {
604:         /*********************************************************************/
605:         /* The norm of the direction is zero; there is nothing to follow.    */
606:         /*********************************************************************/
607:       }
608:       break;
609:     }
610:   }
611:   return(0);
612: #endif
613: }

617: PetscErrorCode KSPSetUp_STCG(KSP ksp)
618: {

621:   /***************************************************************************/
622:   /* Set work vectors needed by conjugate gradient method and allocate       */
623:   /***************************************************************************/

626:   KSPSetWorkVecs(ksp, 3);
627:   return(0);
628: }

632: PetscErrorCode KSPDestroy_STCG(KSP ksp)
633: {

636:   /***************************************************************************/
637:   /* Clear composed functions                                                */
638:   /***************************************************************************/

641:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",NULL);
642:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",NULL);
643:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",NULL);

645:   /***************************************************************************/
646:   /* Destroy KSP object.                                                     */
647:   /***************************************************************************/

649:   KSPDestroyDefault(ksp);
650:   return(0);
651: }

655: static PetscErrorCode  KSPSTCGSetRadius_STCG(KSP ksp, PetscReal radius)
656: {
657:   KSP_STCG *cg = (KSP_STCG*)ksp->data;

660:   cg->radius = radius;
661:   return(0);
662: }

666: static PetscErrorCode  KSPSTCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
667: {
668:   KSP_STCG *cg = (KSP_STCG*)ksp->data;

671:   *norm_d = cg->norm_d;
672:   return(0);
673: }

677: static PetscErrorCode  KSPSTCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
678: {
679:   KSP_STCG *cg = (KSP_STCG*)ksp->data;

682:   *o_fcn = cg->o_fcn;
683:   return(0);
684: }

688: PetscErrorCode KSPSetFromOptions_STCG(KSP ksp)
689: {
691:   KSP_STCG       *cg = (KSP_STCG*)ksp->data;

694:   PetscOptionsHead("KSP STCG options");
695:   PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, NULL);
696:   PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
697:   PetscOptionsTail();
698:   return(0);
699: }

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

706:    Options Database Keys:
707: .      -ksp_stcg_radius <r> - Trust Region Radius

709:    Notes: This is rarely used directly

711:   Use preconditioned conjugate gradient to compute
712:   an approximate minimizer of the quadratic function

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

716:    subject to the trust region constraint

718:             || s || <= delta,

720:    where

722:      delta is the trust region radius,
723:      g is the gradient vector,
724:      H is the Hessian approximation, and
725:      M is the positive definite preconditioner matrix.

727:    KSPConvergedReason may be
728: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
729: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
730: $  other KSP converged/diverged reasons

732:   Notes:
733:   The preconditioner supplied should be symmetric and positive definite.

735:    Level: developer

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

742: PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
743: {
745:   KSP_STCG       *cg;

748:   PetscNewLog(ksp,&cg);

750:   cg->radius = 0.0;
751:   cg->dtype  = STCG_UNPRECONDITIONED_DIRECTION;

753:   ksp->data = (void*) cg;
754:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
755:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
756:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);

758:   /***************************************************************************/
759:   /* Sets the functions that are associated with this data structure         */
760:   /* (in C++ this is the same as defining virtual functions).                */
761:   /***************************************************************************/

763:   ksp->ops->setup          = KSPSetUp_STCG;
764:   ksp->ops->solve          = KSPSolve_STCG;
765:   ksp->ops->destroy        = KSPDestroy_STCG;
766:   ksp->ops->setfromoptions = KSPSetFromOptions_STCG;
767:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
768:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
769:   ksp->ops->view           = 0;

771:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGSetRadius_C",KSPSTCGSetRadius_STCG);
772:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetNormD_C",KSPSTCGGetNormD_STCG);
773:   PetscObjectComposeFunction((PetscObject)ksp,"KSPSTCGGetObjFcn_C",KSPSTCGGetObjFcn_STCG);
774:   return(0);
775: }