Actual source code: qcg.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <petsc-private/kspimpl.h>             /*I "petscksp.h" I*/
  3: #include <../src/ksp/ksp/impls/qcg/qcgimpl.h>

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

  9: /*@
 10:     KSPQCGSetTrustRegionRadius - Sets the radius of the trust region.

 12:     Logically Collective on KSP

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

 18:     Options Database Key:
 19: .   -ksp_qcg_trustregionradius <delta>

 21:     Level: advanced

 23: .keywords: KSP, QCG, set, trust region radius
 24: @*/
 25: PetscErrorCode  KSPQCGSetTrustRegionRadius(KSP ksp,PetscReal delta)
 26: {

 31:   if (delta < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
 32:   PetscTryMethod(ksp,"KSPQCGSetTrustRegionRadius_C",(KSP,PetscReal),(ksp,delta));
 33:   return(0);
 34: }

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

 42:     Not Collective

 44:     Input Parameter:
 45: .   ksp - the iterative context

 47:     Output Parameter:
 48: .   tsnorm - the norm

 50:     Level: advanced
 51: @*/
 52: PetscErrorCode  KSPQCGGetTrialStepNorm(KSP ksp,PetscReal *tsnorm)
 53: {

 58:   PetscUseMethod(ksp,"KSPQCGGetTrialStepNorm_C",(KSP,PetscReal*),(ksp,tsnorm));
 59:   return(0);
 60: }

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

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

 69:     which satisfies the Euclidian Norm trust region constraint

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

 73:     where

 75:      delta is the trust region radius,
 76:      g is the gradient vector, and
 77:      H is Hessian matrix,
 78:      D is a scaling matrix.

 80:     Collective on KSP

 82:     Input Parameter:
 83: .   ksp - the iterative context

 85:     Output Parameter:
 86: .   quadratic - the quadratic function evaluated at the new iterate

 88:     Level: advanced
 89: @*/
 90: PetscErrorCode  KSPQCGGetQuadratic(KSP ksp,PetscReal *quadratic)
 91: {

 96:   PetscUseMethod(ksp,"KSPQCGGetQuadratic_C",(KSP,PetscReal*),(ksp,quadratic));
 97:   return(0);
 98: }


103: PetscErrorCode KSPSolve_QCG(KSP ksp)
104: {
105: /*
106:    Correpondence with documentation above:
107:       B = g = gradient,
108:       X = s = step
109:    Note:  This is not coded correctly for complex arithmetic!
110:  */

112:   KSP_QCG        *pcgP = (KSP_QCG*)ksp->data;
113:   Mat            Amat,Pmat;
114:   Vec            W,WA,WA2,R,P,ASP,BS,X,B;
115:   PetscScalar    scal,beta,rntrn,step;
116:   PetscReal      q1,q2,xnorm,step1,step2,rnrm,btx,xtax;
117:   PetscReal      ptasp,rtr,wtasp,bstp;
118:   PetscReal      dzero = 0.0,bsnrm;
120:   PetscInt       i,maxit;
121:   PC             pc = ksp->pc;
122:   PCSide         side;
123:   PetscBool      diagonalscale;

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

130:   ksp->its = 0;
131:   maxit    = ksp->max_it;
132:   WA       = ksp->work[0];
133:   R        = ksp->work[1];
134:   P        = ksp->work[2];
135:   ASP      = ksp->work[3];
136:   BS       = ksp->work[4];
137:   W        = ksp->work[5];
138:   WA2      = ksp->work[6];
139:   X        = ksp->vec_sol;
140:   B        = ksp->vec_rhs;

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

146:   /* Initialize variables */
147:   VecSet(W,0.0);  /* W = 0 */
148:   VecSet(X,0.0);  /* X = 0 */
149:   PCGetOperators(pc,&Amat,&Pmat);

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

154:   VecNorm(BS,NORM_2,&bsnrm);
155:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
156:   ksp->its   = 0;
157:   ksp->rnorm = bsnrm;
158:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
159:   KSPLogResidualHistory(ksp,bsnrm);
160:   KSPMonitor(ksp,0,bsnrm);
161:   (*ksp->converged)(ksp,0,bsnrm,&ksp->reason,ksp->cnvP);
162:   if (ksp->reason) return(0);

164:   /* Compute the initial scaled direction and scaled residual */
165:   VecCopy(BS,R);
166:   VecScale(R,-1.0);
167:   VecCopy(R,P);
168:   VecDotRealPart(R,R,&rtr);

170:   for (i=0; i<=maxit; i++) {
171:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
172:     ksp->its++;
173:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

175:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
176:     PCApplySymmetricRight(pc,P,WA);
177:     MatMult(Amat,WA,WA2);
178:     PCApplySymmetricLeft(pc,WA2,ASP);

180:     /* Check for negative curvature */
181:     VecDotRealPart(P,ASP,&ptasp);
182:     if (ptasp <= dzero) {

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

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

214:     } else {
215:       /* Compute step along p */
216:       step = rtr/ptasp;
217:       VecCopy(W,X);        /*  x = w  */
218:       VecAXPY(X,step,P);   /*  x <- step*p + x  */
219:       VecNorm(X,NORM_2,&pcgP->ltsnrm);

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

243:       } else {
244:         /* Evaluate the current step */
245:         VecCopy(X,W);  /* update interior iterate */
246:         VecAXPY(R,-step,ASP); /* r <- -step*asp + r */
247:         VecNorm(R,NORM_2,&rnrm);

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

270:   /* Unscale x */
271:   VecCopy(X,WA2);
272:   PCApplySymmetricRight(pc,WA2,X);

274:   MatMult(Amat,X,WA);
275:   VecDotRealPart(B,X,&btx);
276:   VecDotRealPart(X,WA,&xtax);

278:   pcgP->quadratic = btx + .5*xtax;
279:   return(0);
280: }

284: PetscErrorCode KSPSetUp_QCG(KSP ksp)
285: {

289:   /* Get work vectors from user code */
290:   KSPSetWorkVecs(ksp,7);
291:   return(0);
292: }

296: PetscErrorCode KSPDestroy_QCG(KSP ksp)
297: {

301:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",NULL);
302:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",NULL);
303:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",NULL);
304:   KSPDestroyDefault(ksp);
305:   return(0);
306: }

310: static PetscErrorCode  KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
311: {
312:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

315:   cgP->delta = delta;
316:   return(0);
317: }

321: static PetscErrorCode  KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
322: {
323:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

326:   *ltsnrm = cgP->ltsnrm;
327:   return(0);
328: }

332: static PetscErrorCode  KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
333: {
334:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

337:   *quadratic = cgP->quadratic;
338:   return(0);
339: }

343: PetscErrorCode KSPSetFromOptions_QCG(KSP ksp)
344: {
346:   PetscReal      delta;
347:   KSP_QCG        *cgP = (KSP_QCG*)ksp->data;
348:   PetscBool      flg;

351:   PetscOptionsHead("KSP QCG Options");
352:   PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
353:   if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
354:   PetscOptionsTail();
355:   return(0);
356: }

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

362:    Options Database Keys:
363: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

365:    Notes: This is rarely used directly

367:    Level: developer

369:   Notes:  Use preconditioned conjugate gradient to compute
370:       an approximate minimizer of the quadratic function

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

374:    subject to the Euclidean norm trust region constraint

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

378:    where

380:      delta is the trust region radius,
381:      g is the gradient vector, and
382:      H is Hessian matrix,
383:      D is a scaling matrix.

385:    KSPConvergedReason may be
386: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
387: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
388: $  other KSP converged/diverged reasons


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

398:   References:
399:    The Conjugate Gradient Method and Trust Regions in Large Scale Optimization, Trond Steihaug
400:    SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983), pp. 626-637

402: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
403:            KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
404: M*/

408: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
409: {
411:   KSP_QCG        *cgP;

414:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,3);
415:   PetscNewLog(ksp,&cgP);

417:   ksp->data                = (void*)cgP;
418:   ksp->ops->setup          = KSPSetUp_QCG;
419:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
420:   ksp->ops->solve          = KSPSolve_QCG;
421:   ksp->ops->destroy        = KSPDestroy_QCG;
422:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
423:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
424:   ksp->ops->setfromoptions = 0;
425:   ksp->ops->view           = 0;

427:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",KSPQCGGetQuadratic_QCG);
428:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",KSPQCGGetTrialStepNorm_QCG);
429:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",KSPQCGSetTrustRegionRadius_QCG);
430:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
431:   return(0);
432: }

434: /* ---------------------------------------------------------- */
437: /*
438:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
439:          ||s + step*p|| - delta = 0
440:    such that step1 >= 0 >= step2.
441:    where
442:       delta:
443:         On entry delta must contain scalar delta.
444:         On exit delta is unchanged.
445:       step1:
446:         On entry step1 need not be specified.
447:         On exit step1 contains the non-negative root.
448:       step2:
449:         On entry step2 need not be specified.
450:         On exit step2 contains the non-positive root.
451:    C code is translated from the Fortran version of the MINPACK-2 Project,
452:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
453: */
454: static PetscErrorCode KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal delta,PetscReal *step1,PetscReal *step2)
455: {
456:   PetscReal      dsq,ptp,pts,rad,sts;

460:   VecDotRealPart(p,s,&pts);
461:   VecDotRealPart(p,p,&ptp);
462:   VecDotRealPart(s,s,&sts);
463:   dsq  = delta*delta;
464:   rad  = PetscSqrtReal((pts*pts) - ptp*(sts - dsq));
465:   if (pts > 0.0) {
466:     *step2 = -(pts + rad)/ptp;
467:     *step1 = (sts - dsq)/(ptp * *step2);
468:   } else {
469:     *step1 = -(pts - rad)/ptp;
470:     *step2 = (sts - dsq)/(ptp * *step1);
471:   }
472:   return(0);
473: }