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>

 18:     Level: advanced

 20: @*/
 21: PetscErrorCode  KSPQCGSetTrustRegionRadius(KSP ksp,PetscReal delta)
 22: {

 27:   if (delta < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
 28:   PetscTryMethod(ksp,"KSPQCGSetTrustRegionRadius_C",(KSP,PetscReal),(ksp,delta));
 29:   return(0);
 30: }

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

 36:     Not Collective

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

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

 44:     Level: advanced
 45: @*/
 46: PetscErrorCode  KSPQCGGetTrialStepNorm(KSP ksp,PetscReal *tsnorm)
 47: {

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

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

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

 61:     which satisfies the Euclidian Norm trust region constraint

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

 65:     where

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

 72:     Collective on ksp

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

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

 80:     Level: advanced
 81: @*/
 82: PetscErrorCode  KSPQCGGetQuadratic(KSP ksp,PetscReal *quadratic)
 83: {

 88:   PetscUseMethod(ksp,"KSPQCGGetQuadratic_C",(KSP,PetscReal*),(ksp,quadratic));
 89:   return(0);
 90: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

266:   /* Unscale x */
267:   VecCopy(X,WA2);
268:   PCApplySymmetricRight(pc,WA2,X);

270:   KSP_MatMult(ksp,Amat,X,WA);
271:   VecDotRealPart(B,X,&btx);
272:   VecDotRealPart(X,WA,&xtax);

274:   pcgP->quadratic = btx + .5*xtax;
275:   return(0);
276: }

278: PetscErrorCode KSPSetUp_QCG(KSP ksp)
279: {

283:   /* Get work vectors from user code */
284:   KSPSetWorkVecs(ksp,7);
285:   return(0);
286: }

288: PetscErrorCode KSPDestroy_QCG(KSP ksp)
289: {

293:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",NULL);
294:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",NULL);
295:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",NULL);
296:   KSPDestroyDefault(ksp);
297:   return(0);
298: }

300: static PetscErrorCode  KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
301: {
302:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

305:   cgP->delta = delta;
306:   return(0);
307: }

309: static PetscErrorCode  KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
310: {
311:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

314:   *ltsnrm = cgP->ltsnrm;
315:   return(0);
316: }

318: static PetscErrorCode  KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
319: {
320:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

323:   *quadratic = cgP->quadratic;
324:   return(0);
325: }

327: PetscErrorCode KSPSetFromOptions_QCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
328: {
330:   PetscReal      delta;
331:   KSP_QCG        *cgP = (KSP_QCG*)ksp->data;
332:   PetscBool      flg;

335:   PetscOptionsHead(PetscOptionsObject,"KSP QCG Options");
336:   PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
337:   if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
338:   PetscOptionsTail();
339:   return(0);
340: }

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

346:    Options Database Keys:
347: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

349:    Notes:
350:     This is rarely used directly

352:    Level: developer

354:   Notes:
355:     Use preconditioned conjugate gradient to compute
356:       an approximate minimizer of the quadratic function

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

360:    subject to the Euclidean norm trust region constraint

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

364:    where

366:      delta is the trust region radius,
367:      g is the gradient vector, and
368:      H is Hessian matrix,
369:      D is a scaling matrix.

371:    KSPConvergedReason may be
372: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
373: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
374: $  other KSP converged/diverged reasons

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

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

387: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
388:            KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
389: M*/

391: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
392: {
394:   KSP_QCG        *cgP;

397:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,3);
398:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_SYMMETRIC,1);
399:   PetscNewLog(ksp,&cgP);

401:   ksp->data                = (void*)cgP;
402:   ksp->ops->setup          = KSPSetUp_QCG;
403:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
404:   ksp->ops->solve          = KSPSolve_QCG;
405:   ksp->ops->destroy        = KSPDestroy_QCG;
406:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
407:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
408:   ksp->ops->view           = NULL;

410:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",KSPQCGGetQuadratic_QCG);
411:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",KSPQCGGetTrialStepNorm_QCG);
412:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",KSPQCGSetTrustRegionRadius_QCG);
413:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
414:   return(0);
415: }

417: /* ---------------------------------------------------------- */
418: /*
419:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
420:          ||s + step*p|| - delta = 0
421:    such that step1 >= 0 >= step2.
422:    where
423:       delta:
424:         On entry delta must contain scalar delta.
425:         On exit delta is unchanged.
426:       step1:
427:         On entry step1 need not be specified.
428:         On exit step1 contains the non-negative root.
429:       step2:
430:         On entry step2 need not be specified.
431:         On exit step2 contains the non-positive root.
432:    C code is translated from the Fortran version of the MINPACK-2 Project,
433:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
434: */
435: static PetscErrorCode KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal delta,PetscReal *step1,PetscReal *step2)
436: {
437:   PetscReal      dsq,ptp,pts,rad,sts;

441:   VecDotRealPart(p,s,&pts);
442:   VecDotRealPart(p,p,&ptp);
443:   VecDotRealPart(s,s,&sts);
444:   dsq  = delta*delta;
445:   rad  = PetscSqrtReal((pts*pts) - ptp*(sts - dsq));
446:   if (pts > 0.0) {
447:     *step2 = -(pts + rad)/ptp;
448:     *step1 = (sts - dsq)/(ptp * *step2);
449:   } else {
450:     *step1 = -(pts - rad)/ptp;
451:     *step2 = (sts - dsq)/(ptp * *step1);
452:   }
453:   return(0);
454: }