Actual source code: nash.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  2: #include <petsc-private/kspimpl.h>             /*I "petscksp.h" I*/
  3: #include <../src/ksp/ksp/impls/cg/nash/nashimpl.h>

  5: #define NASH_PRECONDITIONED_DIRECTION   0
  6: #define NASH_UNPRECONDITIONED_DIRECTION 1
  7: #define NASH_DIRECTION_TYPES            2

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

 13: /*@
 14:     KSPNASHSetRadius - 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_nash_radius <r>

 25:     Level: advanced

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

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

 43: /*@
 44:     KSPNASHGetNormD - 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, NASH, get, norm direction
 55: @*/
 56: PetscErrorCode  KSPNASHGetNormD(KSP ksp, PetscReal *norm_d)
 57: {

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

 68: /*@
 69:     KSPNASHGetObjFcn - 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, NASH, get, objective function
 80: @*/
 81: PetscErrorCode  KSPNASHGetObjFcn(KSP ksp, PetscReal *o_fcn)
 82: {

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


 94: PetscErrorCode KSPSolve_NASH(KSP ksp)
 95: {
 96: #if defined(PETSC_USE_COMPLEX)
 97:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "NASH is not available for complex systems");
 98: #else
 99:   KSP_NASH       *cg = (KSP_NASH*)ksp->data;
101:   Mat            Qmat, Mmat;
102:   Vec            r, z, p, d;
103:   PC             pc;

105:   PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
106:   PetscReal alpha, beta, kappa, rz, rzm1;
107:   PetscReal rr, r2, step;

109:   PetscInt max_cg_its;

111:   PetscBool diagonalscale;

114:   /***************************************************************************/
115:   /* Check the arguments and parameters.                                     */
116:   /***************************************************************************/

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

122:   /***************************************************************************/
123:   /* Get the workspace vectors and initialize variables                      */
124:   /***************************************************************************/

126:   r2 = cg->radius * cg->radius;
127:   r  = ksp->work[0];
128:   z  = ksp->work[1];
129:   p  = ksp->work[2];
130:   d  = ksp->vec_sol;
131:   pc = ksp->pc;

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

135:   VecGetSize(d, &max_cg_its);
136:   max_cg_its = PetscMin(max_cg_its, ksp->max_it);
137:   ksp->its   = 0;

139:   /***************************************************************************/
140:   /* Initialize objective function and direction.                            */
141:   /***************************************************************************/

143:   cg->o_fcn = 0.0;

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

148:   /***************************************************************************/
149:   /* Begin the conjugate gradient method.  Check the right-hand side for     */
150:   /* numerical problems.  The check for not-a-number and infinite values     */
151:   /* need be performed only once.                                            */
152:   /***************************************************************************/

154:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
155:   VecDot(r, r, &rr);               /* rr = r^T r        */
156:   if (PetscIsInfOrNanScalar(rr)) {
157:     /*************************************************************************/
158:     /* The right-hand side contains not-a-number or an infinite value.       */
159:     /* The gradient step does not work; return a zero value for the step.    */
160:     /*************************************************************************/

162:     ksp->reason = KSP_DIVERGED_NANORINF;
163:     PetscInfo1(ksp, "KSPSolve_NASH: bad right-hand side: rr=%g\n", (double)rr);
164:     return(0);
165:   }

167:   /***************************************************************************/
168:   /* Check the preconditioner for numerical problems and for positive        */
169:   /* definiteness.  The check for not-a-number and infinite values need be   */
170:   /* performed only once.                                                    */
171:   /***************************************************************************/

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

181:     ksp->reason = KSP_DIVERGED_NANORINF;
182:     PetscInfo1(ksp, "KSPSolve_NASH: bad preconditioner: rz=%g\n", (double)rz);

184:     if (cg->radius) {
185:       if (r2 >= rr) {
186:         alpha      = 1.0;
187:         cg->norm_d = PetscSqrtReal(rr);
188:       } else {
189:         alpha      = PetscSqrtReal(r2 / rr);
190:         cg->norm_d = cg->radius;
191:       }

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

195:       /***********************************************************************/
196:       /* Compute objective function.                                         */
197:       /***********************************************************************/

199:       KSP_MatMult(ksp, Qmat, d, z);
200:       VecAYPX(z, -0.5, ksp->vec_rhs);
201:       VecDot(d, z, &cg->o_fcn);
202:       cg->o_fcn = -cg->o_fcn;
203:       ++ksp->its;
204:     }
205:     return(0);
206:   }

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

216:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
217:     PetscInfo1(ksp, "KSPSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz);

219:     if (cg->radius) {
220:       if (r2 >= rr) {
221:         alpha      = 1.0;
222:         cg->norm_d = PetscSqrtReal(rr);
223:       } else {
224:         alpha      = PetscSqrtReal(r2 / rr);
225:         cg->norm_d = cg->radius;
226:       }

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

230:       /***********************************************************************/
231:       /* Compute objective function.                                         */
232:       /***********************************************************************/

234:       KSP_MatMult(ksp, Qmat, d, z);
235:       VecAYPX(z, -0.5, ksp->vec_rhs);
236:       VecDot(d, z, &cg->o_fcn);
237:       cg->o_fcn = -cg->o_fcn;
238:       ++ksp->its;
239:     }
240:     return(0);
241:   }

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

250:   switch (ksp->normtype) {
251:   case KSP_NORM_PRECONDITIONED:
252:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
253:     break;

255:   case KSP_NORM_UNPRECONDITIONED:
256:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
257:     break;

259:   case KSP_NORM_NATURAL:
260:     norm_r = PetscSqrtReal(rz);                                 /* norm_r = |r|_M    */
261:     break;

263:   default:
264:     norm_r = 0.0;
265:     break;
266:   }

268:   KSPLogResidualHistory(ksp, norm_r);
269:   KSPMonitor(ksp, ksp->its, norm_r);
270:   ksp->rnorm = norm_r;

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

274:   /***************************************************************************/
275:   /* Compute the first direction and update the iteration.                   */
276:   /***************************************************************************/

278:   VecCopy(z, p);                   /* p = z             */
279:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
280:   ++ksp->its;

282:   /***************************************************************************/
283:   /* Check the matrix for numerical problems.                                */
284:   /***************************************************************************/

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

294:     ksp->reason = KSP_DIVERGED_NANORINF;
295:     PetscInfo1(ksp, "KSPSolve_NASH: bad matrix: kappa=%g\n", (double)kappa);

297:     if (cg->radius) {
298:       if (r2 >= rr) {
299:         alpha      = 1.0;
300:         cg->norm_d = PetscSqrtReal(rr);
301:       } else {
302:         alpha      = PetscSqrtReal(r2 / rr);
303:         cg->norm_d = cg->radius;
304:       }

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

308:       /***********************************************************************/
309:       /* Compute objective function.                                         */
310:       /***********************************************************************/

312:       KSP_MatMult(ksp, Qmat, d, z);
313:       VecAYPX(z, -0.5, ksp->vec_rhs);
314:       VecDot(d, z, &cg->o_fcn);
315:       cg->o_fcn = -cg->o_fcn;
316:       ++ksp->its;
317:     }
318:     return(0);
319:   }

