Actual source code: qcg.c


  2: #include <../src/ksp/ksp/impls/qcg/qcgimpl.h>

  4: static PetscErrorCode KSPQCGQuadraticRoots(Vec, Vec, PetscReal, PetscReal *, PetscReal *);

  6: /*@
  7:     KSPQCGSetTrustRegionRadius - Sets the radius of the trust region for `KSPQCG`

  9:     Logically Collective

 11:     Input Parameters:
 12: +   ksp   - the iterative context
 13: -   delta - the trust region radius (Infinity is the default)

 15:     Options Database Key:
 16: .   -ksp_qcg_trustregionradius <delta> - trust region radius

 18:     Level: advanced

 20: @*/
 21: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
 22: {
 25:   PetscTryMethod(ksp, "KSPQCGSetTrustRegionRadius_C", (KSP, PetscReal), (ksp, delta));
 26:   return 0;
 27: }

 29: /*@
 30:     KSPQCGGetTrialStepNorm - Gets the norm of a trial step vector in `KSPQCG`.  The WCG step may be
 31:     constrained, so this is not necessarily the length of the ultimate step taken in `KSPQCG`.

 33:     Not Collective

 35:     Input Parameter:
 36: .   ksp - the iterative context

 38:     Output Parameter:
 39: .   tsnorm - the norm

 41:     Level: advanced
 42: @*/
 43: PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
 44: {
 46:   PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
 47:   return 0;
 48: }

 50: /*@
 51:     KSPQCGGetQuadratic - Gets the value of the quadratic function, evaluated at the new iterate:

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

 55:     which satisfies the Euclidean Norm trust region constraint

 57:        || D * s || <= delta,

 59:     where

 61:      delta is the trust region radius,
 62:      g is the gradient vector, and
 63:      H is Hessian matrix,
 64:      D is a scaling matrix.

 66:     Collective

 68:     Input Parameter:
 69: .   ksp - the iterative context

 71:     Output Parameter:
 72: .   quadratic - the quadratic function evaluated at the new iterate

 74:     Level: advanced

 76: .seealso: [](chapter_ksp), `KSPQCG`
 77: @*/
 78: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
 79: {
 81:   PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
 82:   return 0;
 83: }

 85: PetscErrorCode KSPSolve_QCG(KSP ksp)
 86: {
 87:   /*
 88:    Correpondence with documentation above:
 89:       B = g = gradient,
 90:       X = s = step
 91:    Note:  This is not coded correctly for complex arithmetic!
 92:  */

 94:   KSP_QCG    *pcgP = (KSP_QCG *)ksp->data;
 95:   Mat         Amat, Pmat;
 96:   Vec         W, WA, WA2, R, P, ASP, BS, X, B;
 97:   PetscScalar scal, beta, rntrn, step;
 98:   PetscReal   q1, q2, xnorm, step1, step2, rnrm = 0.0, btx, xtax;
 99:   PetscReal   ptasp, rtr, wtasp, bstp;
100:   PetscReal   dzero = 0.0, bsnrm = 0.0;
101:   PetscInt    i, maxit;
102:   PC          pc = ksp->pc;
103:   PetscBool   diagonalscale;

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

109:   ksp->its = 0;
110:   maxit    = ksp->max_it;
111:   WA       = ksp->work[0];
112:   R        = ksp->work[1];
113:   P        = ksp->work[2];
114:   ASP      = ksp->work[3];
115:   BS       = ksp->work[4];
116:   W        = ksp->work[5];
117:   WA2      = ksp->work[6];
118:   X        = ksp->vec_sol;
119:   B        = ksp->vec_rhs;


123:   /* Initialize variables */
124:   VecSet(W, 0.0); /* W = 0 */
125:   VecSet(X, 0.0); /* X = 0 */
126:   PCGetOperators(pc, &Amat, &Pmat);

128:   /* Compute:  BS = D^{-1} B */
129:   PCApplySymmetricLeft(pc, B, BS);

131:   if (ksp->normtype != KSP_NORM_NONE) {
132:     VecNorm(BS, NORM_2, &bsnrm);
133:     KSPCheckNorm(ksp, bsnrm);
134:   }
135:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
136:   ksp->its   = 0;
137:   ksp->rnorm = bsnrm;
138:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
139:   KSPLogResidualHistory(ksp, bsnrm);
140:   KSPMonitor(ksp, 0, bsnrm);
141:   (*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP);
142:   if (ksp->reason) return 0;

144:   /* Compute the initial scaled direction and scaled residual */
145:   VecCopy(BS, R);
146:   VecScale(R, -1.0);
147:   VecCopy(R, P);
148:   VecDotRealPart(R, R, &rtr);

150:   for (i = 0; i <= maxit; i++) {
151:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
152:     ksp->its++;
153:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

155:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
156:     PCApplySymmetricRight(pc, P, WA);
157:     KSP_MatMult(ksp, Amat, WA, WA2);
158:     PCApplySymmetricLeft(pc, WA2, ASP);

160:     /* Check for negative curvature */
161:     VecDotRealPart(P, ASP, &ptasp);
162:     if (ptasp <= dzero) {
163:       /* Scaled negative curvature direction:  Compute a step so that
164:         ||w + step*p|| = delta and QS(w + step*p) is least */

166:       if (!i) {
167:         VecCopy(P, X);
168:         VecNorm(X, NORM_2, &xnorm);
169:         KSPCheckNorm(ksp, xnorm);
170:         scal = pcgP->delta / xnorm;
171:         VecScale(X, scal);
172:       } else {
173:         /* Compute roots of quadratic */
174:         KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2);
175:         VecDotRealPart(W, ASP, &wtasp);
176:         VecDotRealPart(BS, P, &bstp);
177:         VecCopy(W, X);
178:         q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
179:         q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
180:         if (q1 <= q2) {
181:           VecAXPY(X, step1, P);
182:         } else {
183:           VecAXPY(X, step2, P);
184:         }
185:       }
186:       pcgP->ltsnrm = pcgP->delta;                /* convergence in direction of */
187:       ksp->reason  = KSP_CONVERGED_CG_NEG_CURVE; /* negative curvature */
188:       if (!i) {
189:         PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta);
190:       } else {
191:         PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta);
192:       }

