Actual source code: qcg.c

petsc-3.9.2 2018-05-20
Report Typos and Errors
```
2:  #include <../src/ksp/ksp/impls/qcg/qcgimpl.h>

6: /*@

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:

20: .keywords: KSP, QCG, set, trust region radius
21: @*/
23: {

28:   if (delta < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
30:   return(0);
31: }

33: /*@
34:     KSPQCGGetTrialStepNorm - Gets the norm of a trial step vector.  The WCG step may be
35:     constrained, so this is not necessarily the length of the ultimate step taken in QCG.

37:     Not Collective

39:     Input Parameter:
40: .   ksp - the iterative context

42:     Output Parameter:
43: .   tsnorm - the norm

46: @*/
47: PetscErrorCode  KSPQCGGetTrialStepNorm(KSP ksp,PetscReal *tsnorm)
48: {

53:   PetscUseMethod(ksp,"KSPQCGGetTrialStepNorm_C",(KSP,PetscReal*),(ksp,tsnorm));
54:   return(0);
55: }

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

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

62:     which satisfies the Euclidian Norm trust region constraint

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

66:     where

68:      delta is the trust region radius,
69:      g is the gradient vector, and
70:      H is Hessian matrix,
71:      D is a scaling matrix.

73:     Collective on KSP

75:     Input Parameter:
76: .   ksp - the iterative context

78:     Output Parameter:
79: .   quadratic - the quadratic function evaluated at the new iterate

82: @*/
84: {

90:   return(0);
91: }

94: PetscErrorCode KSPSolve_QCG(KSP ksp)
95: {
96: /*
97:    Correpondence with documentation above:
98:       B = g = gradient,
99:       X = s = step
100:    Note:  This is not coded correctly for complex arithmetic!
101:  */

103:   KSP_QCG        *pcgP = (KSP_QCG*)ksp->data;
104:   Mat            Amat,Pmat;
105:   Vec            W,WA,WA2,R,P,ASP,BS,X,B;
106:   PetscScalar    scal,beta,rntrn,step;
107:   PetscReal      q1,q2,xnorm,step1,step2,rnrm,btx,xtax;
108:   PetscReal      ptasp,rtr,wtasp,bstp;
109:   PetscReal      dzero = 0.0,bsnrm;
111:   PetscInt       i,maxit;
112:   PC             pc = ksp->pc;
113:   PCSide         side;
114:   PetscBool      diagonalscale;

117:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
118:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
119:   if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Currently does not support transpose solve");

121:   ksp->its = 0;
122:   maxit    = ksp->max_it;
123:   WA       = ksp->work[0];
124:   R        = ksp->work[1];
125:   P        = ksp->work[2];
126:   ASP      = ksp->work[3];
127:   BS       = ksp->work[4];
128:   W        = ksp->work[5];
129:   WA2      = ksp->work[6];
130:   X        = ksp->vec_sol;
131:   B        = ksp->vec_rhs;

133:   if (pcgP->delta <= dzero) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Input error: delta <= 0");
134:   KSPGetPCSide(ksp,&side);
135:   if (side != PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requires symmetric preconditioner!");

137:   /* Initialize variables */
138:   VecSet(W,0.0);  /* W = 0 */
139:   VecSet(X,0.0);  /* X = 0 */
140:   PCGetOperators(pc,&Amat,&Pmat);

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

145:   VecNorm(BS,NORM_2,&bsnrm);
146:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
147:   ksp->its   = 0;
148:   ksp->rnorm = bsnrm;
149:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
150:   KSPLogResidualHistory(ksp,bsnrm);
151:   KSPMonitor(ksp,0,bsnrm);
152:   (*ksp->converged)(ksp,0,bsnrm,&ksp->reason,ksp->cnvP);
153:   if (ksp->reason) return(0);

155:   /* Compute the initial scaled direction and scaled residual */
156:   VecCopy(BS,R);
157:   VecScale(R,-1.0);
158:   VecCopy(R,P);
159:   VecDotRealPart(R,R,&rtr);

161:   for (i=0; i<=maxit; i++) {
162:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
163:     ksp->its++;
164:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

166:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
167:     PCApplySymmetricRight(pc,P,WA);
168:     KSP_MatMult(ksp,Amat,WA,WA2);
169:     PCApplySymmetricLeft(pc,WA2,ASP);

171:     /* Check for negative curvature */
172:     VecDotRealPart(P,ASP,&ptasp);
173:     if (ptasp <= dzero) {

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

178:       if (!i) {
179:         VecCopy(P,X);
180:         VecNorm(X,NORM_2,&xnorm);
181:         scal = pcgP->delta / xnorm;
182:         VecScale(X,scal);
183:       } else {
184:         /* Compute roots of quadratic */
186:         VecDotRealPart(W,ASP,&wtasp);
187:         VecDotRealPart(BS,P,&bstp);
188:         VecCopy(W,X);
189:         q1   = step1*(bstp + wtasp + .5*step1*ptasp);
190:         q2   = step2*(bstp + wtasp + .5*step2*ptasp);
191:         if (q1 <= q2) {
192:           VecAXPY(X,step1,P);
193:         } else {
194:           VecAXPY(X,step2,P);
195:         }
196:       }
197:       pcgP->ltsnrm = pcgP->delta;                       /* convergence in direction of */
198:       ksp->reason  = KSP_CONVERGED_CG_NEG_CURVE;  /* negative curvature */
199:       if (!i) {
200:         PetscInfo1(ksp,"negative curvature: delta=%g\n",(double)pcgP->delta);
201:       } else {
202:         PetscInfo3(ksp,"negative curvature: step1=%g, step2=%g, delta=%g\n",(double)step1,(double)step2,(double)pcgP->delta);
203:       }

205:     } else {
206:       /* Compute step along p */
207:       step = rtr/ptasp;
208:       VecCopy(W,X);        /*  x = w  */
209:       VecAXPY(X,step,P);   /*  x <- step*p + x  */
210:       VecNorm(X,NORM_2,&pcgP->ltsnrm);

212:       if (pcgP->ltsnrm > pcgP->delta) {
213:         /* Since the trial iterate is outside the trust region,
214:             evaluate a constrained step along p so that
215:                     ||w + step*p|| = delta
216:           The positive step is always better in this case. */
217:         if (!i) {
218:           scal = pcgP->delta / pcgP->ltsnrm;
219:           VecScale(X,scal);
220:         } else {
221:           /* Compute roots of quadratic */
223:           VecCopy(W,X);
224:           VecAXPY(X,step1,P);  /*  x <- step1*p + x  */
225:         }
226:         pcgP->ltsnrm = pcgP->delta;
227:         ksp->reason  = KSP_CONVERGED_CG_CONSTRAINED; /* convergence along constrained step */
228:         if (!i) {
229:           PetscInfo1(ksp,"constrained step: delta=%g\n",(double)pcgP->delta);
230:         } else {
231:           PetscInfo3(ksp,"constrained step: step1=%g, step2=%g, delta=%g\n",(double)step1,(double)step2,(double)pcgP->delta);
232:         }

234:       } else {
235:         /* Evaluate the current step */
236:         VecCopy(X,W);  /* update interior iterate */
237:         VecAXPY(R,-step,ASP); /* r <- -step*asp + r */
238:         VecNorm(R,NORM_2,&rnrm);

240:         PetscObjectSAWsTakeAccess((PetscObject)ksp);
241:         ksp->rnorm = rnrm;
242:         PetscObjectSAWsGrantAccess((PetscObject)ksp);
243:         KSPLogResidualHistory(ksp,rnrm);
244:         KSPMonitor(ksp,i+1,rnrm);
245:         (*ksp->converged)(ksp,i+1,rnrm,&ksp->reason,ksp->cnvP);
246:         if (ksp->reason) {                 /* convergence for */
247:           PetscInfo3(ksp,"truncated step: step=%g, rnrm=%g, delta=%g\n",(double)PetscRealPart(step),(double)rnrm,(double)pcgP->delta);
248:         }
249:       }
250:     }
251:     if (ksp->reason) break;  /* Convergence has been attained */
252:     else {                   /* Compute a new AS-orthogonal direction */
253:       VecDot(R,R,&rntrn);
254:       beta = rntrn/rtr;
255:       VecAYPX(P,beta,R);  /*  p <- r + beta*p  */
256:       rtr  = PetscRealPart(rntrn);
257:     }
258:   }
259:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

261:   /* Unscale x */
262:   VecCopy(X,WA2);
263:   PCApplySymmetricRight(pc,WA2,X);

265:   KSP_MatMult(ksp,Amat,X,WA);
266:   VecDotRealPart(B,X,&btx);
267:   VecDotRealPart(X,WA,&xtax);

269:   pcgP->quadratic = btx + .5*xtax;
270:   return(0);
271: }

273: PetscErrorCode KSPSetUp_QCG(KSP ksp)
274: {

278:   /* Get work vectors from user code */
279:   KSPSetWorkVecs(ksp,7);
280:   return(0);
281: }

283: PetscErrorCode KSPDestroy_QCG(KSP ksp)
284: {

289:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",NULL);
291:   KSPDestroyDefault(ksp);
292:   return(0);
293: }

295: static PetscErrorCode  KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
296: {
297:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

300:   cgP->delta = delta;
301:   return(0);
302: }

304: static PetscErrorCode  KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
305: {
306:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

309:   *ltsnrm = cgP->ltsnrm;
310:   return(0);
311: }

314: {
315:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

319:   return(0);
320: }

322: PetscErrorCode KSPSetFromOptions_QCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
323: {
325:   PetscReal      delta;
326:   KSP_QCG        *cgP = (KSP_QCG*)ksp->data;
327:   PetscBool      flg;

332:   if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
333:   PetscOptionsTail();
334:   return(0);
335: }

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

341:    Options Database Keys:

344:    Notes: This is rarely used directly

346:    Level: developer

348:   Notes:  Use preconditioned conjugate gradient to compute
349:       an approximate minimizer of the quadratic function

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

353:    subject to the Euclidean norm trust region constraint

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

357:    where

359:      delta is the trust region radius,
360:      g is the gradient vector, and
361:      H is Hessian matrix,
362:      D is a scaling matrix.

364:    KSPConvergedReason may be
365: \$  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
366: \$  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
367: \$  other KSP converged/diverged reasons

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

377:   References:
378: .  1. - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
379:    SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).

381: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
383: M*/

385: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
386: {
388:   KSP_QCG        *cgP;

391:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,3);
392:   PetscNewLog(ksp,&cgP);

394:   ksp->data                = (void*)cgP;
395:   ksp->ops->setup          = KSPSetUp_QCG;
396:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
397:   ksp->ops->solve          = KSPSolve_QCG;
398:   ksp->ops->destroy        = KSPDestroy_QCG;
399:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
400:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
401:   ksp->ops->setfromoptions = 0;
402:   ksp->ops->view           = 0;

405:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",KSPQCGGetTrialStepNorm_QCG);
407:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
408:   return(0);
409: }

411: /* ---------------------------------------------------------- */
412: /*
414:          ||s + step*p|| - delta = 0
415:    such that step1 >= 0 >= step2.
416:    where
417:       delta:
418:         On entry delta must contain scalar delta.
419:         On exit delta is unchanged.
420:       step1:
421:         On entry step1 need not be specified.
422:         On exit step1 contains the non-negative root.
423:       step2:
424:         On entry step2 need not be specified.
425:         On exit step2 contains the non-positive root.
426:    C code is translated from the Fortran version of the MINPACK-2 Project,
427:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
428: */
429: static PetscErrorCode KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal delta,PetscReal *step1,PetscReal *step2)
430: {

435:   VecDotRealPart(p,s,&pts);
436:   VecDotRealPart(p,p,&ptp);
437:   VecDotRealPart(s,s,&sts);
438:   dsq  = delta*delta;
439:   rad  = PetscSqrtReal((pts*pts) - ptp*(sts - dsq));
440:   if (pts > 0.0) {
441:     *step2 = -(pts + rad)/ptp;
442:     *step1 = (sts - dsq)/(ptp * *step2);
443:   } else {
444:     *step1 = -(pts - rad)/ptp;
445:     *step2 = (sts - dsq)/(ptp * *step1);
446:   }
447:   return(0);
448: }
```