321:   /***************************************************************************/
322:   /* Initialize variables for calculating the norm of the direction.         */
323:   /***************************************************************************/

325:   dMp    = 0.0;
326:   norm_d = 0.0;
327:   switch (cg->dtype) {
328:   case NASH_PRECONDITIONED_DIRECTION:
329:     norm_p = rz;
330:     break;

332:   default:
333:     VecDot(p, p, &norm_p);
334:     break;
335:   }

337:   /***************************************************************************/
338:   /* Check for negative curvature.                                           */
339:   /***************************************************************************/

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

348:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
349:     PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);

351:     if (cg->radius && norm_p > 0.0) {
352:       /***********************************************************************/
353:       /* Follow direction of negative curvature to the boundary of the       */
354:       /* trust region.                                                       */
355:       /***********************************************************************/

357:       step       = PetscSqrtReal(r2 / norm_p);
358:       cg->norm_d = cg->radius;

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

362:       /***********************************************************************/
363:       /* Update objective function.                                          */
364:       /***********************************************************************/

366:       cg->o_fcn += step * (0.5 * step * kappa - rz);
367:     } else if (cg->radius) {
368:       /***********************************************************************/
369:       /* The norm of the preconditioned direction is zero; use the gradient  */
370:       /* step.                                                               */
371:       /***********************************************************************/

373:       if (r2 >= rr) {
374:         alpha      = 1.0;
375:         cg->norm_d = PetscSqrtReal(rr);
376:       } else {
377:         alpha      = PetscSqrtReal(r2 / rr);
378:         cg->norm_d = cg->radius;
379:       }

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

383:       /***********************************************************************/
384:       /* Compute objective function.                                         */
385:       /***********************************************************************/

387:       KSP_MatMult(ksp, Qmat, d, z);
388:       VecAYPX(z, -0.5, ksp->vec_rhs);
389:       VecDot(d, z, &cg->o_fcn);
390:       cg->o_fcn = -cg->o_fcn;
391:       ++ksp->its;
392:     }
393:     return(0);
394:   }

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

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

