Actual source code: qcg.c
petsc-3.5.4 2015-05-23
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: }