Actual source code: nash.c


  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;
 16:   Mat         Qmat, Mmat;
 17:   Vec         r, z, p, d;
 18:   PC          pc;

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

 24:   PetscInt max_cg_its;

 26:   PetscBool diagonalscale;

 28:   /***************************************************************************/
 29:   /* Check the arguments and parameters.                                     */
 30:   /***************************************************************************/

 32:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

 36:   /***************************************************************************/
 37:   /* Get the workspace vectors and initialize variables                      */
 38:   /***************************************************************************/

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

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

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

 53:   /***************************************************************************/
 54:   /* Initialize objective function and direction.                            */
 55:   /***************************************************************************/

 57:   cg->o_fcn = 0.0;

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

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

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

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

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

 86:     ksp->reason = KSP_DIVERGED_NANORINF;
 87:     PetscInfo(ksp, "KSPCGSolve_NASH: bad preconditioner: rz=%g\n", (double)rz);

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

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

100:       /***********************************************************************/
101:       /* Compute objective function.                                         */
102:       /***********************************************************************/

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

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

121:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
122:     PetscInfo(ksp, "KSPCGSolve_NASH: indefinite preconditioner: rz=%g\n", (double)rz);

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

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

135:       /***********************************************************************/
136:       /* Compute objective function.                                         */
137:       /***********************************************************************/

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

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

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

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

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

168:   default:
169:     norm_r = 0.0;
170:     break;
171:   }

173:   KSPLogResidualHistory(ksp, norm_r);
174:   KSPMonitor(ksp, ksp->its, norm_r);
175:   ksp->rnorm = norm_r;

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

179:   /***************************************************************************/
180:   /* Compute the first direction and update the iteration.                   */
181:   /***************************************************************************/

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

187:   /***************************************************************************/
188:   /* Check the matrix for numerical problems.                                */
189:   /***************************************************************************/

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

199:     ksp->reason = KSP_DIVERGED_NANORINF;
200:     PetscInfo(ksp, "KSPCGSolve_NASH: bad matrix: kappa=%g\n", (double)kappa);

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

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

213:       /***********************************************************************/
214:       /* Compute objective function.                                         */
215:       /***********************************************************************/

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

226:   /***************************************************************************/
227:   /* Initialize variables for calculating the norm of the direction.         */
228:   /***************************************************************************/

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

237:   default:
238:     VecDot(p, p, &norm_p);
239:     break;
240:   }

242:   /***************************************************************************/
243:   /* Check for negative curvature.                                           */
244:   /***************************************************************************/

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

253:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
254:     PetscInfo(ksp, "KSPCGSolve_NASH: negative curvature: kappa=%g\n", (double)kappa);

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

262:       step       = PetscSqrtReal(r2 / norm_p);
263:       cg->norm_d = cg->radius;

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

267:       /***********************************************************************/
268:       /* Update objective function.                                          */
269:       /***********************************************************************/

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

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

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

288:       /***********************************************************************/
289:       /* Compute objective function.                                         */
290:       /***********************************************************************/

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

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

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

313:     alpha = rz / kappa;

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

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

327:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
328:       PetscInfo(ksp, "KSPCGSolve_NASH: constrained step: radius=%g\n", (double)cg->radius);

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

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

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

340:         /*********************************************************************/
341:         /* Update objective function.                                        */
342:         /*********************************************************************/

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

353:     /*************************************************************************/
354:     /* Now we can update the direction and residual.                         */
355:     /*************************************************************************/

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

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

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

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

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

378:     /*************************************************************************/
379:     /* Check that the preconditioner appears positive semidefinite.          */
380:     /*************************************************************************/

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

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

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

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

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

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

412:     default:
413:       norm_r = 0.;
414:       break;
415:     }

417:     KSPLogResidualHistory(ksp, norm_r);
418:     KSPMonitor(ksp, ksp->its, norm_r);
419:     ksp->rnorm = norm_r;

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

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

431:     /*************************************************************************/
432:     /* We have not converged yet.  Check for breakdown.                      */
433:     /*************************************************************************/

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

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

446:     /*************************************************************************/
447:     /* Check iteration limit.                                                */
448:     /*************************************************************************/

