Actual source code: qcg.c
petsc-3.11.4 2019-09-28
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: .keywords: KSP, QCG, set, trust region radius
21: @*/
22: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp,PetscReal delta)
23: {
28: if (delta < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
29: PetscTryMethod(ksp,"KSPQCGSetTrustRegionRadius_C",(KSP,PetscReal),(ksp,delta));
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
45: Level: advanced
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
81: Level: advanced
82: @*/
83: PetscErrorCode KSPQCGGetQuadratic(KSP ksp,PetscReal *quadratic)
84: {
89: PetscUseMethod(ksp,"KSPQCGGetQuadratic_C",(KSP,PetscReal*),(ksp,quadratic));
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: KSPCheckNorm(ksp,bsnrm);
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: VecNorm(R,NORM_2,&rnrm);
242: KSPCheckNorm(ksp,rnrm);
243: PetscObjectSAWsTakeAccess((PetscObject)ksp);
244: ksp->rnorm = rnrm;
245: PetscObjectSAWsGrantAccess((PetscObject)ksp);
246: KSPLogResidualHistory(ksp,rnrm);
247: KSPMonitor(ksp,i+1,rnrm);
248: (*ksp->converged)(ksp,i+1,rnrm,&ksp->reason,ksp->cnvP);
249: if (ksp->reason) { /* convergence for */
250: PetscInfo3(ksp,"truncated step: step=%g, rnrm=%g, delta=%g\n",(double)PetscRealPart(step),(double)rnrm,(double)pcgP->delta);
251: }
252: }
253: }
254: if (ksp->reason) break; /* Convergence has been attained */
255: else { /* Compute a new AS-orthogonal direction */
256: VecDot(R,R,&rntrn);
257: beta = rntrn/rtr;
258: VecAYPX(P,beta,R); /* p <- r + beta*p */
259: rtr = PetscRealPart(rntrn);
260: }
261: }
262: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
264: /* Unscale x */
265: VecCopy(X,WA2);
266: PCApplySymmetricRight(pc,WA2,X);
268: KSP_MatMult(ksp,Amat,X,WA);
269: VecDotRealPart(B,X,&btx);
270: VecDotRealPart(X,WA,&xtax);
272: pcgP->quadratic = btx + .5*xtax;
273: return(0);
274: }
276: PetscErrorCode KSPSetUp_QCG(KSP ksp)
277: {
281: /* Get work vectors from user code */
282: KSPSetWorkVecs(ksp,7);
283: return(0);
284: }
286: PetscErrorCode KSPDestroy_QCG(KSP ksp)
287: {
291: PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",NULL);
292: PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",NULL);
293: PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",NULL);
294: KSPDestroyDefault(ksp);
295: return(0);
296: }
298: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
299: {
300: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
303: cgP->delta = delta;
304: return(0);
305: }
307: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
308: {
309: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
312: *ltsnrm = cgP->ltsnrm;
313: return(0);
314: }
316: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
317: {
318: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
321: *quadratic = cgP->quadratic;
322: return(0);
323: }
325: PetscErrorCode KSPSetFromOptions_QCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
326: {
328: PetscReal delta;
329: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
330: PetscBool flg;
333: PetscOptionsHead(PetscOptionsObject,"KSP QCG Options");
334: PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
335: if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
336: PetscOptionsTail();
337: return(0);
338: }
340: /*MC
341: KSPQCG - Code to run conjugate gradient method subject to a constraint
342: on the solution norm. This is used in Trust Region methods for nonlinear equations, SNESNEWTONTR
344: Options Database Keys:
345: . -ksp_qcg_trustregionradius <r> - Trust Region Radius
347: Notes:
348: This is rarely used directly
350: Level: developer
352: Notes:
353: Use preconditioned conjugate gradient to compute
354: an approximate minimizer of the quadratic function
356: q(s) = g^T * s + .5 * s^T * H * s
358: subject to the Euclidean norm trust region constraint
360: || D * s || <= delta,
362: where
364: delta is the trust region radius,
365: g is the gradient vector, and
366: H is Hessian matrix,
367: D is a scaling matrix.
369: KSPConvergedReason may be
370: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
371: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
372: $ other KSP converged/diverged reasons
375: Notes:
376: Currently we allow symmetric preconditioning with the following scaling matrices:
377: PCNONE: D = Identity matrix
378: PCJACOBI: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
379: PCICC: D = L^T, implemented with forward and backward solves.
380: Here L is an incomplete Cholesky factor of H.
382: References:
383: . 1. - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
384: SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).
386: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
387: KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
388: M*/
390: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
391: {
393: KSP_QCG *cgP;
396: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,3);
397: PetscNewLog(ksp,&cgP);
399: ksp->data = (void*)cgP;
400: ksp->ops->setup = KSPSetUp_QCG;
401: ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
402: ksp->ops->solve = KSPSolve_QCG;
403: ksp->ops->destroy = KSPDestroy_QCG;
404: ksp->ops->buildsolution = KSPBuildSolutionDefault;
405: ksp->ops->buildresidual = KSPBuildResidualDefault;
406: ksp->ops->setfromoptions = 0;
407: ksp->ops->view = 0;
409: PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",KSPQCGGetQuadratic_QCG);
410: PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",KSPQCGGetTrialStepNorm_QCG);
411: PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",KSPQCGSetTrustRegionRadius_QCG);
412: cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
413: return(0);
414: }
416: /* ---------------------------------------------------------- */
417: /*
418: KSPQCGQuadraticRoots - Computes the roots of the quadratic,
419: ||s + step*p|| - delta = 0
420: such that step1 >= 0 >= step2.
421: where
422: delta:
423: On entry delta must contain scalar delta.
424: On exit delta is unchanged.
425: step1:
426: On entry step1 need not be specified.
427: On exit step1 contains the non-negative root.
428: step2:
429: On entry step2 need not be specified.
430: On exit step2 contains the non-positive root.
431: C code is translated from the Fortran version of the MINPACK-2 Project,
432: Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
433: */
434: static PetscErrorCode KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal delta,PetscReal *step1,PetscReal *step2)
435: {
436: PetscReal dsq,ptp,pts,rad,sts;
440: VecDotRealPart(p,s,&pts);
441: VecDotRealPart(p,p,&ptp);
442: VecDotRealPart(s,s,&sts);
443: dsq = delta*delta;
444: rad = PetscSqrtReal((pts*pts) - ptp*(sts - dsq));
445: if (pts > 0.0) {
446: *step2 = -(pts + rad)/ptp;
447: *step1 = (sts - dsq)/(ptp * *step2);
448: } else {
449: *step1 = -(pts - rad)/ptp;
450: *step2 = (sts - dsq)/(ptp * *step1);
451: }
452: return(0);
453: }