408:     alpha = rz / kappa;

410:     /*************************************************************************/
411:     /* Compute the steplength and check for intersection with the trust      */
412:     /* region.                                                               */
413:     /*************************************************************************/

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

422:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
423:       PetscInfo1(ksp, "KSPSolve_NASH: constrained step: radius=%g\n", (double)cg->radius);

425:       if (norm_p > 0.0) {
426:         /*********************************************************************/
427:         /* Follow the direction to the boundary of the trust region.         */
428:         /*********************************************************************/

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

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

435:         /*********************************************************************/
436:         /* Update objective function.                                        */
437:         /*********************************************************************/

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

448:     /*************************************************************************/
449:     /* Now we can update the direction and residual.                         */
450:     /*************************************************************************/

452:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
453:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
454:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

456:     switch (cg->dtype) {
457:     case NASH_PRECONDITIONED_DIRECTION:
458:       norm_d = norm_dp1;
459:       break;

461:     default:
462:       VecDot(d, d, &norm_d);
463:       break;
464:     }
465:     cg->norm_d = PetscSqrtReal(norm_d);

467:     /*************************************************************************/
468:     /* Update objective function.                                            */
469:     /*************************************************************************/

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

473:     /*************************************************************************/
474:     /* Check that the preconditioner appears positive semidefinite.          */
475:     /*************************************************************************/

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

484:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
485:       PetscInfo1(ksp, "KSPSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz);
486:       break;
487:     }

489:     /*************************************************************************/
490:     /* As far as we know, the preconditioner is positive semidefinite.       */
491:     /* Compute the residual and check for convergence.                       */
492:     /*************************************************************************/

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

499:     case KSP_NORM_UNPRECONDITIONED:
500:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
501:       break;

503:     case KSP_NORM_NATURAL:
504:       norm_r = PetscSqrtReal(rz);                               /* norm_r = |r|_M    */
505:       break;

507:     default:
508:       norm_r = 0.;
509:       break;
510:     }

512:     KSPLogResidualHistory(ksp, norm_r);
513:     KSPMonitor(ksp, ksp->its, norm_r);
514:     ksp->rnorm = norm_r;

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

522:       PetscInfo2(ksp, "KSPSolve_NASH: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
523:       break;
524:     }

526:     /*************************************************************************/
527:     /* We have not converged yet.  Check for breakdown.                      */
528:     /*************************************************************************/

530:     beta = rz / rzm1;
531:     if (PetscAbsReal(beta) <= 0.0) {
532:       /***********************************************************************/
533:       /* Conjugate gradients has broken down.                                */
534:       /***********************************************************************/

536:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
537:       PetscInfo1(ksp, "KSPSolve_NASH: breakdown: beta=%g\n", (double)beta);
538:       break;
539:     }

541:     /*************************************************************************/
542:     /* Check iteration limit.                                                */
543:     /*************************************************************************/

545:     if (ksp->its >= max_cg_its) {
546:       ksp->reason = KSP_DIVERGED_ITS;
547:       PetscInfo1(ksp, "KSPSolve_NASH: iterlim: its=%D\n", ksp->its);
548:       break;
549:     }

551:     /*************************************************************************/
552:     /* Update p and the norms.                                               */
553:     /*************************************************************************/

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

557:     switch (cg->dtype) {
558:     case NASH_PRECONDITIONED_DIRECTION:
559:       dMp    = beta*(dMp + alpha*norm_p);
560:       norm_p = beta*(rzm1 + beta*norm_p);
561:       break;

563:     default:
564:       VecDot(d, p, &dMp);
565:       VecDot(p, p, &norm_p);
566:       break;
567:     }

569:     /*************************************************************************/
570:     /* Compute the new direction and update the iteration.                   */
571:     /*************************************************************************/

573:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
574:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
575:     ++ksp->its;

577:     /*************************************************************************/
578:     /* Check for negative curvature.                                         */
579:     /*************************************************************************/