450:     if (ksp->its >= max_cg_its) {
451:       ksp->reason = KSP_DIVERGED_ITS;
452:       PetscInfo(ksp, "KSPCGSolve_NASH: iterlim: its=%" PetscInt_FMT "\n", ksp->its);
453:       break;
454:     }

456:     /*************************************************************************/
457:     /* Update p and the norms.                                               */
458:     /*************************************************************************/

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

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

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

474:     /*************************************************************************/
475:     /* Compute the new direction and update the iteration.                   */
476:     /*************************************************************************/

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

482:     /*************************************************************************/
483:     /* Check for negative curvature.                                         */
484:     /*************************************************************************/

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

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

501: static PetscErrorCode KSPCGSetUp_NASH(KSP ksp)
502: {
503:   /***************************************************************************/
504:   /* Set work vectors needed by conjugate gradient method and allocate       */
505:   /***************************************************************************/

507:   KSPSetWorkVecs(ksp, 3);
508:   return 0;
509: }

511: static PetscErrorCode KSPCGDestroy_NASH(KSP ksp)
512: {
513:   /***************************************************************************/
514:   /* Clear composed functions                                                */
515:   /***************************************************************************/

517:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL);
518:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL);
519:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL);

521:   /***************************************************************************/
522:   /* Destroy KSP object.                                                     */
523:   /***************************************************************************/

525:   KSPDestroyDefault(ksp);
526:   return 0;
527: }

529: static PetscErrorCode KSPCGSetRadius_NASH(KSP ksp, PetscReal radius)
530: {
531:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;

533:   cg->radius = radius;
534:   return 0;
535: }

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

541:   *norm_d = cg->norm_d;
542:   return 0;
543: }

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

549:   *o_fcn = cg->o_fcn;
550:   return 0;
551: }

553: static PetscErrorCode KSPCGSetFromOptions_NASH(KSP ksp, PetscOptionItems *PetscOptionsObject)
554: {
555:   KSPCG_NASH *cg = (KSPCG_NASH *)ksp->data;

557:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG NASH options");

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

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

563:   PetscOptionsHeadEnd();
564:   return 0;
565: }

567: /*MC
568:      KSPNASH -   Code to run conjugate gradient method subject to a constraint on the solution norm.

570:    Options Database Keys:
571: .      -ksp_cg_radius <r> - Trust Region Radius

573:    Level: developer

575:    Notes:
576:     This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`

578:   Uses preconditioned conjugate gradient to compute
579:   an approximate minimizer of the quadratic function

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

583:    subject to the trust region constraint

585:             || s || <= delta,

587:    where

589:      delta is the trust region radius,
590:      g is the gradient vector,
591:      H is the Hessian approximation, and
592:      M is the positive definite preconditioner matrix.

594:    `KSPConvergedReason` may be
595: .vb
596:   KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
597:   KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
598: .ve
599:   other `KSP` converged/diverged reasons

601:   The preconditioner supplied should be symmetric and positive definite.

603:   Reference:
604:    Nash, Stephen G. Newton-type minimization via the Lanczos method. SIAM Journal on Numerical Analysis 21, no. 4 (1984): 770-788.

606: .seealso: [](chapter_ksp), `KSPQCG`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`
607: M*/

609: PETSC_EXTERN PetscErrorCode KSPCreate_NASH(KSP ksp)
610: {
611:   KSPCG_NASH *cg;

613:   PetscNew(&cg);
614:   cg->radius = 0.0;
615:   cg->dtype  = NASH_UNPRECONDITIONED_DIRECTION;

617:   ksp->data = (void *)cg;
618:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
619:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
620:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
621:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

623:   /***************************************************************************/
624:   /* Sets the functions that are associated with this data structure         */
625:   /* (in C++ this is the same as defining virtual functions).                */
626:   /***************************************************************************/

628:   ksp->ops->setup          = KSPCGSetUp_NASH;
629:   ksp->ops->solve          = KSPCGSolve_NASH;
630:   ksp->ops->destroy        = KSPCGDestroy_NASH;
631:   ksp->ops->setfromoptions = KSPCGSetFromOptions_NASH;
632:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
633:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
634:   ksp->ops->view           = NULL;

636:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_NASH);
637:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_NASH);
638:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_NASH);
639:   return 0;
640: }