Actual source code: nash.c

petsc-master 2017-04-20
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"};

 10: static PetscErrorCode KSPCGSolve_NASH(KSP ksp)
 11: {
 12: #if defined(PETSC_USE_COMPLEX)
 13:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "NASH is not available for complex systems");
 14: #else
 15:   KSPCG_NASH     *cg = (KSPCG_NASH*)ksp->data;
 17:   Mat            Qmat, Mmat;
 18:   Vec            r, z, p, d;
 19:   PC             pc;

 21:   PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
 22:   PetscReal alpha, beta, kappa, rz, rzm1;
 23:   PetscReal rr, r2, step;

 25:   PetscInt max_cg_its;

 27:   PetscBool diagonalscale;

 30:   /***************************************************************************/
 31:   /* Check the arguments and parameters.                                     */
 32:   /***************************************************************************/

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

 38:   /***************************************************************************/
 39:   /* Get the workspace vectors and initialize variables                      */
 40:   /***************************************************************************/

 42:   r2 = cg->radius * cg->radius;
 43:   r  = ksp->work[0];
 44:   z  = ksp->work[1];
 45:   p  = ksp->work[2];
 46:   d  = ksp->vec_sol;
 47:   pc = ksp->pc;

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

 51:   VecGetSize(d, &max_cg_its);
 52:   max_cg_its = PetscMin(max_cg_its, ksp->max_it);
 53:   ksp->its   = 0;

 55:   /***************************************************************************/
 56:   /* Initialize objective function and direction.                            */
 57:   /***************************************************************************/

 59:   cg->o_fcn = 0.0;

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

 64:   /***************************************************************************/
 65:   /* Begin the conjugate gradient method.  Check the right-hand side for     */
 66:   /* numerical problems.  The check for not-a-number and infinite values     */
 67:   /* need be performed only once.                                            */
 68:   /***************************************************************************/

 70:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
 71:   VecDot(r, r, &rr);               /* rr = r^T r        */
 72:   KSPCheckDot(ksp,rr);

 74:   /***************************************************************************/
 75:   /* Check the preconditioner for numerical problems and for positive        */
 76:   /* definiteness.  The check for not-a-number and infinite values need be   */
 77:   /* performed only once.                                                    */
 78:   /***************************************************************************/

 80:   KSP_PCApply(ksp, r, z);          /* z = inv(M) r      */
 81:   VecDot(r, z, &rz);               /* rz = r^T inv(M) r */
 82:   if (PetscIsInfOrNanScalar(rz)) {
 83:     /*************************************************************************/
 84:     /* The preconditioner contains not-a-number or an infinite value.        */
 85:     /* Return the gradient direction intersected with the trust region.      */
 86:     /*************************************************************************/

 88:     ksp->reason = KSP_DIVERGED_NANORINF;
 89:     PetscInfo1(ksp, "KSPCGSolve_NASH: bad preconditioner: rz=%g\n", (double)rz);

 91:     if (cg->radius) {
 92:       if (r2 >= rr) {
 93:         alpha      = 1.0;
 94:         cg->norm_d = PetscSqrtReal(rr);
 95:       } else {
 96:         alpha      = PetscSqrtReal(r2 / rr);
 97:         cg->norm_d = cg->radius;
 98:       }

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

102:       /***********************************************************************/
103:       /* Compute objective function.                                         */
104:       /***********************************************************************/

106:       KSP_MatMult(ksp, Qmat, d, z);
107:       VecAYPX(z, -0.5, ksp->vec_rhs);
108:       VecDot(d, z, &cg->o_fcn);
109:       cg->o_fcn = -cg->o_fcn;
110:       ++ksp->its;
111:     }
112:     return(0);
113:   }

