Actual source code: nash.c

petsc-master 2016-08-27
Report Typos and Errors
 2:  #include <../src/ksp/ksp/impls/cg/nash/nashimpl.h>

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

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

 12: /*@
 13:     KSPNASHSetRadius - Sets the radius of the trust region.

 15:     Logically Collective on KSP

 17:     Input Parameters:
 18: +   ksp    - the iterative context
 19: -   radius - the trust region radius (Infinity is the default)

 21:     Options Database Key:
 22: .   -ksp_nash_radius <r>

 24:     Level: advanced

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

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

 42: /*@
 43:     KSPNASHGetNormD - Got norm of the direction.

 45:     Collective on KSP

 47:     Input Parameters:
 48: +   ksp    - the iterative context
 49: -   norm_d - the norm of the direction

 51:     Level: advanced

 53: .keywords: KSP, NASH, get, norm direction
 54: @*/
 55: PetscErrorCode  KSPNASHGetNormD(KSP ksp, PetscReal *norm_d)
 56: {

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

 67: /*@
 68:     KSPNASHGetObjFcn - Get objective function value.

 70:     Collective on KSP

 72:     Input Parameters:
 73: +   ksp   - the iterative context
 74: -   o_fcn - the objective function value

 76:     Level: advanced

 78: .keywords: KSP, NASH, get, objective function
 79: @*/
 80: PetscErrorCode  KSPNASHGetObjFcn(KSP ksp, PetscReal *o_fcn)
 81: {

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

 92: static PetscErrorCode KSPSolve_NASH(KSP ksp)
 93: {
 94: #if defined(PETSC_USE_COMPLEX)
 95:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "NASH is not available for complex systems");
 96: #else
 97:   KSP_NASH       *cg = (KSP_NASH*)ksp->data;
 99:   Mat            Qmat, Mmat;
100:   Vec            r, z, p, d;
101:   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;

107:   PetscInt max_cg_its;

109:   PetscBool diagonalscale;

112:   /***************************************************************************/
113:   /* Check the arguments and parameters.                                     */
114:   /***************************************************************************/

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

120:   /***************************************************************************/
121:   /* Get the workspace vectors and initialize variables                      */
122:   /***************************************************************************/

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

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

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

137:   /***************************************************************************/
138:   /* Initialize objective function and direction.                            */
139:   /***************************************************************************/

141:   cg->o_fcn = 0.0;

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

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

152:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
153:   VecDot(r, r, &rr);               /* rr = r^T r        */
154:   KSPCheckDot(ksp,rr);

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

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

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

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

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

184:       /***********************************************************************/
185:       /* Compute objective function.                                         */
186:       /***********************************************************************/

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

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

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

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

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

219:       /***********************************************************************/
220:       /* Compute objective function.                                         */
221:       /***********************************************************************/

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

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

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

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

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

252:   default:
253:     norm_r = 0.0;
254:     break;
255:   }

257:   KSPLogResidualHistory(ksp, norm_r);
258:   KSPMonitor(ksp, ksp->its, norm_r);
259:   ksp->rnorm = norm_r;

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

263:   /***************************************************************************/
264:   /* Compute the first direction and update the iteration.                   */
265:   /***************************************************************************/

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

271:   /***************************************************************************/
272:   /* Check the matrix for numerical problems.                                */
273:   /***************************************************************************/

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

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

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

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

297:       /***********************************************************************/
298:       /* Compute objective function.                                         */
299:       /***********************************************************************/

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

310:   /***************************************************************************/
311:   /* Initialize variables for calculating the norm of the direction.         */
312:   /***************************************************************************/

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

321:   default:
322:     VecDot(p, p, &norm_p);
323:     break;
324:   }

326:   /***************************************************************************/
327:   /* Check for negative curvature.                                           */
328:   /***************************************************************************/

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

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

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

346:       step       = PetscSqrtReal(r2 / norm_p);
347:       cg->norm_d = cg->radius;

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

351:       /***********************************************************************/
352:       /* Update objective function.                                          */
353:       /***********************************************************************/

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

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

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

372:       /***********************************************************************/
373:       /* Compute objective function.                                         */
374:       /***********************************************************************/

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

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

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

397:     alpha = rz / kappa;

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

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

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

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

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

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

424:         /*********************************************************************/
425:         /* Update objective function.                                        */
426:         /*********************************************************************/

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

437:     /*************************************************************************/
438:     /* Now we can update the direction and residual.                         */
439:     /*************************************************************************/

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

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

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

456:     /*************************************************************************/
457:     /* Update objective function.                                            */
458:     /*************************************************************************/

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

462:     /*************************************************************************/
463:     /* Check that the preconditioner appears positive semidefinite.          */
464:     /*************************************************************************/

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

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

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

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

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

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

496:     default:
497:       norm_r = 0.;
498:       break;
499:     }

501:     KSPLogResidualHistory(ksp, norm_r);
502:     KSPMonitor(ksp, ksp->its, norm_r);
503:     ksp->rnorm = norm_r;

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

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

515:     /*************************************************************************/
516:     /* We have not converged yet.  Check for breakdown.                      */
517:     /*************************************************************************/

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

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

530:     /*************************************************************************/
531:     /* Check iteration limit.                                                */
532:     /*************************************************************************/

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

540:     /*************************************************************************/
541:     /* Update p and the norms.                                               */
542:     /*************************************************************************/

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

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

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

558:     /*************************************************************************/
559:     /* Compute the new direction and update the iteration.                   */
560:     /*************************************************************************/

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

566:     /*************************************************************************/
567:     /* Check for negative curvature.                                         */
568:     /*************************************************************************/

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

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

587: static PetscErrorCode KSPSetUp_NASH(KSP ksp)
588: {

592:   /***************************************************************************/
593:   /* Set work vectors needed by conjugate gradient method and allocate       */
594:   /***************************************************************************/

596:   KSPSetWorkVecs(ksp, 3);
597:   return(0);
598: }

602: static PetscErrorCode KSPDestroy_NASH(KSP ksp)
603: {

607:   /***************************************************************************/
608:   /* Clear composed functions                                                */
609:   /***************************************************************************/

611:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHSetRadius_C",NULL);
612:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetNormD_C",NULL);
613:   PetscObjectComposeFunction((PetscObject)ksp,"KSPNASHGetObjFcn_C",NULL);

615:   /***************************************************************************/
616:   /* Destroy KSP object.                                                     */
617:   /***************************************************************************/

619:   KSPDestroyDefault(ksp);
620:   return(0);
621: }

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

630:   cg->radius = radius;
631:   return(0);
632: }

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

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

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

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

658: static PetscErrorCode KSPSetFromOptions_NASH(PetscOptionItems *PetscOptionsObject,KSP ksp)
659: {
661:   KSP_NASH       *cg = (KSP_NASH*)ksp->data;

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

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

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

670:   PetscOptionsTail();
671:   return(0);
672: }

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

679:    Options Database Keys:
680: .      -ksp_nash_radius <r> - Trust Region Radius

682:    Notes: This is rarely used directly

684:    Level: developer

686:   Use preconditioned conjugate gradient to compute
687:   an approximate minimizer of the quadratic function

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

691:    subject to the trust region constraint

693:             || s || <= delta,

695:    where

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

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

707:   Notes:
708:   The preconditioner supplied should be symmetric and positive definite.

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

715: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
716: {
718:   KSP_NASH       *cg;

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

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

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

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

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