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: }