194:     } else {
195:       /* Compute step along p */
196:       step = rtr / ptasp;
197:       VecCopy(W, X);       /*  x = w  */
198:       VecAXPY(X, step, P); /*  x <- step*p + x  */
199:       VecNorm(X, NORM_2, &pcgP->ltsnrm);
200:       KSPCheckNorm(ksp, pcgP->ltsnrm);

202:       if (pcgP->ltsnrm > pcgP->delta) {
203:         /* Since the trial iterate is outside the trust region,
204:             evaluate a constrained step along p so that
205:                     ||w + step*p|| = delta
206:           The positive step is always better in this case. */
207:         if (!i) {
208:           scal = pcgP->delta / pcgP->ltsnrm;
209:           VecScale(X, scal);
210:         } else {
211:           /* Compute roots of quadratic */
212:           KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2);
213:           VecCopy(W, X);
214:           VecAXPY(X, step1, P); /*  x <- step1*p + x  */
215:         }
216:         pcgP->ltsnrm = pcgP->delta;
217:         ksp->reason  = KSP_CONVERGED_CG_CONSTRAINED; /* convergence along constrained step */
218:         if (!i) {
219:           PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta);
220:         } else {
221:           PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta);
222:         }

224:       } else {
225:         /* Evaluate the current step */
226:         VecCopy(X, W);          /* update interior iterate */
227:         VecAXPY(R, -step, ASP); /* r <- -step*asp + r */
228:         if (ksp->normtype != KSP_NORM_NONE) {
229:           VecNorm(R, NORM_2, &rnrm);
230:           KSPCheckNorm(ksp, rnrm);
231:         }
232:         PetscObjectSAWsTakeAccess((PetscObject)ksp);
233:         ksp->rnorm = rnrm;
234:         PetscObjectSAWsGrantAccess((PetscObject)ksp);
235:         KSPLogResidualHistory(ksp, rnrm);
236:         KSPMonitor(ksp, i + 1, rnrm);
237:         (*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP);
238:         if (ksp->reason) { /* convergence for */
239:           PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta);
240:         }
241:       }
242:     }
243:     if (ksp->reason) break; /* Convergence has been attained */
244:     else { /* Compute a new AS-orthogonal direction */ VecDot(R, R, &rntrn);
245:       beta = rntrn / rtr;
246:       VecAYPX(P, beta, R); /*  p <- r + beta*p  */
247:       rtr = PetscRealPart(rntrn);
248:     }
249:   }
250:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

252:   /* Unscale x */
253:   VecCopy(X, WA2);
254:   PCApplySymmetricRight(pc, WA2, X);

