Actual source code: stcg.c

  1: #define PETSCKSP_DLL

 3:  #include include/private/kspimpl.h
 4:  #include src/ksp/ksp/impls/cg/stcg/stcg.h

  8: /*@
  9:     KSPSTCGSetRadius - Sets the radius of the trust region.

 11:     Collective on KSP

 13:     Input Parameters:
 14: +   ksp    - the iterative context
 15: -   radius - the trust region radius (Infinity is the default)

 17:     Options Database Key:
 18: .   -ksp_stcg_radius <r>

 20:     Level: advanced

 22: .keywords: KSP, STCG, set, trust region radius
 23: @*/
 24: PetscErrorCode  KSPSTCGSetRadius(KSP ksp,PetscReal radius)
 25: {
 26:   PetscErrorCode ierr, (*f)(KSP, PetscReal);

 30:   if (radius < 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be positive");
 31:   PetscObjectQueryFunction((PetscObject)ksp, "KSPSTCGSetRadius_C", (void (**)(void))&f);
 32:   if (f) {
 33:     (*f)(ksp, radius);
 34:   }
 35:   return(0);
 36: }

 40: /*@
 41:     KSPSTCGGetNormD - Got norm of the direction.

 43:     Collective on KSP

 45:     Input Parameters:
 46: +   ksp    - the iterative context
 47: -   norm_d - the norm of the direction

 49:     Level: advanced

 51: .keywords: KSP, STCG, get, norm direction
 52: @*/
 53: PetscErrorCode  KSPSTCGGetNormD(KSP ksp,PetscReal *norm_d)
 54: {
 55:   PetscErrorCode ierr, (*f)(KSP, PetscReal *);

 59:   PetscObjectQueryFunction((PetscObject)ksp, "KSPSTCGGetNormD_C", (void (**)(void))&f);
 60:   if (f) {
 61:     (*f)(ksp, norm_d);
 62:   }
 63:   return(0);
 64: }

 68: /*
 69:   KSPSolve_STCG - Use preconditioned conjugate gradient to compute
 70:   an approximate minimizer of the quadratic function

 72:             q(s) = g^T * s + .5 * s^T * H * s

 74:    subject to the trust region constraint

 76:             || s ||_M <= delta,

 78:    where

 80:      delta is the trust region radius,
 81:      g is the gradient vector, and
 82:      H is the Hessian matrix,
 83:      M is the positive definite preconditioner matrix.

 85:    KSPConvergedReason may be
 86: $  KSP_CONVERGED_STCG_NEG_CURVE if convergence is reached along a negative curvature direction,
 87: $  KSP_CONVERGED_STCG_CONSTRAINED if convergence is reached along a constrained step,
 88: $  other KSP converged/diverged reasons


 91:   Notes:
 92:   The preconditioner supplied should be symmetric and positive definite.
 93: */

 95: #ifdef PETSC_USE_COMPLEX
 96: PetscErrorCode STCG_VecDot(Vec x, Vec y, PetscReal *a)
 97: {
 98:   PetscScalar ca;

101:   VecDot(x,y,&ca);
102:   *a = PetscRealPart(ca);
103:   return ierr;
104: }
105: #else
106: PetscErrorCode STCG_VecDot(Vec x, Vec y, PetscReal *a)
107: {
108:   return VecDot(x,y,a);
109: }
110: #endif

112: PetscErrorCode KSPSolve_STCG(KSP ksp)
113: {
115:   MatStructure  pflag;
116:   Mat Qmat, Mmat;
117:   Vec r, z, p, d;
118:   PC  pc;
119:   KSP_STCG *cg;
120:   PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
121:   PetscReal alpha, beta, kappa, rz, rzm1;
122:   PetscReal r2;
123:   PetscInt  i, maxit;
124:   PetscTruth diagonalscale;

127:   PCDiagonalScale(ksp->pc, &diagonalscale);
128:   if (diagonalscale) SETERRQ1(PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ksp->type_name);

130:   cg       = (KSP_STCG *)ksp->data;
131:   r2       = cg->radius * cg->radius;
132:   maxit    = ksp->max_it;
133:   r        = ksp->work[0];
134:   z        = ksp->work[1];
135:   p        = ksp->work[2];
136:   d        = ksp->vec_sol;
137:   pc       = ksp->pc;

139:   ksp->its = 0;
140:   if (cg->radius < 0.0) {
141:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
142:   }

144:   /* Initialize variables */
145:   PCGetOperators(pc, &Qmat, &Mmat, &pflag);

147:   VecSet(d, 0.0);                         /* d = 0        */
148:   VecCopy(ksp->vec_rhs, r);         /* r = -grad    */
149:   KSP_PCApply(ksp, r, z);                 /* z = M_{-1} r */
150:   cg->norm_d = 0.0;

152:   /* Check for numerical problems with the preconditioner */
153:   STCG_VecDot(r, z, &rz);                 /* rz = r^T z   */
154:   if ((rz != rz) || (rz && (rz / rz != rz / rz))) {
155:     ksp->reason = KSP_DIVERGED_NAN;
156:     PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g\n", rz);

158:     /* In this case, the preconditioner produced not a number or an         */
159:     /* infinite value.  We just take the gradient step.                     */
160:     /* Only needs to be checked once.                                       */

162:     if (cg->radius) {
163:       STCG_VecDot(r, r, &rz);         /* rz = r^T r   */

165:       alpha = sqrt(r2 / rz);
166:       VecAXPY(d, alpha, r);         /* d = d + alpha r */
167:       cg->norm_d = cg->radius;
168:     }
169:     return(0);
170:   }

172:   /* Check that the preconditioner is positive definite */
173:   if (rz <= 0.0) {
174:     VecNorm(r, NORM_2, &norm_r);
175:     if (rz < 0.0 || norm_r > 0.0) {
176:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
177:       PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", rz);

179:       /* In this case, the preconditioner is indefinite, so we cannot measure */
180:       /* the direction in the preconditioned norm.  Therefore, we must use    */
181:       /* an unpreconditioned calculation.  The direction in this case is      */
182:       /* uses the right hand side, which should be the negative gradient      */
183:       /* intersected with the trust region.                                   */

185:       if (cg->radius) {
186:         STCG_VecDot(r, r, &rz);         /* rz = r^T r   */
187: 
188:         alpha = sqrt(r2 / rz);
189:         VecAXPY(d, alpha, r);         /* d = d + alpha r */
190:         cg->norm_d = cg->radius;
191:       }
192:       return(0);
193:     }
194:   }

196:   /* As far as we know, the preconditioner is positive definite.  Compute   */
197:   /* the appropriate residual depending on what the user has set.           */
198:   if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
199:     VecNorm(z, NORM_2, &norm_r);         /* norm_r = |z| */
200:   }
201:   else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
202:     VecNorm(r, NORM_2, &norm_r);         /* norm_r = |r| */
203:   }
204:   else if (ksp->normtype == KSP_NATURAL_NORM) {
205:     norm_r = sqrt(rz);                                        /* norm_r = |r|_B */
206:   }
207:   else {
208:     norm_r = 0;
209:   }

