Actual source code: nash.c

petsc-3.4.4 2014-03-13
  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:   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(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
121:   if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),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_NANORINF;
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_NANORINF;
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:       } else {
190:         alpha      = PetscSqrtReal(r2 / rr);
191:         cg->norm_d = cg->radius;
192:       }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

409:     alpha = rz / kappa;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

599: PetscErrorCode KSPSetUp_NASH(KSP ksp)
600: {

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

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

614: PetscErrorCode KSPDestroy_NASH(KSP ksp)
615: {

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

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

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

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

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

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

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

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

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

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

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

676:   PetscOptionsHead("KSP NASH options");

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

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

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

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

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

694:    Notes: This is rarely used directly

696:    Level: developer

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

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

703:    subject to the trust region constraint

705:             || s || <= delta,

707:    where

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

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

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

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

727: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
728: {
730:   KSP_NASH       *cg;

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

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

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

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

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