Actual source code: qcg.c
petsc-3.3-p7 2013-05-11
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(((PetscObject)ksp)->comm,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: MatStructure pflag;
114: Mat Amat,Pmat;
115: Vec W,WA,WA2,R,P,ASP,BS,X,B;
116: PetscScalar scal,btx,xtax,beta,rntrn,step;
117: PetscReal ptasp,q1,q2,wtasp,bstp,rtr,xnorm,step1,step2,rnrm,p5 = 0.5;
118: PetscReal dzero = 0.0,bsnrm;
120: PetscInt i,maxit;
121: PC pc = ksp->pc;
122: PCSide side;
123: #if defined(PETSC_USE_COMPLEX)
124: PetscScalar cstep1,cstep2,cbstp,crtr,cwtasp,cptasp;
125: #endif
126: PetscBool diagonalscale;
129: PCGetDiagonalScale(ksp->pc,&diagonalscale);
130: if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
131: if (ksp->transpose_solve) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Currently does not support transpose solve");
133: ksp->its = 0;
134: maxit = ksp->max_it;
135: WA = ksp->work[0];
136: R = ksp->work[1];
137: P = ksp->work[2];
138: ASP = ksp->work[3];
139: BS = ksp->work[4];
140: W = ksp->work[5];
141: WA2 = ksp->work[6];
142: X = ksp->vec_sol;
143: B = ksp->vec_rhs;
145: if (pcgP->delta <= dzero) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Input error: delta <= 0");
146: KSPGetPCSide(ksp,&side);
147: if (side != PC_SYMMETRIC) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Requires symmetric preconditioner!");
149: /* Initialize variables */
150: VecSet(W,0.0); /* W = 0 */
151: VecSet(X,0.0); /* X = 0 */
152: PCGetOperators(pc,&Amat,&Pmat,&pflag);
154: /* Compute: BS = D^{-1} B */
155: PCApplySymmetricLeft(pc,B,BS);
157: VecNorm(BS,NORM_2,&bsnrm);
158: PetscObjectTakeAccess(ksp);
159: ksp->its = 0;
160: ksp->rnorm = bsnrm;
161: PetscObjectGrantAccess(ksp);
162: KSPLogResidualHistory(ksp,bsnrm);
163: KSPMonitor(ksp,0,bsnrm);
164: (*ksp->converged)(ksp,0,bsnrm,&ksp->reason,ksp->cnvP);
165: if (ksp->reason) return(0);
167: /* Compute the initial scaled direction and scaled residual */
168: VecCopy(BS,R);
169: VecScale(R,-1.0);
170: VecCopy(R,P);
171: #if defined(PETSC_USE_COMPLEX)
172: VecDot(R,R,&crtr); rtr = PetscRealPart(crtr);
173: #else
174: VecDot(R,R,&rtr);
175: #endif
177: for (i=0; i<=maxit; i++) {
178: PetscObjectTakeAccess(ksp);
179: ksp->its++;
180: PetscObjectGrantAccess(ksp);
182: /* Compute: asp = D^{-T}*A*D^{-1}*p */
183: PCApplySymmetricRight(pc,P,WA);
184: MatMult(Amat,WA,WA2);
185: PCApplySymmetricLeft(pc,WA2,ASP);
187: /* Check for negative curvature */
188: #if defined(PETSC_USE_COMPLEX)
189: VecDot(P,ASP,&cptasp);
190: ptasp = PetscRealPart(cptasp);
191: #else
192: VecDot(P,ASP,&ptasp); /* ptasp = p^T asp */
193: #endif
194: if (ptasp <= dzero) {
196: /* Scaled negative curvature direction: Compute a step so that
197: ||w + step*p|| = delta and QS(w + step*p) is least */
199: if (!i) {
200: VecCopy(P,X);
201: VecNorm(X,NORM_2,&xnorm);
202: scal = pcgP->delta / xnorm;
203: VecScale(X,scal);
204: } else {
205: /* Compute roots of quadratic */
206: KSPQCGQuadraticRoots(W,P,&pcgP->delta,&step1,&step2);
207: #if defined(PETSC_USE_COMPLEX)
208: VecDot(W,ASP,&cwtasp); wtasp = PetscRealPart(cwtasp);
209: VecDot(BS,P,&cbstp); bstp = PetscRealPart(cbstp);
210: #else
211: VecDot(W,ASP,&wtasp);
212: VecDot(BS,P,&bstp);
213: #endif
214: VecCopy(W,X);
215: q1 = step1*(bstp + wtasp + p5*step1*ptasp);
216: q2 = step2*(bstp + wtasp + p5*step2*ptasp);
217: #if defined(PETSC_USE_COMPLEX)
218: if (q1 <= q2) {
219: cstep1 = step1; VecAXPY(X,cstep1,P);
220: } else {
221: cstep2 = step2; VecAXPY(X,cstep2,P);
222: }
223: #else
224: if (q1 <= q2) {VecAXPY(X,step1,P);}
225: else {VecAXPY(X,step2,P);}
226: #endif
227: }
228: pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
229: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE; /* negative curvature */
230: if (!i) {
231: PetscInfo1(ksp,"negative curvature: delta=%G\n",pcgP->delta);
232: } else {
233: PetscInfo3(ksp,"negative curvature: step1=%G, step2=%G, delta=%G\n",step1,step2,pcgP->delta);
234: }
235:
236: } else {
237:
238: /* Compute step along p */
240: step = rtr/ptasp;
241: VecCopy(W,X); /* x = w */
242: VecAXPY(X,step,P); /* x <- step*p + x */
243: VecNorm(X,NORM_2,&pcgP->ltsnrm);
245: if (pcgP->ltsnrm > pcgP->delta) {
247: /* Since the trial iterate is outside the trust region,
248: evaluate a constrained step along p so that
249: ||w + step*p|| = delta
250: The positive step is always better in this case. */
252: if (!i) {
253: scal = pcgP->delta / pcgP->ltsnrm;
254: VecScale(X,scal);
255: } else {
256: /* Compute roots of quadratic */
257: KSPQCGQuadraticRoots(W,P,&pcgP->delta,&step1,&step2);
258: VecCopy(W,X);
259: #if defined(PETSC_USE_COMPLEX)
260: cstep1 = step1; VecAXPY(X,cstep1,P);
261: #else
262: VecAXPY(X,step1,P); /* x <- step1*p + x */
263: #endif
264: }
265: pcgP->ltsnrm = pcgP->delta;
266: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED; /* convergence along constrained step */
267: if (!i) {
268: PetscInfo1(ksp,"constrained step: delta=%G\n",pcgP->delta);
269: } else {
270: PetscInfo3(ksp,"constrained step: step1=%G, step2=%G, delta=%G\n",step1,step2,pcgP->delta);
271: }
273: } else {
275: /* Evaluate the current step */
277: VecCopy(X,W); /* update interior iterate */
278: VecAXPY(R,-step,ASP); /* r <- -step*asp + r */
279: VecNorm(R,NORM_2,&rnrm);
281: PetscObjectTakeAccess(ksp);
282: ksp->rnorm = rnrm;
283: PetscObjectGrantAccess(ksp);
284: KSPLogResidualHistory(ksp,rnrm);
285: KSPMonitor(ksp,i+1,rnrm);
286: (*ksp->converged)(ksp,i+1,rnrm,&ksp->reason,ksp->cnvP);
287: if (ksp->reason) { /* convergence for */
288: #if defined(PETSC_USE_COMPLEX)
289: PetscInfo3(ksp,"truncated step: step=%G, rnrm=%G, delta=%G\n",PetscRealPart(step),rnrm,pcgP->delta);
290: #else
291: PetscInfo3(ksp,"truncated step: step=%G, rnrm=%G, delta=%G\n",step,rnrm,pcgP->delta);
292: #endif
293: }
294: }
295: }
296: if (ksp->reason) break; /* Convergence has been attained */
297: else { /* Compute a new AS-orthogonal direction */
298: VecDot(R,R,&rntrn);
299: beta = rntrn/rtr;
300: VecAYPX(P,beta,R); /* p <- r + beta*p */
301: #if defined(PETSC_USE_COMPLEX)
302: rtr = PetscRealPart(rntrn);
303: #else
304: rtr = rntrn;
305: #endif
306: }
307: }
308: if (!ksp->reason) {
309: ksp->reason = KSP_DIVERGED_ITS;
310: }
312: /* Unscale x */
313: VecCopy(X,WA2);
314: PCApplySymmetricRight(pc,WA2,X);
316: MatMult(Amat,X,WA);
317: VecDot(B,X,&btx);
318: VecDot(X,WA,&xtax);
319: #if defined(PETSC_USE_COMPLEX)
320: pcgP->quadratic = PetscRealPart(btx) + p5* PetscRealPart(xtax);
321: #else
322: pcgP->quadratic = btx + p5*xtax; /* Compute q(x) */
323: #endif
324: return(0);
325: }
329: PetscErrorCode KSPSetUp_QCG(KSP ksp)
330: {
334: /* Get work vectors from user code */
335: KSPDefaultGetWork(ksp,7);
336: return(0);
337: }
341: PetscErrorCode KSPDestroy_QCG(KSP ksp)
342: {
346: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGGetQuadratic_C","",PETSC_NULL);
347: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C","",PETSC_NULL);
348: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C","",PETSC_NULL);
349: KSPDefaultDestroy(ksp);
350: return(0);
351: }
353: EXTERN_C_BEGIN
356: PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
357: {
358: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
361: cgP->delta = delta;
362: return(0);
363: }
364: EXTERN_C_END
366: EXTERN_C_BEGIN
369: PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
370: {
371: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
374: *ltsnrm = cgP->ltsnrm;
375: return(0);
376: }
377: EXTERN_C_END
379: EXTERN_C_BEGIN
382: PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
383: {
384: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
387: *quadratic = cgP->quadratic;
388: return(0);
389: }
390: EXTERN_C_END
394: PetscErrorCode KSPSetFromOptions_QCG(KSP ksp)
395: {
397: PetscReal delta;
398: KSP_QCG *cgP = (KSP_QCG*)ksp->data;
399: PetscBool flg;
402: PetscOptionsHead("KSP QCG Options");
403: PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
404: if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
405: PetscOptionsTail();
406: return(0);
407: }
409: /*MC
410: KSPQCG - Code to run conjugate gradient method subject to a constraint
411: on the solution norm. This is used in Trust Region methods for nonlinear equations, SNESTR
413: Options Database Keys:
414: . -ksp_qcg_trustregionradius <r> - Trust Region Radius
416: Notes: This is rarely used directly
418: Level: developer
420: Notes: Use preconditioned conjugate gradient to compute
421: an approximate minimizer of the quadratic function
423: q(s) = g^T * s + .5 * s^T * H * s
425: subject to the Euclidean norm trust region constraint
427: || D * s || <= delta,
429: where
431: delta is the trust region radius,
432: g is the gradient vector, and
433: H is Hessian matrix,
434: D is a scaling matrix.
436: KSPConvergedReason may be
437: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
438: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
439: $ other KSP converged/diverged reasons
442: Notes:
443: Currently we allow symmetric preconditioning with the following scaling matrices:
444: PCNONE: D = Identity matrix
445: PCJACOBI: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
446: PCICC: D = L^T, implemented with forward and backward solves.
447: Here L is an incomplete Cholesky factor of H.
449: References:
450: The Conjugate Gradient Method and Trust Regions in Large Scale Optimization, Trond Steihaug
451: SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983), pp. 626-637
453: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
454: KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
455: M*/
457: EXTERN_C_BEGIN
460: PetscErrorCode KSPCreate_QCG(KSP ksp)
461: {
463: KSP_QCG *cgP;
466: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,2);
467: PetscNewLog(ksp,KSP_QCG,&cgP);
468: ksp->data = (void*)cgP;
469: ksp->ops->setup = KSPSetUp_QCG;
470: ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
471: ksp->ops->solve = KSPSolve_QCG;
472: ksp->ops->destroy = KSPDestroy_QCG;
473: ksp->ops->buildsolution = KSPDefaultBuildSolution;
474: ksp->ops->buildresidual = KSPDefaultBuildResidual;
475: ksp->ops->setfromoptions = 0;
476: ksp->ops->view = 0;
478: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGGetQuadratic_C",
479: "KSPQCGGetQuadratic_QCG",
480: KSPQCGGetQuadratic_QCG);
481: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",
482: "KSPQCGGetTrialStepNorm_QCG",
483: KSPQCGGetTrialStepNorm_QCG);
484: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",
485: "KSPQCGSetTrustRegionRadius_QCG",
486: KSPQCGSetTrustRegionRadius_QCG);
487: cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
488: return(0);
489: }
490: EXTERN_C_END
492: /* ---------------------------------------------------------- */
495: /*
496: KSPQCGQuadraticRoots - Computes the roots of the quadratic,
497: ||s + step*p|| - delta = 0
498: such that step1 >= 0 >= step2.
499: where
500: delta:
501: On entry delta must contain scalar delta.
502: On exit delta is unchanged.
503: step1:
504: On entry step1 need not be specified.
505: On exit step1 contains the non-negative root.
506: step2:
507: On entry step2 need not be specified.
508: On exit step2 contains the non-positive root.
509: C code is translated from the Fortran version of the MINPACK-2 Project,
510: Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
511: */
512: static PetscErrorCode KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal *delta,PetscReal *step1,PetscReal *step2)
513: {
514: PetscReal dsq,ptp,pts,rad,sts;
516: #if defined(PETSC_USE_COMPLEX)
517: PetscScalar cptp,cpts,csts;
518: #endif
521: #if defined(PETSC_USE_COMPLEX)
522: VecDot(p,s,&cpts); pts = PetscRealPart(cpts);
523: VecDot(p,p,&cptp); ptp = PetscRealPart(cptp);
524: VecDot(s,s,&csts); sts = PetscRealPart(csts);
525: #else
526: VecDot(p,s,&pts);
527: VecDot(p,p,&ptp);
528: VecDot(s,s,&sts);
529: #endif
530: dsq = (*delta)*(*delta);
531: rad = PetscSqrtReal((pts*pts) - ptp*(sts - dsq));
532: if (pts > 0.0) {
533: *step2 = -(pts + rad)/ptp;
534: *step1 = (sts - dsq)/(ptp * *step2);
535: } else {
536: *step1 = -(pts - rad)/ptp;
537: *step2 = (sts - dsq)/(ptp * *step1);
538: }
539: return(0);
540: }