Actual source code: nash.c

petsc-master 2015-07-31
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:   KSPCheckDot(ksp,rr);

158:   /***************************************************************************/
159:   /* Check the preconditioner for numerical problems and for positive        */
160:   /* definiteness.  The check for not-a-number and infinite values need be   */
161:   /* performed only once.                                                    */
162:   /***************************************************************************/

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

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

175:     if (cg->radius) {
176:       if (r2 >= rr) {
177:         alpha      = 1.0;
178:         cg->norm_d = PetscSqrtReal(rr);
179:       } else {
180:         alpha      = PetscSqrtReal(r2 / rr);
181:         cg->norm_d = cg->radius;
182:       }

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

186:       /***********************************************************************/
187:       /* Compute objective function.                                         */
188:       /***********************************************************************/

190:       KSP_MatMult(ksp, Qmat, d, z);
191:       VecAYPX(z, -0.5, ksp->vec_rhs);
192:       VecDot(d, z, &cg->o_fcn);
193:       cg->o_fcn = -cg->o_fcn;
194:       ++ksp->its;
195:     }
196:     return(0);
197:   }

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

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

210:     if (cg->radius) {
211:       if (r2 >= rr) {
212:         alpha      = 1.0;
213:         cg->norm_d = PetscSqrtReal(rr);
214:       } else {
215:         alpha      = PetscSqrtReal(r2 / rr);
216:         cg->norm_d = cg->radius;
217:       }

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

221:       /***********************************************************************/
222:       /* Compute objective function.                                         */
223:       /***********************************************************************/

225:       KSP_MatMult(ksp, Qmat, d, z);
226:       VecAYPX(z, -0.5, ksp->vec_rhs);
227:       VecDot(d, z, &cg->o_fcn);
228:       cg->o_fcn = -cg->o_fcn;
229:       ++ksp->its;
230:     }
231:     return(0);
232:   }

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

241:   switch (ksp->normtype) {
242:   case KSP_NORM_PRECONDITIONED:
243:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
244:     break;

246:   case KSP_NORM_UNPRECONDITIONED:
247:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
248:     break;

250:   case KSP_NORM_NATURAL:
251:     norm_r = PetscSqrtReal(rz);                                 /* norm_r = |r|_M    */
252:     break;

254:   default:
255:     norm_r = 0.0;
256:     break;
257:   }

259:   KSPLogResidualHistory(ksp, norm_r);
260:   KSPMonitor(ksp, ksp->its, norm_r);
261:   ksp->rnorm = norm_r;

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

265:   /***************************************************************************/
266:   /* Compute the first direction and update the iteration.                   */
267:   /***************************************************************************/

269:   VecCopy(z, p);                   /* p = z             */
270:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
271:   ++ksp->its;

273:   /***************************************************************************/
274:   /* Check the matrix for numerical problems.                                */
275:   /***************************************************************************/

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

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

288:     if (cg->radius) {
289:       if (r2 >= rr) {
290:         alpha      = 1.0;
291:         cg->norm_d = PetscSqrtReal(rr);
292:       } else {
293:         alpha      = PetscSqrtReal(r2 / rr);
294:         cg->norm_d = cg->radius;
295:       }

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

299:       /***********************************************************************/
300:       /* Compute objective function.                                         */
301:       /***********************************************************************/

303:       KSP_MatMult(ksp, Qmat, d, z);
304:       VecAYPX(z, -0.5, ksp->vec_rhs);
305:       VecDot(d, z, &cg->o_fcn);
306:       cg->o_fcn = -cg->o_fcn;
307:       ++ksp->its;
308:     }
309:     return(0);
310:   }

312:   /***************************************************************************/
313:   /* Initialize variables for calculating the norm of the direction.         */
314:   /***************************************************************************/

316:   dMp    = 0.0;
317:   norm_d = 0.0;
318:   switch (cg->dtype) {
319:   case NASH_PRECONDITIONED_DIRECTION:
320:     norm_p = rz;
321:     break;

323:   default:
324:     VecDot(p, p, &norm_p);
325:     break;
326:   }

328:   /***************************************************************************/
329:   /* Check for negative curvature.                                           */
330:   /***************************************************************************/

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

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