115:   if (rz < 0.0) {
116:     /*************************************************************************/
117:     /* The preconditioner is indefinite.  Because this is the first          */
118:     /* and we do not have a direction yet, we use the gradient step.  Note   */
119:     /* that we cannot use the preconditioned norm when computing the step    */
120:     /* because the matrix is indefinite.                                     */
121:     /*************************************************************************/

123:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
124:     PetscInfo1(ksp, "KSPCGSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz);

126:     if (cg->radius) {
127:       if (r2 >= rr) {
128:         alpha      = 1.0;
129:         cg->norm_d = PetscSqrtReal(rr);
130:       } else {
131:         alpha      = PetscSqrtReal(r2 / rr);
132:         cg->norm_d = cg->radius;
133:       }

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

137:       /***********************************************************************/
138:       /* Compute objective function.                                         */
139:       /***********************************************************************/

141:       KSP_MatMult(ksp, Qmat, d, z);
142:       VecAYPX(z, -0.5, ksp->vec_rhs);
143:       VecDot(d, z, &cg->o_fcn);
144:       cg->o_fcn = -cg->o_fcn;
145:       ++ksp->its;
146:     }
147:     return(0);
148:   }

150:   /***************************************************************************/
151:   /* As far as we know, the preconditioner is positive semidefinite.         */
152:   /* Compute and log the residual.  Check convergence because this           */
153:   /* initializes things, but do not terminate until at least one conjugate   */
154:   /* gradient iteration has been performed.                                  */
155:   /***************************************************************************/

157:   switch (ksp->normtype) {
158:   case KSP_NORM_PRECONDITIONED:
159:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
160:     break;

162:   case KSP_NORM_UNPRECONDITIONED:
163:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
164:     break;

166:   case KSP_NORM_NATURAL:
167:     norm_r = PetscSqrtReal(rz);                                 /* norm_r = |r|_M    */
168:     break;

170:   default:
171:     norm_r = 0.0;
172:     break;
173:   }

175:   KSPLogResidualHistory(ksp, norm_r);
176:   KSPMonitor(ksp, ksp->its, norm_r);
177:   ksp->rnorm = norm_r;

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

181:   /***************************************************************************/
182:   /* Compute the first direction and update the iteration.                   */
183:   /***************************************************************************/

185:   VecCopy(z, p);                   /* p = z             */
186:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
187:   ++ksp->its;

189:   /***************************************************************************/
190:   /* Check the matrix for numerical problems.                                */
191:   /***************************************************************************/

193:   VecDot(p, z, &kappa);            /* kappa = p^T Q p   */
194:   if (PetscIsInfOrNanScalar(kappa)) {
195:     /*************************************************************************/
196:     /* The matrix produced not-a-number or an infinite value.  In this case, */
197:     /* we must stop and use the gradient direction.  This condition need     */
198:     /* only be checked once.                                                 */
199:     /*************************************************************************/

201:     ksp->reason = KSP_DIVERGED_NANORINF;
202:     PetscInfo1(ksp, "KSPCGSolve_NASH: bad matrix: kappa=%g\n", (double)kappa);

204:     if (cg->radius) {
205:       if (r2 >= rr) {
206:         alpha      = 1.0;
207:         cg->norm_d = PetscSqrtReal(rr);
208:       } else {
209:         alpha      = PetscSqrtReal(r2 / rr);
210:         cg->norm_d = cg->radius;
211:       }

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

215:       /***********************************************************************/
216:       /* Compute objective function.                                         */
217:       /***********************************************************************/

219:       KSP_MatMult(ksp, Qmat, d, z);
220:       VecAYPX(z, -0.5, ksp->vec_rhs);
221:       VecDot(d, z, &cg->o_fcn);
222:       cg->o_fcn = -cg->o_fcn;
223:       ++ksp->its;
224:     }
225:     return(0);
226:   }

228:   /***************************************************************************/
229:   /* Initialize variables for calculating the norm of the direction.         */
230:   /***************************************************************************/

232:   dMp    = 0.0;
233:   norm_d = 0.0;
234:   switch (cg->dtype) {
235:   case NASH_PRECONDITIONED_DIRECTION:
236:     norm_p = rz;
237:     break;

239:   default:
240:     VecDot(p, p, &norm_p);
241:     break;
242:   }

244:   /***************************************************************************/
245:   /* Check for negative curvature.                                           */
246:   /***************************************************************************/

248:   if (kappa <= 0.0) {
249:     /*************************************************************************/
250:     /* In this case, the matrix is indefinite and we have encountered a      */
251:     /* direction of negative curvature.  Because negative curvature occurs   */
252:     /* during the first step, we must follow a direction.                    */
253:     /*************************************************************************/

255:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
256:     PetscInfo1(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);

258:     if (cg->radius && norm_p > 0.0) {
259:       /***********************************************************************/
260:       /* Follow direction of negative curvature to the boundary of the       */
261:       /* trust region.                                                       */
262:       /***********************************************************************/

264:       step       = PetscSqrtReal(r2 / norm_p);
265:       cg->norm_d = cg->radius;

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

269:       /***********************************************************************/
270:       /* Update objective function.                                          */
271:       /***********************************************************************/

273:       cg->o_fcn += step * (0.5 * step * kappa - rz);
274:     } else if (cg->radius) {
275:       /***********************************************************************/
276:       /* The norm of the preconditioned direction is zero; use the gradient  */
277:       /* step.                                                               */
278:       /***********************************************************************/

280:       if (r2 >= rr) {
281:         alpha      = 1.0;
282:         cg->norm_d = PetscSqrtReal(rr);
283:       } else {
284:         alpha      = PetscSqrtReal(r2 / rr);
285:         cg->norm_d = cg->radius;
286:       }

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

290:       /***********************************************************************/
291:       /* Compute objective function.                                         */
292:       /***********************************************************************/

294:       KSP_MatMult(ksp, Qmat, d, z);
295:       VecAYPX(z, -0.5, ksp->vec_rhs);
296:       VecDot(d, z, &cg->o_fcn);
297:       cg->o_fcn = -cg->o_fcn;
298:       ++ksp->its;
299:     }
300:     return(0);
301:   }

303:   /***************************************************************************/
304:   /* Run the conjugate gradient method until either the problem is solved,   */
305:   /* we encounter the boundary of the trust region, or the conjugate         */
306:   /* gradient method breaks down.                                            */
307:   /***************************************************************************/

309:   while (1) {
310:     /*************************************************************************/
311:     /* Know that kappa is nonzero, because we have not broken down, so we    */
312:     /* can compute the steplength.                                           */
313:     /*************************************************************************/

315:     alpha = rz / kappa;

317:     /*************************************************************************/
318:     /* Compute the steplength and check for intersection with the trust      */
319:     /* region.                                                               */
320:     /*************************************************************************/

322:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
323:     if (cg->radius && norm_dp1 >= r2) {
324:       /***********************************************************************/
325:       /* In this case, the matrix is positive definite as far as we know.    */
326:       /* However, the full step goes beyond the trust region.                */
327:       /***********************************************************************/

329:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
330:       PetscInfo1(ksp, "KSPCGSolve_NASH: constrained step: radius=%g\n", (double)cg->radius);

332:       if (norm_p > 0.0) {
333:         /*********************************************************************/
334:         /* Follow the direction to the boundary of the trust region.         */
335:         /*********************************************************************/

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

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

342:         /*********************************************************************/
343:         /* Update objective function.                                        */
344:         /*********************************************************************/

346:         cg->o_fcn += step * (0.5 * step * kappa - rz);
347:       } else {
348:         /*********************************************************************/
349:         /* The norm of the direction is zero; there is nothing to follow.    */
350:         /*********************************************************************/
351:       }
352:       break;
353:     }

355:     /*************************************************************************/
356:     /* Now we can update the direction and residual.                         */
357:     /*************************************************************************/

359:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
360:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
361:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

363:     switch (cg->dtype) {
364:     case NASH_PRECONDITIONED_DIRECTION:
365:       norm_d = norm_dp1;
366:       break;

368:     default:
369:       VecDot(d, d, &norm_d);
370:       break;
371:     }
372:     cg->norm_d = PetscSqrtReal(norm_d);

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

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

380:     /*************************************************************************/
381:     /* Check that the preconditioner appears positive semidefinite.          */
382:     /*************************************************************************/

384:     rzm1 = rz;
385:     VecDot(r, z, &rz);             /* rz = r^T z        */
386:     if (rz < 0.0) {
387:       /***********************************************************************/
388:       /* The preconditioner is indefinite.                                   */
389:       /***********************************************************************/

391:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
392:       PetscInfo1(ksp, "KSPCGSolve_NASH: cg indefinite preconditioner: rz=%g\n", (double)rz);
393:       break;
394:     }

396:     /*************************************************************************/
397:     /* As far as we know, the preconditioner is positive semidefinite.       */
398:     /* Compute the residual and check for convergence.                       */
399:     /*************************************************************************/

401:     switch (ksp->normtype) {
402:     case KSP_NORM_PRECONDITIONED:
403:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
404:       break;

406:     case KSP_NORM_UNPRECONDITIONED:
407:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
408:       break;

410:     case KSP_NORM_NATURAL:
411:       norm_r = PetscSqrtReal(rz);                               /* norm_r = |r|_M    */
412:       break;

414:     default:
415:       norm_r = 0.;
416:       break;
417:     }

419:     KSPLogResidualHistory(ksp, norm_r);
420:     KSPMonitor(ksp, ksp->its, norm_r);
421:     ksp->rnorm = norm_r;

423:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
424:     if (ksp->reason) {
425:       /***********************************************************************/
426:       /* The method has converged.                                           */
427:       /***********************************************************************/

429:       PetscInfo2(ksp, "KSPCGSolve_NASH: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
430:       break;
431:     }

433:     /*************************************************************************/
434:     /* We have not converged yet.  Check for breakdown.                      */
435:     /*************************************************************************/

437:     beta = rz / rzm1;
438:     if (PetscAbsReal(beta) <= 0.0) {
439:       /***********************************************************************/
440:       /* Conjugate gradients has broken down.                                */
441:       /***********************************************************************/

443:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
444:       PetscInfo1(ksp, "KSPCGSolve_NASH: breakdown: beta=%g\n", (double)beta);
445:       break;
446:     }

448:     /*************************************************************************/
449:     /* Check iteration limit.                                                */
450:     /*************************************************************************/

452:     if (ksp->its >= max_cg_its) {
453:       ksp->reason = KSP_DIVERGED_ITS;
454:       PetscInfo1(ksp, "KSPCGSolve_NASH: iterlim: its=%D\n", ksp->its);
455:       break;
456:     }

458:     /*************************************************************************/
459:     /* Update p and the norms.                                               */
460:     /*************************************************************************/

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

464:     switch (cg->dtype) {
465:     case NASH_PRECONDITIONED_DIRECTION:
466:       dMp    = beta*(dMp + alpha*norm_p);
467:       norm_p = beta*(rzm1 + beta*norm_p);
468:       break;

470:     default:
471:       VecDot(d, p, &dMp);
472:       VecDot(p, p, &norm_p);
473:       break;
474:     }

476:     /*************************************************************************/
477:     /* Compute the new direction and update the iteration.                   */
478:     /*************************************************************************/

480:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
481:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
482:     ++ksp->its;

484:     /*************************************************************************/
485:     /* Check for negative curvature.                                         */
486:     /*************************************************************************/

488:     if (kappa <= 0.0) {
489:       /***********************************************************************/
490:       /* In this case, the matrix is indefinite and we have encountered      */
491:       /* a direction of negative curvature.  Stop at the base.               */
492:       /***********************************************************************/

494:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
495:       PetscInfo1(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);
496:       break;
497:     }
498:   }
499:   return(0);
500: #endif
501: }

503: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
504: {

507:   /***************************************************************************/
508:   /* Set work vectors needed by conjugate gradient method and allocate       */
509:   /***************************************************************************/

512:   KSPSetWorkVecs(ksp,3);
513:   return(0);
514: }

516: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
517: {

521:   /***************************************************************************/
522:   /* Clear composed functions                                                */
523:   /***************************************************************************/

525:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
526:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
527:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);

529:   /***************************************************************************/
530:   /* Destroy KSP object.                                                     */
531:   /***************************************************************************/

533:   KSPDestroyDefault(ksp);
534:   return(0);
535: }

