Actual source code: nash.c

petsc-3.3-p7 2013-05-11
  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(((PetscObject)ksp)->comm,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: #ifdef PETSC_USE_COMPLEX
 97:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP, "NASH is not available for complex systems");
 98: #else
 99:   KSP_NASH       *cg = (KSP_NASH *)ksp->data;
101:   MatStructure   pflag;
102:   Mat            Qmat, Mmat;
103:   Vec            r, z, p, d;
104:   PC             pc;

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

110:   PetscInt       max_cg_its;

112:   PetscBool      diagonalscale;

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

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

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

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

134:   PCGetOperators(pc, &Qmat, &Mmat, &pflag);

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

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

144:   cg->o_fcn = 0.0;

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

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

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

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

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

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

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

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

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

197:       /***********************************************************************/
198:       /* Compute objective function.                                         */
199:       /***********************************************************************/

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

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

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

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

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

233:       /***********************************************************************/
234:       /* Compute objective function.                                         */
235:       /***********************************************************************/

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

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

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

258:   case KSP_NORM_UNPRECONDITIONED:
259:     norm_r = PetscSqrtReal(rr);                                        /* norm_r = |r|      */
260:     break;

262:   case KSP_NORM_NATURAL:
263:     norm_r = PetscSqrtReal(rz);                                        /* norm_r = |r|_M    */
264:     break;

266:   default:
267:     norm_r = 0.0;
268:     break;
269:   }

271:   KSPLogResidualHistory(ksp, norm_r);
272:   KSPMonitor(ksp, ksp->its, norm_r);
273:   ksp->rnorm = norm_r;

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

277:   /***************************************************************************/
278:   /* Compute the first direction and update the iteration.                   */
279:   /***************************************************************************/

281:   VecCopy(z, p);                        /* p = z             */
282:   KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
283:   ++ksp->its;

285:   /***************************************************************************/
286:   /* Check the matrix for numerical problems.                                */
287:   /***************************************************************************/

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

297:     ksp->reason = KSP_DIVERGED_NAN;
298:     PetscInfo1(ksp, "KSPSolve_NASH: bad matrix: kappa=%g\n", kappa);

300:     if (cg->radius) {
301:       if (r2 >= rr) {
302:         alpha = 1.0;
303:         cg->norm_d = PetscSqrtReal(rr);
304:       }
305:       else {
306:         alpha = PetscSqrtReal(r2 / rr);
307:         cg->norm_d = cg->radius;
308:       }

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

312:       /***********************************************************************/
313:       /* Compute objective function.                                         */
314:       /***********************************************************************/

316:       KSP_MatMult(ksp, Qmat, d, z);
317:       VecAYPX(z, -0.5, ksp->vec_rhs);
318:       VecDot(d, z, &cg->o_fcn);
319:       cg->o_fcn = -cg->o_fcn;
320:       ++ksp->its;
321:     }
322:     return(0);
323:   }

325:   /***************************************************************************/
326:   /* Initialize variables for calculating the norm of the direction.         */
327:   /***************************************************************************/

329:   dMp = 0.0;
330:   norm_d = 0.0;
331:   switch(cg->dtype) {
332:   case NASH_PRECONDITIONED_DIRECTION:
333:     norm_p = rz;
334:     break;

336:   default:
337:     VecDot(p, p, &norm_p);
338:     break;
339:   }

341:   /***************************************************************************/
342:   /* Check for negative curvature.                                           */
343:   /***************************************************************************/

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

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

355:     if (cg->radius && norm_p > 0.0) {
356:       /***********************************************************************/
357:       /* Follow direction of negative curvature to the boundary of the       */
358:       /* trust region.                                                       */
359:       /***********************************************************************/

361:       step = PetscSqrtReal(r2 / norm_p);
362:       cg->norm_d = cg->radius;

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

366:       /***********************************************************************/
367:       /* Update objective function.                                          */
368:       /***********************************************************************/

370:       cg->o_fcn += step * (0.5 * step * kappa - rz);
371:     }
372:     else if (cg->radius) {
373:       /***********************************************************************/
374:       /* The norm of the preconditioned direction is zero; use the gradient  */
375:       /* step.                                                               */
376:       /***********************************************************************/

378:       if (r2 >= rr) {
379:         alpha = 1.0;
380:         cg->norm_d = PetscSqrtReal(rr);
381:       }
382:       else {
383:         alpha = PetscSqrtReal(r2 / rr);
384:         cg->norm_d = cg->radius;
385:       }

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

389:       /***********************************************************************/
390:       /* Compute objective function.                                         */
391:       /***********************************************************************/

393:       KSP_MatMult(ksp, Qmat, d, z);
394:       VecAYPX(z, -0.5, ksp->vec_rhs);
395:       VecDot(d, z, &cg->o_fcn);
396:       cg->o_fcn = -cg->o_fcn;
397:       ++ksp->its;
398:     }
399:     return(0);
400:   }