342:     if (cg->radius && norm_p > 0.0) {
343:       /***********************************************************************/
344:       /* Follow direction of negative curvature to the boundary of the       */
345:       /* trust region.                                                       */
346:       /***********************************************************************/

348:       step       = PetscSqrtReal(r2 / norm_p);
349:       cg->norm_d = cg->radius;

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

353:       /***********************************************************************/
354:       /* Update objective function.                                          */
355:       /***********************************************************************/

357:       cg->o_fcn += step * (0.5 * step * kappa - rz);
358:     } else if (cg->radius) {
359:       /***********************************************************************/
360:       /* The norm of the preconditioned direction is zero; use the gradient  */
361:       /* step.                                                               */
362:       /***********************************************************************/

364:       if (r2 >= rr) {
365:         alpha      = 1.0;
366:         cg->norm_d = PetscSqrtReal(rr);
367:       } else {
368:         alpha      = PetscSqrtReal(r2 / rr);
369:         cg->norm_d = cg->radius;
370:       }

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

374:       /***********************************************************************/
375:       /* Compute objective function.                                         */
376:       /***********************************************************************/

378:       KSP_MatMult(ksp, Qmat, d, z);
379:       VecAYPX(z, -0.5, ksp->vec_rhs);
380:       VecDot(d, z, &cg->o_fcn);
381:       cg->o_fcn = -cg->o_fcn;
382:       ++ksp->its;
383:     }
384:     return(0);
385:   }

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

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

399:     alpha = rz / kappa;

401:     /*************************************************************************/
402:     /* Compute the steplength and check for intersection with the trust      */
403:     /* region.                                                               */
404:     /*************************************************************************/

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

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

416:       if (norm_p > 0.0) {
417:         /*********************************************************************/
418:         /* Follow the direction to the boundary of the trust region.         */
419:         /*********************************************************************/

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

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

426:         /*********************************************************************/
427:         /* Update objective function.                                        */
428:         /*********************************************************************/

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

439:     /*************************************************************************/
440:     /* Now we can update the direction and residual.                         */
441:     /*************************************************************************/

443:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
444:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
445:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

447:     switch (cg->dtype) {
448:     case NASH_PRECONDITIONED_DIRECTION:
449:       norm_d = norm_dp1;
450:       break;

452:     default:
453:       VecDot(d, d, &norm_d);
454:       break;
455:     }
456:     cg->norm_d = PetscSqrtReal(norm_d);

458:     /*************************************************************************/
459:     /* Update objective function.                                            */
460:     /*************************************************************************/

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

464:     /*************************************************************************/
465:     /* Check that the preconditioner appears positive semidefinite.          */
466:     /*************************************************************************/

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

475:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
476:       PetscInfo1(ksp, "KSPSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz);
477:       break;
478:     }

480:     /*************************************************************************/
481:     /* As far as we know, the preconditioner is positive semidefinite.       */
482:     /* Compute the residual and check for convergence.                       */
483:     /*************************************************************************/

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

490:     case KSP_NORM_UNPRECONDITIONED:
491:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
492:       break;

494:     case KSP_NORM_NATURAL:
495:       norm_r = PetscSqrtReal(rz);                               /* norm_r = |r|_M    */
496:       break;

498:     default:
499:       norm_r = 0.;
500:       break;
501:     }

503:     KSPLogResidualHistory(ksp, norm_r);
504:     KSPMonitor(ksp, ksp->its, norm_r);
505:     ksp->rnorm = norm_r;

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

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

517:     /*************************************************************************/
518:     /* We have not converged yet.  Check for breakdown.                      */
519:     /*************************************************************************/

521:     beta = rz / rzm1;
522:     if (PetscAbsReal(beta) <= 0.0) {
523:       /***********************************************************************/
524:       /* Conjugate gradients has broken down.                                */
525:       /***********************************************************************/

527:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
528:       PetscInfo1(ksp, "KSPSolve_NASH: breakdown: beta=%g\n", (double)beta);
529:       break;
530:     }

532:     /*************************************************************************/
533:     /* Check iteration limit.                                                */
534:     /*************************************************************************/

536:     if (ksp->its >= max_cg_its) {
537:       ksp->reason = KSP_DIVERGED_ITS;
538:       PetscInfo1(ksp, "KSPSolve_NASH: iterlim: its=%D\n", ksp->its);
539:       break;
540:     }

542:     /*************************************************************************/
543:     /* Update p and the norms.                                               */
544:     /*************************************************************************/

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

548:     switch (cg->dtype) {
549:     case NASH_PRECONDITIONED_DIRECTION:
550:       dMp    = beta*(dMp + alpha*norm_p);
551:       norm_p = beta*(rzm1 + beta*norm_p);
552:       break;

554:     default:
555:       VecDot(d, p, &dMp);
556:       VecDot(p, p, &norm_p);
557:       break;
558:     }