537: static PetscErrorCode  KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
538: {
539:   KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;

542:   cg->radius = radius;
543:   return(0);
544: }

546: static PetscErrorCode  KSPCGGetNormD_NASH(KSP ksp, PetscReal *norm_d)
547: {
548:   KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;

551:   *norm_d = cg->norm_d;
552:   return(0);
553: }

555: static PetscErrorCode  KSPCGGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
556: {
557:   KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;

560:   *o_fcn = cg->o_fcn;
561:   return(0);
562: }

564: static PetscErrorCode KSPCGSetFromOptions_NASH(PetscOptionItems *PetscOptionsObject,KSP ksp)
565: {
567:   KSPCG_NASH     *cg = (KSPCG_NASH*)ksp->data;

570:   PetscOptionsHead(PetscOptionsObject,"KSPCG NASH options");

572:   PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);

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

576:   PetscOptionsTail();
577:   return(0);
578: }

580: /*MC
581:      KSPCGNASH -   Code to run conjugate gradient method subject to a constraint
582:          on the solution norm. This is used in Trust Region methods for
583:          nonlinear equations, SNESNEWTONTR

585:    Options Database Keys:
586: .      -ksp_cg_radius <r> - Trust Region Radius

588:    Notes: This is rarely used directly

590:    Level: developer

592:   Use preconditioned conjugate gradient to compute
593:   an approximate minimizer of the quadratic function

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

597:    subject to the trust region constraint

599:             || s || <= delta,

601:    where

603:      delta is the trust region radius,
604:      g is the gradient vector,
605:      H is the Hessian approximation, and
606:      M is the positive definite preconditioner matrix.

608:    KSPConvergedReason may be
609: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
610: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
611: $  other KSP converged/diverged reasons

613:   Notes:
614:   The preconditioner supplied should be symmetric and positive definite.

616: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn()
617: M*/

619: PETSC_EXTERN PetscErrorCode KSPCreate_CGNASH(KSP ksp)
620: {
622:   KSPCG_NASH     *cg;

625:   PetscNewLog(ksp,&cg);
626:   cg->radius = 0.0;
627:   cg->dtype  = NASH_UNPRECONDITIONED_DIRECTION;

629:   ksp->data = (void*) cg;
630:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
631:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
632:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

634:   /***************************************************************************/
635:   /* Sets the functions that are associated with this data structure         */
636:   /* (in C++ this is the same as defining virtual functions).                */
637:   /***************************************************************************/

639:   ksp->ops->setup          = KSPCGSetUp_NASH;
640:   ksp->ops->solve          = KSPCGSolve_NASH;
641:   ksp->ops->destroy        = KSPCGDestroy_NASH;
642:   ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
643:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
644:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
645:   ksp->ops->view           = 0;

647:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_NASH);
648:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",KSPCGGetNormD_NASH);
649:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_NASH);
650:   return(0);
651: }