402:   /***************************************************************************/
403:   /* Run the conjugate gradient method until either the problem is solved,   */
404:   /* we encounter the boundary of the trust region, or the conjugate         */
405:   /* gradient method breaks down.                                            */
406:   /***************************************************************************/

408:   while(1) {
409:     /*************************************************************************/
410:     /* Know that kappa is nonzero, because we have not broken down, so we    */
411:     /* can compute the steplength.                                           */
412:     /*************************************************************************/

414:     alpha = rz / kappa;

416:     /*************************************************************************/
417:     /* Compute the steplength and check for intersection with the trust      */
418:     /* region.                                                               */
419:     /*************************************************************************/

421:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
422:     if (cg->radius && norm_dp1 >= r2) {
423:       /***********************************************************************/
424:       /* In this case, the matrix is positive definite as far as we know.    */
425:       /* However, the full step goes beyond the trust region.                */
426:       /***********************************************************************/

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

431:       if (norm_p > 0.0) {
432:         /*********************************************************************/
433:         /* Follow the direction to the boundary of the trust region.         */
434:         /*********************************************************************/

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

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

441:         /*********************************************************************/
442:         /* Update objective function.                                        */
443:         /*********************************************************************/

445:         cg->o_fcn += step * (0.5 * step * kappa - rz);
446:       }
447:       else {
448:         /*********************************************************************/
449:         /* The norm of the direction is zero; there is nothing to follow.    */
450:         /*********************************************************************/
451:       }
452:       break;
453:     }

455:     /*************************************************************************/
456:     /* Now we can update the direction and residual.                         */
457:     /*************************************************************************/

459:     VecAXPY(d, alpha, p);                /* d = d + alpha p   */
460:     VecAXPY(r, -alpha, z);                        /* r = r - alpha Q p */
461:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

463:     switch(cg->dtype) {
464:     case NASH_PRECONDITIONED_DIRECTION:
465:       norm_d = norm_dp1;
466:       break;

468:     default:
469:       VecDot(d, d, &norm_d);
470:       break;
471:     }
472:     cg->norm_d = PetscSqrtReal(norm_d);

474:     /*************************************************************************/
475:     /* Update objective function.                                            */
476:     /*************************************************************************/

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

480:     /*************************************************************************/
481:     /* Check that the preconditioner appears positive semidefinite.          */
482:     /*************************************************************************/

484:     rzm1 = rz;
485:     VecDot(r, z, &rz);                /* rz = r^T z        */
486:     if (rz < 0.0) {
487:       /***********************************************************************/
488:       /* The preconditioner is indefinite.                                   */
489:       /***********************************************************************/

491:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
492:       PetscInfo1(ksp, "KSPSolve_NASH: cg indefinite preconditioner: rz=%g\n", rz);
493:       break;
494:     }

496:     /*************************************************************************/
497:     /* As far as we know, the preconditioner is positive semidefinite.       */
498:     /* Compute the residual and check for convergence.                       */
499:     /*************************************************************************/

501:     switch(ksp->normtype) {
502:     case KSP_NORM_PRECONDITIONED:
503:       VecNorm(z, NORM_2, &norm_r);/* norm_r = |z|      */
504:       break;

506:     case KSP_NORM_UNPRECONDITIONED:
507:       VecNorm(r, NORM_2, &norm_r);/* norm_r = |r|      */
508:       break;

510:     case KSP_NORM_NATURAL:
511:       norm_r = PetscSqrtReal(rz);                                /* norm_r = |r|_M    */
512:       break;

514:     default:
515:       norm_r = 0.;
516:       break;
517:     }

519:     KSPLogResidualHistory(ksp, norm_r);
520:     KSPMonitor(ksp, ksp->its, norm_r);
521:     ksp->rnorm = norm_r;
522: 
523:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
524:     if (ksp->reason) {
525:       /***********************************************************************/
526:       /* The method has converged.                                           */
527:       /***********************************************************************/

529:       PetscInfo2(ksp, "KSPSolve_NASH: truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
530:       break;
531:     }

533:     /*************************************************************************/
534:     /* We have not converged yet.  Check for breakdown.                      */
535:     /*************************************************************************/

537:     beta = rz / rzm1;
538:     if (fabs(beta) <= 0.0) {
539:       /***********************************************************************/
540:       /* Conjugate gradients has broken down.                                */
541:       /***********************************************************************/

543:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
544:       PetscInfo1(ksp, "KSPSolve_NASH: breakdown: beta=%g\n", beta);
545:       break;
546:     }

548:     /*************************************************************************/
549:     /* Check iteration limit.                                                */
550:     /*************************************************************************/

552:     if (ksp->its >= max_cg_its) {
553:       ksp->reason = KSP_DIVERGED_ITS;
554:       PetscInfo1(ksp, "KSPSolve_NASH: iterlim: its=%d\n", ksp->its);
555:       break;
556:     }

558:     /*************************************************************************/
559:     /* Update p and the norms.                                               */
560:     /*************************************************************************/

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

564:     switch(cg->dtype) {
565:     case NASH_PRECONDITIONED_DIRECTION:
566:       dMp = beta*(dMp + alpha*norm_p);
567:       norm_p = beta*(rzm1 + beta*norm_p);
568:       break;

570:     default:
571:       VecDot(d, p, &dMp);
572:       VecDot(p, p, &norm_p);
573:       break;
574:     }

576:     /*************************************************************************/
577:     /* Compute the new direction and update the iteration.                   */
578:     /*************************************************************************/

580:     KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
581:     VecDot(p, z, &kappa);                /* kappa = p^T Q p   */
582:     ++ksp->its;

584:     /*************************************************************************/
585:     /* Check for negative curvature.                                         */
586:     /*************************************************************************/

588:     if (kappa <= 0.0) {
589:       /***********************************************************************/
590:       /* In this case, the matrix is indefinite and we have encountered      */
591:       /* a direction of negative curvature.  Stop at the base.               */
592:       /***********************************************************************/

594:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
595:       PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", kappa);
596:       break;
597:     }
598:   }

600:   return(0);
601: #endif
602: }