256:   KSP_MatMult(ksp, Amat, X, WA);
257:   VecDotRealPart(B, X, &btx);
258:   VecDotRealPart(X, WA, &xtax);

260:   pcgP->quadratic = btx + .5 * xtax;
261:   return 0;
262: }

264: PetscErrorCode KSPSetUp_QCG(KSP ksp)
265: {
266:   /* Get work vectors from user code */
267:   KSPSetWorkVecs(ksp, 7);
268:   return 0;
269: }

271: PetscErrorCode KSPDestroy_QCG(KSP ksp)
272: {
273:   PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL);
274:   PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL);
275:   PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL);
276:   KSPDestroyDefault(ksp);
277:   return 0;
278: }

280: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
281: {
282:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

284:   cgP->delta = delta;
285:   return 0;
286: }

288: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
289: {
290:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

292:   *ltsnrm = cgP->ltsnrm;
293:   return 0;
294: }

296: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
297: {
298:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

300:   *quadratic = cgP->quadratic;
301:   return 0;
302: }

304: PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
305: {
306:   PetscReal delta;
307:   KSP_QCG  *cgP = (KSP_QCG *)ksp->data;
308:   PetscBool flg;

310:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
311:   PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg);
312:   if (flg) KSPQCGSetTrustRegionRadius(ksp, delta);
313:   PetscOptionsHeadEnd();
314:   return 0;
315: }

317: /*MC
318:      KSPQCG -   Code to run conjugate gradient method subject to a constraint on the solution norm.

320:    Options Database Key:
321: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

323:    Level: developer

325:    Notes:
326:     This is rarely used directly, ir is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`

328:     Uses preconditioned conjugate gradient to compute
329:       an approximate minimizer of the quadratic function

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

333:    subject to the Euclidean norm trust region constraint

335:             || D * s || <= delta,

337:    where

339:      delta is the trust region radius,
340:      g is the gradient vector, and
341:      H is Hessian matrix,
342:      D is a scaling matrix.

344:    `KSPConvergedReason` may be
345: .vb
346:    KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
347:    KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
348: .ve
349:    and other `KSP` converged/diverged reasons

351:   Notes:
352:   Currently we allow symmetric preconditioning with the following scaling matrices:
353: .vb
354:       `PCNONE`:   D = Identity matrix
355:       `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
356:       `PCICC`:    D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
357: .ve

359:   References:
360: . * - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
361:    SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).

363: .seealso: [](chapter_ksp), 'KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`
364:           `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
365: M*/

367: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
368: {
369:   KSP_QCG *cgP;

371:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3);
372:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1);
373:   PetscNew(&cgP);

375:   ksp->data                = (void *)cgP;
376:   ksp->ops->setup          = KSPSetUp_QCG;
377:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
378:   ksp->ops->solve          = KSPSolve_QCG;
379:   ksp->ops->destroy        = KSPDestroy_QCG;
380:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
381:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
382:   ksp->ops->view           = NULL;

384:   PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG);
385:   PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG);
386:   PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG);
387:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
388:   return 0;
389: }

391: /* ---------------------------------------------------------- */
392: /*
393:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
394:          ||s + step*p|| - delta = 0
395:    such that step1 >= 0 >= step2.
396:    where
397:       delta:
398:         On entry delta must contain scalar delta.
399:         On exit delta is unchanged.
400:       step1:
401:         On entry step1 need not be specified.
402:         On exit step1 contains the non-negative root.
403:       step2:
404:         On entry step2 need not be specified.
405:         On exit step2 contains the non-positive root.
406:    C code is translated from the Fortran version of the MINPACK-2 Project,
407:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
408: */
409: static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
410: {
411:   PetscReal dsq, ptp, pts, rad, sts;

413:   VecDotRealPart(p, s, &pts);
414:   VecDotRealPart(p, p, &ptp);
415:   VecDotRealPart(s, s, &sts);
416:   dsq = delta * delta;
417:   rad = PetscSqrtReal((pts * pts) - ptp * (sts - dsq));
418:   if (pts > 0.0) {
419:     *step2 = -(pts + rad) / ptp;
420:     *step1 = (sts - dsq) / (ptp * *step2);
421:   } else {
422:     *step1 = -(pts - rad) / ptp;
423:     *step2 = (sts - dsq) / (ptp * *step1);
424:   }
425:   return 0;
426: }