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 for `KSPQCG`
9: Logically Collective
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 in `KSPQCG`. The WCG step may be
31: constrained, so this is not necessarily the length of the ultimate step taken in `KSPQCG`.
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 Euclidean 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
68: Input Parameter:
69: . ksp - the iterative context
71: Output Parameter:
72: . quadratic - the quadratic function evaluated at the new iterate
74: Level: advanced
76: .seealso: [](chapter_ksp), `KSPQCG`
77: @*/
78: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
79: {
81: PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
82: return 0;
83: }
85: PetscErrorCode KSPSolve_QCG(KSP ksp)
86: {
87: /*
88: Correpondence with documentation above:
89: B = g = gradient,
90: X = s = step
91: Note: This is not coded correctly for complex arithmetic!
92: */
94: KSP_QCG *pcgP = (KSP_QCG *)ksp->data;
95: Mat Amat, Pmat;
96: Vec W, WA, WA2, R, P, ASP, BS, X, B;
97: PetscScalar scal, beta, rntrn, step;
98: PetscReal q1, q2, xnorm, step1, step2, rnrm = 0.0, btx, xtax;
99: PetscReal ptasp, rtr, wtasp, bstp;
100: PetscReal dzero = 0.0, bsnrm = 0.0;
101: PetscInt i, maxit;
102: PC pc = ksp->pc;
103: PetscBool diagonalscale;
105: PCGetDiagonalScale(ksp->pc, &diagonalscale);
109: ksp->its = 0;
110: maxit = ksp->max_it;
111: WA = ksp->work[0];
112: R = ksp->work[1];
113: P = ksp->work[2];
114: ASP = ksp->work[3];
115: BS = ksp->work[4];
116: W = ksp->work[5];
117: WA2 = ksp->work[6];
118: X = ksp->vec_sol;
119: B = ksp->vec_rhs;
123: /* Initialize variables */
124: VecSet(W, 0.0); /* W = 0 */
125: VecSet(X, 0.0); /* X = 0 */
126: PCGetOperators(pc, &Amat, &Pmat);
128: /* Compute: BS = D^{-1} B */
129: PCApplySymmetricLeft(pc, B, BS);
131: if (ksp->normtype != KSP_NORM_NONE) {
132: VecNorm(BS, NORM_2, &bsnrm);
133: KSPCheckNorm(ksp, bsnrm);
134: }
135: PetscObjectSAWsTakeAccess((PetscObject)ksp);
136: ksp->its = 0;
137: ksp->rnorm = bsnrm;
138: PetscObjectSAWsGrantAccess((PetscObject)ksp);
139: KSPLogResidualHistory(ksp, bsnrm);
140: KSPMonitor(ksp, 0, bsnrm);
141: (*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP);
142: if (ksp->reason) return 0;
144: /* Compute the initial scaled direction and scaled residual */
145: VecCopy(BS, R);
146: VecScale(R, -1.0);
147: VecCopy(R, P);
148: VecDotRealPart(R, R, &rtr);
150: for (i = 0; i <= maxit; i++) {
151: PetscObjectSAWsTakeAccess((PetscObject)ksp);
152: ksp->its++;
153: PetscObjectSAWsGrantAccess((PetscObject)ksp);
155: /* Compute: asp = D^{-T}*A*D^{-1}*p */
156: PCApplySymmetricRight(pc, P, WA);
157: KSP_MatMult(ksp, Amat, WA, WA2);
158: PCApplySymmetricLeft(pc, WA2, ASP);
160: /* Check for negative curvature */
161: VecDotRealPart(P, ASP, &ptasp);
162: if (ptasp <= dzero) {
163: /* Scaled negative curvature direction: Compute a step so that
164: ||w + step*p|| = delta and QS(w + step*p) is least */
166: if (!i) {
167: VecCopy(P, X);
168: VecNorm(X, NORM_2, &xnorm);
169: KSPCheckNorm(ksp, xnorm);
170: scal = pcgP->delta / xnorm;
171: VecScale(X, scal);
172: } else {
173: /* Compute roots of quadratic */
174: KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2);
175: VecDotRealPart(W, ASP, &wtasp);
176: VecDotRealPart(BS, P, &bstp);
177: VecCopy(W, X);
178: q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
179: q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
180: if (q1 <= q2) {
181: VecAXPY(X, step1, P);
182: } else {
183: VecAXPY(X, step2, P);
184: }
185: }
186: pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
187: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE; /* negative curvature */
188: if (!i) {
189: PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta);
190: } else {
191: PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta);
192: }
194: } else {
195: /* Compute step along p */
196: step = rtr / ptasp;
197: VecCopy(W, X); /* x = w */
198: VecAXPY(X, step, P); /* x <- step*p + x */
199: VecNorm(X, NORM_2, &pcgP->ltsnrm);
200: KSPCheckNorm(ksp, pcgP->ltsnrm);
202: if (pcgP->ltsnrm > pcgP->delta) {
203: /* Since the trial iterate is outside the trust region,
204: evaluate a constrained step along p so that
205: ||w + step*p|| = delta
206: The positive step is always better in this case. */
207: if (!i) {
208: scal = pcgP->delta / pcgP->ltsnrm;
209: VecScale(X, scal);
210: } else {
211: /* Compute roots of quadratic */
212: KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2);
213: VecCopy(W, X);
214: VecAXPY(X, step1, P); /* x <- step1*p + x */
215: }
216: pcgP->ltsnrm = pcgP->delta;
217: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED; /* convergence along constrained step */
218: if (!i) {
219: PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta);
220: } else {
221: PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta);
222: }
224: } else {
225: /* Evaluate the current step */
226: VecCopy(X, W); /* update interior iterate */
227: VecAXPY(R, -step, ASP); /* r <- -step*asp + r */
228: if (ksp->normtype != KSP_NORM_NONE) {
229: VecNorm(R, NORM_2, &rnrm);
230: KSPCheckNorm(ksp, rnrm);
231: }
232: PetscObjectSAWsTakeAccess((PetscObject)ksp);
233: ksp->rnorm = rnrm;
234: PetscObjectSAWsGrantAccess((PetscObject)ksp);
235: KSPLogResidualHistory(ksp, rnrm);
236: KSPMonitor(ksp, i + 1, rnrm);
237: (*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP);
238: if (ksp->reason) { /* convergence for */
239: PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta);
240: }
241: }
242: }
243: if (ksp->reason) break; /* Convergence has been attained */
244: else { /* Compute a new AS-orthogonal direction */ VecDot(R, R, &rntrn);
245: beta = rntrn / rtr;
246: VecAYPX(P, beta, R); /* p <- r + beta*p */
247: rtr = PetscRealPart(rntrn);
248: }
249: }
250: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
252: /* Unscale x */
253: VecCopy(X, WA2);
254: PCApplySymmetricRight(pc, WA2, X);
256: KSP_MatMult(ksp, Amat, X, WA);
257: VecDotRealPart(B, X, &btx);
258: VecDotRealPart(X, WA, &xtax);
260: pcgP->quadratic = btx + .5 * xtax;
261: return 0;
262: }
264: PetscErrorCode KSPSetUp_QCG(KSP ksp)
265: {
266: /* Get work vectors from user code */
267: KSPSetWorkVecs(ksp, 7);
268: return 0;
269: }
271: PetscErrorCode KSPDestroy_QCG(KSP ksp)
272: {
273: PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL);
274: PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL);
275: PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL);
276: KSPDestroyDefault(ksp);
277: return 0;
278: }
280: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
281: {
282: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
284: cgP->delta = delta;
285: return 0;
286: }
288: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
289: {
290: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
292: *ltsnrm = cgP->ltsnrm;
293: return 0;
294: }
296: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
297: {
298: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
300: *quadratic = cgP->quadratic;
301: return 0;
302: }
304: PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
305: {
306: PetscReal delta;
307: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
308: PetscBool flg;
310: PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
311: PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg);
312: if (flg) KSPQCGSetTrustRegionRadius(ksp, delta);
313: PetscOptionsHeadEnd();
314: return 0;
315: }
317: /*MC
318: KSPQCG - Code to run conjugate gradient method subject to a constraint on the solution norm.
320: Options Database Key:
321: . -ksp_qcg_trustregionradius <r> - Trust Region Radius
323: Level: developer
325: Notes:
326: This is rarely used directly, ir is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
328: Uses preconditioned conjugate gradient to compute
329: an approximate minimizer of the quadratic function
331: q(s) = g^T * s + .5 * s^T * H * s
333: subject to the Euclidean norm trust region constraint
335: || D * s || <= delta,
337: where
339: delta is the trust region radius,
340: g is the gradient vector, and
341: H is Hessian matrix,
342: D is a scaling matrix.
344: `KSPConvergedReason` may be
345: .vb
346: KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
347: KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
348: .ve
349: and other `KSP` converged/diverged reasons
351: Notes:
352: Currently we allow symmetric preconditioning with the following scaling matrices:
353: .vb
354: `PCNONE`: D = Identity matrix
355: `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
356: `PCICC`: D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
357: .ve
359: References:
360: . * - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
361: SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).
363: .seealso: [](chapter_ksp), 'KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`
364: `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
365: M*/
367: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
368: {
369: KSP_QCG *cgP;
371: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3);
372: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1);
373: PetscNew(&cgP);
375: ksp->data = (void *)cgP;
376: ksp->ops->setup = KSPSetUp_QCG;
377: ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
378: ksp->ops->solve = KSPSolve_QCG;
379: ksp->ops->destroy = KSPDestroy_QCG;
380: ksp->ops->buildsolution = KSPBuildSolutionDefault;
381: ksp->ops->buildresidual = KSPBuildResidualDefault;
382: ksp->ops->view = NULL;
384: PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG);
385: PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG);
386: PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG);
387: cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
388: return 0;
389: }
391: /* ---------------------------------------------------------- */
392: /*
393: KSPQCGQuadraticRoots - Computes the roots of the quadratic,
394: ||s + step*p|| - delta = 0
395: such that step1 >= 0 >= step2.
396: where
397: delta:
398: On entry delta must contain scalar delta.
399: On exit delta is unchanged.
400: step1:
401: On entry step1 need not be specified.
402: On exit step1 contains the non-negative root.
403: step2:
404: On entry step2 need not be specified.
405: On exit step2 contains the non-positive root.
406: C code is translated from the Fortran version of the MINPACK-2 Project,
407: Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
408: */
409: static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
410: {
411: PetscReal dsq, ptp, pts, rad, sts;
413: VecDotRealPart(p, s, &pts);
414: VecDotRealPart(p, p, &ptp);
415: VecDotRealPart(s, s, &sts);
416: dsq = delta * delta;
417: rad = PetscSqrtReal((pts * pts) - ptp * (sts - dsq));
418: if (pts > 0.0) {
419: *step2 = -(pts + rad) / ptp;
420: *step1 = (sts - dsq) / (ptp * *step2);
421: } else {
422: *step1 = -(pts - rad) / ptp;
423: *step2 = (sts - dsq) / (ptp * *step1);
424: }
425: return 0;
426: }