560:     /*************************************************************************/
561:     /* Compute the new direction and update the iteration.                   */
562:     /*************************************************************************/

564:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
565:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
566:     ++ksp->its;

568:     /*************************************************************************/
569:     /* Check for negative curvature.                                         */
570:     /*************************************************************************/

572:     if (kappa <= 0.0) {
573:       /***********************************************************************/
574:       /* In this case, the matrix is indefinite and we have encountered      */
575:       /* a direction of negative curvature.  Stop at the base.               */
576:       /***********************************************************************/

578:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
579:       PetscInfo1(ksp, "KSPSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
580:       break;
581:     }
582:   }
583:   return(0);
584: #endif
585: }

589: PetscErrorCode KSPSetUp_NASH(KSP ksp)
590: {

594:   /***************************************************************************/
595:   /* Set work vectors needed by conjugate gradient method and allocate       */
596:   /***************************************************************************/

598:   KSPSetWorkVecs(ksp, 3);
599:   return(0);
600: }

604: PetscErrorCode KSPDestroy_NASH(KSP ksp)
605: {

609:   /***************************************************************************/
610:   /* Clear composed functions                                                */
611:   /***************************************************************************/

613:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",NULL);
614:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",NULL);
615:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",NULL);

617:   /***************************************************************************/
618:   /* Destroy KSP object.                                                     */
619:   /***************************************************************************/

621:   KSPDestroyDefault(ksp);
622:   return(0);
623: }

627: static PetscErrorCode  KSPNASHSetRadius_NASH(KSP ksp, PetscReal radius)
628: {
629:   KSP_NASH *cg = (KSP_NASH*)ksp->data;

632:   cg->radius = radius;
633:   return(0);
634: }

638: static PetscErrorCode  KSPNASHGetNormD_NASH(KSP ksp, PetscReal *norm_d)
639: {
640:   KSP_NASH *cg = (KSP_NASH*)ksp->data;

643:   *norm_d = cg->norm_d;
644:   return(0);
645: }

649: static PetscErrorCode  KSPNASHGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
650: {
651:   KSP_NASH *cg = (KSP_NASH*)ksp->data;

654:   *o_fcn = cg->o_fcn;
655:   return(0);
656: }

660: PetscErrorCode KSPSetFromOptions_NASH(PetscOptions *PetscOptionsObject,KSP ksp)
661: {
663:   KSP_NASH       *cg = (KSP_NASH*)ksp->data;

666:   PetscOptionsHead(PetscOptionsObject,"KSP NASH options");

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

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

672:   PetscOptionsTail();
673:   return(0);
674: }

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

681:    Options Database Keys:
682: .      -ksp_nash_radius <r> - Trust Region Radius

684:    Notes: This is rarely used directly

686:    Level: developer

688:   Use preconditioned conjugate gradient to compute
689:   an approximate minimizer of the quadratic function

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

693:    subject to the trust region constraint

695:             || s || <= delta,

697:    where

699:      delta is the trust region radius,
700:      g is the gradient vector,
701:      H is the Hessian approximation, and
702:      M is the positive definite preconditioner matrix.

704:    KSPConvergedReason may be
705: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
706: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
707: $  other KSP converged/diverged reasons

709:   Notes:
710:   The preconditioner supplied should be symmetric and positive definite.

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

717: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
718: {
720:   KSP_NASH       *cg;

723:   PetscNewLog(ksp,&cg);
724:   cg->radius = 0.0;
725:   cg->dtype  = NASH_UNPRECONDITIONED_DIRECTION;

727:   ksp->data = (void*) cg;
728:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
729:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
730:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

732:   /***************************************************************************/
733:   /* Sets the functions that are associated with this data structure         */
734:   /* (in C++ this is the same as defining virtual functions).                */
735:   /***************************************************************************/

737:   ksp->ops->setup          = KSPSetUp_NASH;
738:   ksp->ops->solve          = KSPSolve_NASH;
739:   ksp->ops->destroy        = KSPDestroy_NASH;
740:   ksp->ops->setfromoptions = KSPSetFromOptions_NASH;
741:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
742:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
743:   ksp->ops->view           = 0;

745:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",KSPNASHSetRadius_NASH);
746:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",KSPNASHGetNormD_NASH);
747:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",KSPNASHGetObjFcn_NASH);
748:   return(0);
749: }