211:   /* Log the residual and call any registered moitor routines */
212:   KSPLogResidualHistory(ksp, norm_r);
213:   KSPMonitor(ksp, 0, norm_r);
214:   ksp->rnorm = norm_r;

216:   /* Test for convergence */
217:   (*ksp->converged)(ksp, 0, norm_r, &ksp->reason, ksp->cnvP);
218:   if (ksp->reason) {
219:     return(0);
220:   }

222:   /* Compute the initial vectors and variables for trust-region computations */
223:   VecCopy(z, p);                    /* p = z       */

225:   if (STCG_PRECONDITIONED_DIRECTION == cg->dtype) {
226:     dMp = 0.0;
227:     norm_p = rz;
228:     norm_d = 0.0;
229:   }
230:   else {
231:     dMp = 0.0;
232:     STCG_VecDot(p, p, &norm_p);
233:     norm_d = 0.0;
234:   }

236:   /* Compute the direction */
237:   KSP_MatMult(ksp, Qmat, p, z);         /* Qp = Q * p   */
238:   STCG_VecDot(p, z, &kappa);         /* kappa = p^T Qp */

240:   if ((kappa != kappa) || (kappa && (kappa / kappa != kappa / kappa))) {
241:     ksp->reason = KSP_DIVERGED_NAN;
242:     PetscInfo1(ksp, "KSPSolve_STCG: bad matrix: kappa=%g\n", kappa);

244:     /* In this case, the matrix produced not a number or an infinite value. */
245:     /* We just take the gradient step.  Only needs to be checked once.      */
246:     if (cg->radius) {
247:       STCG_VecDot(r, r, &rz);         /* rz = r^T r   */

249:       alpha = sqrt(r2 / rz);
250:       VecAXPY(d, alpha, r);         /* d = d + alpha r */
251:       cg->norm_d = cg->radius;
252:     }
253:     return(0);
254:   }