606: PetscErrorCode KSPSetUp_NASH(KSP ksp)
607: {

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

615:   KSPDefaultGetWork(ksp, 3);
616:   return(0);
617: }

621: PetscErrorCode KSPDestroy_NASH(KSP ksp)
622: {

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

630:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPNASHSetRadius_C","",PETSC_NULL);
631:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPNASHGetNormD_C","",PETSC_NULL);
632:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPNASHGetObjFcn_C","",PETSC_NULL);

634:   /***************************************************************************/
635:   /* Destroy KSP object.                                                     */
636:   /***************************************************************************/

638:   KSPDefaultDestroy(ksp);
639:   return(0);
640: }

642: EXTERN_C_BEGIN
645: PetscErrorCode  KSPNASHSetRadius_NASH(KSP ksp, PetscReal radius)
646: {
647:   KSP_NASH *cg = (KSP_NASH *)ksp->data;

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

656: PetscErrorCode  KSPNASHGetNormD_NASH(KSP ksp, PetscReal *norm_d)
657: {
658:   KSP_NASH *cg = (KSP_NASH *)ksp->data;

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

667: PetscErrorCode  KSPNASHGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn){
668:   KSP_NASH *cg = (KSP_NASH *)ksp->data;

671:   *o_fcn = cg->o_fcn;
672:   return(0);
673: }
674: EXTERN_C_END

678: PetscErrorCode KSPSetFromOptions_NASH(KSP ksp)
679: {
681:   KSP_NASH       *cg = (KSP_NASH *)ksp->data;

684:   PetscOptionsHead("KSP NASH options");

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

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

690:   PetscOptionsTail();
691:   return(0);
692: }

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

699:    Options Database Keys:
700: .      -ksp_nash_radius <r> - Trust Region Radius

702:    Notes: This is rarely used directly

704:    Level: developer

706:   Use preconditioned conjugate gradient to compute
707:   an approximate minimizer of the quadratic function

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

711:    subject to the trust region constraint

713:             || s || <= delta,

715:    where

717:      delta is the trust region radius,
718:      g is the gradient vector, 
719:      H is the Hessian approximation, and
720:      M is the positive definite preconditioner matrix.

722:    KSPConvergedReason may be
723: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
724: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
725: $  other KSP converged/diverged reasons

727:   Notes:
728:   The preconditioner supplied should be symmetric and positive definite.

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

733: EXTERN_C_BEGIN
736: PetscErrorCode  KSPCreate_NASH(KSP ksp)
737: {
739:   KSP_NASH       *cg;

742:   PetscNewLog(ksp,KSP_NASH, &cg);
743:   cg->radius = 0.0;
744:   cg->dtype = NASH_UNPRECONDITIONED_DIRECTION;

746:   ksp->data = (void *) cg;
747:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
748:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
749:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);

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

756:   ksp->ops->setup                = KSPSetUp_NASH;
757:   ksp->ops->solve                = KSPSolve_NASH;
758:   ksp->ops->destroy              = KSPDestroy_NASH;
759:   ksp->ops->setfromoptions       = KSPSetFromOptions_NASH;
760:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
761:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;
762:   ksp->ops->view                 = 0;

764:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
765:                                            "KSPNASHSetRadius_C",
766:                                            "KSPNASHSetRadius_NASH",
767:                                             KSPNASHSetRadius_NASH);
768:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
769:                                            "KSPNASHGetNormD_C",
770:                                            "KSPNASHGetNormD_NASH",
771:                                             KSPNASHGetNormD_NASH);
772:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
773:                                            "KSPNASHGetObjFcn_C",
774:                                            "KSPNASHGetObjFcn_NASH",
775:                                             KSPNASHGetObjFcn_NASH);
776:   return(0);
777: }
778: EXTERN_C_END