Actual source code: nash.c

petsc-master 2016-12-07
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: static PetscErrorCode KSPCGSolve_NASH(KSP ksp)
 13: {
 14: #if defined(PETSC_USE_COMPLEX)
 15:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "NASH is not available for complex systems");
 16: #else
 17:   KSPCG_NASH     *cg = (KSPCG_NASH*)ksp->data;
 19:   Mat            Qmat, Mmat;
 20:   Vec            r, z, p, d;
 21:   PC             pc;

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

 27:   PetscInt max_cg_its;

 29:   PetscBool diagonalscale;

 32:   /***************************************************************************/
 33:   /* Check the arguments and parameters.                                     */
 34:   /***************************************************************************/

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

 40:   /***************************************************************************/
 41:   /* Get the workspace vectors and initialize variables                      */
 42:   /***************************************************************************/

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

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

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

 57:   /***************************************************************************/
 58:   /* Initialize objective function and direction.                            */
 59:   /***************************************************************************/

 61:   cg->o_fcn = 0.0;

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

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

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

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

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

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

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

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

104:       /***********************************************************************/
105:       /* Compute objective function.                                         */
106:       /***********************************************************************/

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

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

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

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

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

139:       /***********************************************************************/
140:       /* Compute objective function.                                         */
141:       /***********************************************************************/

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

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

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

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

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

172:   default:
173:     norm_r = 0.0;
174:     break;
175:   }

177:   KSPLogResidualHistory(ksp, norm_r);
178:   KSPMonitor(ksp, ksp->its, norm_r);
179:   ksp->rnorm = norm_r;

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

183:   /***************************************************************************/
184:   /* Compute the first direction and update the iteration.                   */
185:   /***************************************************************************/

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

191:   /***************************************************************************/
192:   /* Check the matrix for numerical problems.                                */
193:   /***************************************************************************/

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

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

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

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

217:       /***********************************************************************/
218:       /* Compute objective function.                                         */
219:       /***********************************************************************/

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

230:   /***************************************************************************/
231:   /* Initialize variables for calculating the norm of the direction.         */
232:   /***************************************************************************/

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

241:   default:
242:     VecDot(p, p, &norm_p);
243:     break;
244:   }

246:   /***************************************************************************/
247:   /* Check for negative curvature.                                           */
248:   /***************************************************************************/

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

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

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

266:       step       = PetscSqrtReal(r2 / norm_p);
267:       cg->norm_d = cg->radius;

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

271:       /***********************************************************************/
272:       /* Update objective function.                                          */
273:       /***********************************************************************/

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

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

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

292:       /***********************************************************************/
293:       /* Compute objective function.                                         */
294:       /***********************************************************************/

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

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

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

317:     alpha = rz / kappa;

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

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

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

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

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

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

344:         /*********************************************************************/
345:         /* Update objective function.                                        */
346:         /*********************************************************************/

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

357:     /*************************************************************************/
358:     /* Now we can update the direction and residual.                         */
359:     /*************************************************************************/

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

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

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

376:     /*************************************************************************/
377:     /* Update objective function.                                            */
378:     /*************************************************************************/

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

382:     /*************************************************************************/
383:     /* Check that the preconditioner appears positive semidefinite.          */
384:     /*************************************************************************/

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

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

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

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

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

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

416:     default:
417:       norm_r = 0.;
418:       break;
419:     }

421:     KSPLogResidualHistory(ksp, norm_r);
422:     KSPMonitor(ksp, ksp->its, norm_r);
423:     ksp->rnorm = norm_r;

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

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

435:     /*************************************************************************/
436:     /* We have not converged yet.  Check for breakdown.                      */
437:     /*************************************************************************/

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

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

450:     /*************************************************************************/
451:     /* Check iteration limit.                                                */
452:     /*************************************************************************/

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

460:     /*************************************************************************/
461:     /* Update p and the norms.                                               */
462:     /*************************************************************************/

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

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

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

478:     /*************************************************************************/
479:     /* Compute the new direction and update the iteration.                   */
480:     /*************************************************************************/

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

486:     /*************************************************************************/
487:     /* Check for negative curvature.                                         */
488:     /*************************************************************************/

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

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

507: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
508: {

511:   /***************************************************************************/
512:   /* Set work vectors needed by conjugate gradient method and allocate       */
513:   /***************************************************************************/

516:   KSPSetWorkVecs(ksp,3);
517:   return(0);
518: }

522: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
523: {

527:   /***************************************************************************/
528:   /* Clear composed functions                                                */
529:   /***************************************************************************/

531:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
532:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
533:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);

535:   /***************************************************************************/
536:   /* Destroy KSP object.                                                     */
537:   /***************************************************************************/

539:   KSPDestroyDefault(ksp);
540:   return(0);
541: }

545: static PetscErrorCode  KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
546: {
547:   KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;

550:   cg->radius = radius;
551:   return(0);
552: }

556: static PetscErrorCode  KSPCGGetNormD_NASH(KSP ksp, PetscReal *norm_d)
557: {
558:   KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;

561:   *norm_d = cg->norm_d;
562:   return(0);
563: }

567: static PetscErrorCode  KSPCGGetObjFcn_NASH(KSP ksp, PetscReal *o_fcn)
568: {
569:   KSPCG_NASH *cg = (KSPCG_NASH*)ksp->data;

572:   *o_fcn = cg->o_fcn;
573:   return(0);
574: }

578: static PetscErrorCode KSPCGSetFromOptions_NASH(PetscOptionItems *PetscOptionsObject,KSP ksp)
579: {
581:   KSPCG_NASH     *cg = (KSPCG_NASH*)ksp->data;

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

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

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

590:   PetscOptionsTail();
591:   return(0);
592: }

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

599:    Options Database Keys:
600: .      -ksp_cg_radius <r> - Trust Region Radius

602:    Notes: This is rarely used directly

604:    Level: developer

606:   Use preconditioned conjugate gradient to compute
607:   an approximate minimizer of the quadratic function

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

611:    subject to the trust region constraint

613:             || s || <= delta,

615:    where

617:      delta is the trust region radius,
618:      g is the gradient vector,
619:      H is the Hessian approximation, and
620:      M is the positive definite preconditioner matrix.

622:    KSPConvergedReason may be
623: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
624: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
625: $  other KSP converged/diverged reasons

627:   Notes:
628:   The preconditioner supplied should be symmetric and positive definite.

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

635: PETSC_EXTERN PetscErrorCode KSPCreate_CGNASH(KSP ksp)
636: {
638:   KSPCG_NASH     *cg;

641:   PetscNewLog(ksp,&cg);
642:   cg->radius = 0.0;
643:   cg->dtype  = NASH_UNPRECONDITIONED_DIRECTION;

645:   ksp->data = (void*) cg;
646:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
647:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
648:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

650:   /***************************************************************************/
651:   /* Sets the functions that are associated with this data structure         */
652:   /* (in C++ this is the same as defining virtual functions).                */
653:   /***************************************************************************/

655:   ksp->ops->setup          = KSPCGSetUp_NASH;
656:   ksp->ops->solve          = KSPCGSolve_NASH;
657:   ksp->ops->destroy        = KSPCGDestroy_NASH;
658:   ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
659:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
660:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
661:   ksp->ops->view           = 0;

663:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_NASH);
664:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",KSPCGGetNormD_NASH);
665:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_NASH);
666:   return(0);
667: }