256:   /* Begin iterating */
257:   for (i = 0; i <= maxit; i++) {
258:     ++ksp->its;

260:     /* Check for negative curvature */
261:     if (kappa <= 0.0) {
262:       ksp->reason = KSP_CONVERGED_STCG_NEG_CURVE;
263:       PetscInfo1(ksp, "KSPSolve_STCG: negative curvature: kappa=%g\n", kappa);

265:       /* In this case, the matrix is indefinite and we have encountered     */
266:       /* a direction of negative curvature.  Follow the direction to the    */
267:       /* boundary of the trust region.                                      */
268:       if (cg->radius) {
269:         alpha = (sqrt(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
270:         VecAXPY(d, alpha, p);         /* d = d + alpha p */
271:         cg->norm_d = cg->radius;
272:       }
273:       break;
274:     }
275:     alpha = rz / kappa;

277:     /* Now we can update the direction and residual.  This is perhaps not   */
278:     /* the best order to perform the breakdown checks, but the optimization */
279:     /* codes need a direction and this is the best we can do.               */

281:     /* First test if the new direction intersects the trust region. */
282:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
283:     if (cg->radius && norm_dp1 >= r2) {
284:       ksp->reason = KSP_CONVERGED_STCG_CONSTRAINED;
285:       PetscInfo1(ksp, "KSPSolve_STCG: constrained step: radius=%g\n", cg->radius);

287:       /* In this case, the matrix is  positive definite as far as we know.  */
288:       /* However, the direction does beyond the trust region.  Follow the   */
289:       /* direction to the boundary of the trust region.                     */

291:       alpha = (sqrt(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
292:       VecAXPY(d, alpha, p);         /* d = d + alpha p */
293:       cg->norm_d = cg->radius;
294:       break;
295:     }

297:     /* Update the direction and residual */
298:     VecAXPY(d, alpha, p);                 /* d = d + alpha p */
299:     VecAXPY(r, -alpha, z);                        /* r = r - alpha Qp */
300:     KSP_PCApply(ksp, r, z);

302:     if (STCG_PRECONDITIONED_DIRECTION == cg->dtype) {
303:       norm_d = norm_dp1;
304:     }
305:     else {
306:       STCG_VecDot(d, d, &norm_d);
307:     }
308:     cg->norm_d = sqrt(norm_d);

310:     /* Check that the preconditioner is positive semidefinite */
311:     rzm1 = rz;
312:     STCG_VecDot(r, z, &rz);         /* rz = r^T z   */
313:     if (rz <= 0.0) {
314:       VecNorm(r, NORM_2, &norm_r);
315:       if (rz < 0.0 || norm_r > 0.0) {
316:         ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
317:         PetscInfo1(ksp, "KSPSolve_STCG: indefinite preconditioner: rz=%g\n", rz);

319:         /* In this case, the preconditioner is indefinite.  We could follow   */
320:         /* the direction to the boundary of the trust region, but it seems    */
321:         /* best to stop at the current point.                                 */
322:         break;
323:       }
324:     }

326:     /* As far as we know, the matrix and preconditioner are positive        */
327:     /* definite.  Compute the appropriate residual depending on what the    */
328:     /* user has set.                                                        */
329:     if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
330:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
331:     }
332:     else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
333:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
334:     }
335:     else if (ksp->normtype == KSP_NATURAL_NORM) {
336:       norm_r = sqrt(rz);                                /* norm_r = |r|_B */
337:     }
338:     else {
339:       norm_r = 0;
340:     }

342:     /* Log the residual and call any registered moitor routines */
343:     KSPLogResidualHistory(ksp, norm_r);
344:     KSPMonitor(ksp, 0, norm_r);
345:     ksp->rnorm = norm_r;
346: 
347:     /* Test for convergence */
348:     (*ksp->converged)(ksp, i+1, norm_r, &ksp->reason, ksp->cnvP);
349:     if (ksp->reason) {                 /* convergence */
350:       PetscInfo2(ksp,"KSPSolve_STCG: truncated step: rnorm=%g, radius=%g\n",norm_r,cg->radius);
351:       break;
352:     }

354:     /* Update p and the norms */
355:     beta = rz / rzm1;
356:     VecAYPX(p, beta, z);                    /* p = z + beta p */

358:     if (STCG_PRECONDITIONED_DIRECTION == cg->dtype) {
359:       dMp = beta*(dMp + alpha*norm_p);
360:       norm_p = beta*(rzm1 + beta*norm_p);
361:     }
362:     else {
363:       STCG_VecDot(d, p, &dMp);
364:       STCG_VecDot(p, p, &norm_p);
365:     }

367:     /* Compute new direction */
368:     KSP_MatMult(ksp, Qmat, p, z); /* Qp = Q * p   */
369:     STCG_VecDot(p, z, &kappa);         /* kappa = p^T Qp */
370:   }

372:   if (!ksp->reason) {
373:     ksp->reason = KSP_DIVERGED_ITS;
374:   }
375:   return(0);
376: }

380: PetscErrorCode KSPSetUp_STCG(KSP ksp)
381: {

385:   /* This implementation of CG only handles left preconditioning
386:    * so generate an error otherwise.
387:    */
388:   if (ksp->pc_side == PC_RIGHT) {
389:     SETERRQ(PETSC_ERR_SUP, "No right preconditioning for KSPSTCG");
390:   }
391:   else if (ksp->pc_side == PC_SYMMETRIC) {
392:     SETERRQ(PETSC_ERR_SUP, "No symmetric preconditioning for KSPSTCG");
393:   }

395:   /* get work vectors needed by CG */
396:   KSPDefaultGetWork(ksp, 3);
397:   return(0);
398: }

402: PetscErrorCode KSPDestroy_STCG(KSP ksp)
403: {
404:   KSP_STCG *cg = (KSP_STCG *)ksp->data;

408:   KSPDefaultFreeWork(ksp);

410:   /* Free the context variable */
411:   PetscFree(cg);
412:   return(0);
413: }

418: PetscErrorCode  KSPSTCGSetRadius_STCG(KSP ksp,PetscReal radius)
419: {
420:   KSP_STCG *cg = (KSP_STCG *)ksp->data;

423:   cg->radius = radius;
424:   return(0);
425: }

429: PetscErrorCode  KSPSTCGGetNormD_STCG(KSP ksp,PetscReal *norm_d)
430: {
431:   KSP_STCG *cg = (KSP_STCG *)ksp->data;

434:   *norm_d = cg->norm_d;
435:   return(0);
436: }

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

445: PetscErrorCode KSPSetFromOptions_STCG(KSP ksp)
446: {
448:   KSP_STCG *cg = (KSP_STCG *)ksp->data;

451:   PetscOptionsHead("KSP STCG options");
452:   PetscOptionsReal("-ksp_stcg_radius", "Trust Region Radius", "KSPSTCGSetRadius", cg->radius, &cg->radius, PETSC_NULL);
453:   PetscOptionsEList("-ksp_stcg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, PETSC_NULL);
454:   PetscOptionsTail();
455:   return(0);
456: }

458: /*MC
459:      KSPSTCG -   Code to run conjugate gradient method subject to a constraint
460:          on the solution norm. This is used in Trust Region methods for
461:          nonlinear equations, SNESTR

463:    Options Database Keys:
464: .      -ksp_stcg_radius <r> - Trust Region Radius

466:    Notes: This is rarely used directly

468:    Level: developer

470: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSTCGSetRadius(), KSPSTCGGetNormD()
471: M*/

476: PetscErrorCode  KSPCreate_STCG(KSP ksp)
477: {
479:   KSP_STCG        *cg;

482:   PetscNew(KSP_STCG, &cg);
483:   PetscLogObjectMemory(ksp, sizeof(KSP_STCG));
484:   cg->radius                     = PETSC_MAX;
485:   cg->dtype                         = STCG_UNPRECONDITIONED_DIRECTION;
486:   ksp->data                      = (void *)cg;
487:   ksp->pc_side                   = PC_LEFT;

489:   /* Sets the functions that are associated with this data structure
490:    * (in C++ this is the same as defining virtual functions)
491:    */
492:   ksp->ops->setup                = KSPSetUp_STCG;
493:   ksp->ops->solve                = KSPSolve_STCG;
494:   ksp->ops->destroy              = KSPDestroy_STCG;
495:   ksp->ops->setfromoptions       = KSPSetFromOptions_STCG;
496:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
497:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;
498:   ksp->ops->view                 = 0;

500:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGSetRadius_C",
501:                                     "KSPSTCGSetRadius_STCG",
502:                                      KSPSTCGSetRadius_STCG);
503:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPSTCGGetNormD_C",
504:                                     "KSPSTCGGetNormD_STCG",
505:                                      KSPSTCGGetNormD_STCG);
506:   return(0);
507: }