581:     if (kappa <= 0.0) {
582:       /***********************************************************************/
583:       /* In this case, the matrix is indefinite and we have encountered      */
584:       /* a direction of negative curvature.  Stop at the base.               */
585:       /***********************************************************************/

587:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
588:       PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
589:       break;
590:     }
591:   }
592:   return(0);
593: #endif
594: }

598: PetscErrorCode KSPSetUp_NASH(KSP ksp)
599: {

603:   /***************************************************************************/
604:   /* Set work vectors needed by conjugate gradient method and allocate       */
605:   /***************************************************************************/

607:   KSPSetWorkVecs(ksp, 3);
608:   return(0);
609: }

613: PetscErrorCode KSPDestroy_NASH(KSP ksp)
614: {

618:   /***************************************************************************/
619:   /* Clear composed functions                                                */
620:   /***************************************************************************/

622:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",NULL);
623:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",NULL);
624:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",NULL);

626:   /***************************************************************************/
627:   /* Destroy KSP object.                                                     */
628:   /***************************************************************************/

630:   KSPDestroyDefault(ksp);
631:   return(0);
632: }

636: static PetscErrorCode  KSPNASHSetRadius_NASH(KSP ksp, PetscReal radius)
637: {
638:   KSP_NASH *cg = (KSP_NASH*)ksp->data;

641:   cg->radius = radius;
642:   return(0);
643: }

647: static PetscErrorCode  KSPNASHGetNormD_NASH(KSP ksp, PetscReal *norm_d)
648: {
649:   KSP_NASH *cg = (KSP_NASH*)ksp->data;

652:   *norm_d = cg->norm_d;
653:   return(0);
654: }

658: static PetscErrorCode  KSPNASHGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
659: {
660:   KSP_NASH *cg = (KSP_NASH*)ksp->data;

663:   *o_fcn = cg->o_fcn;
664:   return(0);
665: }

669: PetscErrorCode KSPSetFromOptions_NASH(KSP ksp)
670: {
672:   KSP_NASH       *cg = (KSP_NASH*)ksp->data;

675:   PetscOptionsHead("KSP NASH options");

677:   PetscOptionsReal("-ksp_nash_radius", "Trust Region Radius", "KSPNASHSetRadius", cg->radius, &cg->radius, NULL);

679:   PetscOptionsEList("-ksp_nash_dtype", "Norm used for direction", "", DType_Table, NASH_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);

681:   PetscOptionsTail();
682:   return(0);
683: }

685: /*MC
686:      KSPNASH -   Code to run conjugate gradient method subject to a constraint
687:          on the solution norm. This is used in Trust Region methods for
688:          nonlinear equations, SNESNEWTONTR

690:    Options Database Keys:
691: .      -ksp_nash_radius <r> - Trust Region Radius

693:    Notes: This is rarely used directly

695:    Level: developer

697:   Use preconditioned conjugate gradient to compute
698:   an approximate minimizer of the quadratic function

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

702:    subject to the trust region constraint

704:             || s || <= delta,

706:    where

708:      delta is the trust region radius,
709:      g is the gradient vector,
710:      H is the Hessian approximation, and
711:      M is the positive definite preconditioner matrix.

713:    KSPConvergedReason may be
714: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
715: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
716: $  other KSP converged/diverged reasons

718:   Notes:
719:   The preconditioner supplied should be symmetric and positive definite.

721: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPNASHSetRadius(), KSPNASHGetNormD(), KSPNASHGetObjFcn()
722: M*/

726: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
727: {
729:   KSP_NASH       *cg;

732:   PetscNewLog(ksp,&cg);
733:   cg->radius = 0.0;
734:   cg->dtype  = NASH_UNPRECONDITIONED_DIRECTION;

736:   ksp->data = (void*) cg;
737:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
738:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
739:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);

741:   /***************************************************************************/
742:   /* Sets the functions that are associated with this data structure         */
743:   /* (in C++ this is the same as defining virtual functions).                */
744:   /***************************************************************************/

746:   ksp->ops->setup          = KSPSetUp_NASH;
747:   ksp->ops->solve          = KSPSolve_NASH;
748:   ksp->ops->destroy        = KSPDestroy_NASH;
749:   ksp->ops->setfromoptions = KSPSetFromOptions_NASH;
750:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
751:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
752:   ksp->ops->view           = 0;

754:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",KSPNASHSetRadius_NASH);
755:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",KSPNASHGetNormD_NASH);
756:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",KSPNASHGetObjFcn_NASH);
757:   return(0);
758: }