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.

  9:     Logically Collective on ksp

 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.  The WCG step may be
 31:     constrained, so this is not necessarily the length of the ultimate step taken in QCG.

 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 Euclidian 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 on ksp

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

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

 74:     Level: advanced
 75: @*/
 76: PetscErrorCode  KSPQCGGetQuadratic(KSP ksp,PetscReal *quadratic)
 77: {
 79:   PetscUseMethod(ksp,"KSPQCGGetQuadratic_C",(KSP,PetscReal*),(ksp,quadratic));
 80:   return 0;
 81: }

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

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

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

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

121:   KSPGetPCSide(ksp,&side);

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

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

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

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

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

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

161:     /* Check for negative curvature */
162:     VecDotRealPart(P,ASP,&ptasp);
163:     if (ptasp <= dzero) {

165:       /* Scaled negative curvature direction:  Compute a step so that
166:         ||w + step*p|| = delta and QS(w + step*p) is least */

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

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

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

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

255:   /* Unscale x */
256:   VecCopy(X,WA2);
257:   PCApplySymmetricRight(pc,WA2,X);

259:   KSP_MatMult(ksp,Amat,X,WA);
260:   VecDotRealPart(B,X,&btx);
261:   VecDotRealPart(X,WA,&xtax);

263:   pcgP->quadratic = btx + .5*xtax;
264:   return 0;
265: }

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

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

283: static PetscErrorCode  KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
284: {
285:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

287:   cgP->delta = delta;
288:   return 0;
289: }

291: static PetscErrorCode  KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
292: {
293:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

295:   *ltsnrm = cgP->ltsnrm;
296:   return 0;
297: }

299: static PetscErrorCode  KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
300: {
301:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

303:   *quadratic = cgP->quadratic;
304:   return 0;
305: }

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

313:   PetscOptionsHead(PetscOptionsObject,"KSP QCG Options");
314:   PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
315:   if (flg) KSPQCGSetTrustRegionRadius(ksp,delta);
316:   PetscOptionsTail();
317:   return 0;
318: }

320: /*MC
321:      KSPQCG -   Code to run conjugate gradient method subject to a constraint
322:          on the solution norm. This is used in Trust Region methods for nonlinear equations, SNESNEWTONTR

324:    Options Database Keys:
325: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

327:    Notes:
328:     This is rarely used directly

330:    Level: developer

332:   Notes:
333:     Use preconditioned conjugate gradient to compute
334:       an approximate minimizer of the quadratic function

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

338:    subject to the Euclidean norm trust region constraint

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

342:    where

344:      delta is the trust region radius,
345:      g is the gradient vector, and
346:      H is Hessian matrix,
347:      D is a scaling matrix.

349:    KSPConvergedReason may be
350: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
351: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
352: $  other KSP converged/diverged reasons

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

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

365: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
366:            KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
367: M*/

369: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
370: {
371:   KSP_QCG        *cgP;

373:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,3);
374:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_SYMMETRIC,1);
375:   PetscNewLog(ksp,&cgP);

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

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

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

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