Actual source code: qcg.c
1: #include <../src/ksp/ksp/impls/qcg/qcgimpl.h>
3: /*
4: KSPQCGQuadraticRoots - Computes the roots of the quadratic,
5: ||s + step*p|| - delta = 0
6: such that step1 >= 0 >= step2.
7: where
8: delta:
9: On entry delta must contain scalar delta.
10: On exit delta is unchanged.
11: step1:
12: On entry step1 need not be specified.
13: On exit step1 contains the non-negative root.
14: step2:
15: On entry step2 need not be specified.
16: On exit step2 contains the non-positive root.
17: C code is translated from the Fortran version of the MINPACK-2 Project,
18: Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
19: */
20: static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
21: {
22: PetscReal dsq, ptp, pts, rad, sts;
24: PetscFunctionBegin;
25: PetscCall(VecDotRealPart(p, s, &pts));
26: PetscCall(VecDotRealPart(p, p, &ptp));
27: PetscCall(VecDotRealPart(s, s, &sts));
28: dsq = delta * delta;
29: rad = PetscSqrtReal((pts * pts) - ptp * (sts - dsq));
30: if (pts > 0.0) {
31: *step2 = -(pts + rad) / ptp;
32: *step1 = (sts - dsq) / (ptp * *step2);
33: } else {
34: *step1 = -(pts - rad) / ptp;
35: *step2 = (sts - dsq) / (ptp * *step1);
36: }
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: /*@
41: KSPQCGSetTrustRegionRadius - Sets the radius of the trust region for `KSPQCG`
43: Logically Collective
45: Input Parameters:
46: + ksp - the iterative context
47: - delta - the trust region radius (Infinity is the default)
49: Options Database Key:
50: . -ksp_qcg_trustregionradius delta - trust region radius
52: Level: advanced
54: Developer Note:
55: `KSPMINRESSetRadius()`, for example, does not have TrustRegion in the name
57: .seealso: [](ch_ksp), `KSPQCG`, `KSPQCGGetTrialStepNorm()`
58: @*/
59: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
60: {
61: PetscFunctionBegin;
63: PetscCheck(delta >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
64: PetscTryMethod(ksp, "KSPQCGSetTrustRegionRadius_C", (KSP, PetscReal), (ksp, delta));
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*@
69: KSPQCGGetTrialStepNorm - Gets the norm of a trial step vector in `KSPQCG`. The WCG step may be
70: constrained, so this is not necessarily the length of the ultimate step taken in `KSPQCG`.
72: Not Collective
74: Input Parameter:
75: . ksp - the iterative context
77: Output Parameter:
78: . tsnorm - the norm
80: Level: advanced
82: .seealso: [](ch_ksp), `KSPQCG`, `KSPQCGSetTrustRegionRadius()`
83: @*/
84: PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
85: {
86: PetscFunctionBegin;
88: PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: /*@
93: KSPQCGGetQuadratic - Gets the value of the quadratic function, evaluated at the new iterate
95: Collective
97: Input Parameter:
98: . ksp - the iterative context
100: Output Parameter:
101: . quadratic - the quadratic function evaluated at the new iterate
103: Level: advanced
105: Note:
106: The quadratic function is
108: $$
109: q(s) = g^T * s + 0.5 * s^T * H * s
110: $$
112: which satisfies the Euclidean Norm trust region constraint
114: $$
115: || D * s || \le delta,
116: $$
118: where
119: .vb
120: delta is the trust region radius,
121: g is the gradient vector, and
122: H is Hessian matrix,
123: D is a scaling matrix.
124: .ve
126: .seealso: [](ch_ksp), `KSPQCG`
127: @*/
128: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
129: {
130: PetscFunctionBegin;
132: PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: static PetscErrorCode KSPSolve_QCG(KSP ksp)
137: {
138: /*
139: Correspondence with documentation above:
140: B = g = gradient,
141: X = s = step
142: Note: This is not coded correctly for complex arithmetic!
143: */
145: KSP_QCG *pcgP = (KSP_QCG *)ksp->data;
146: Mat Amat, Pmat;
147: Vec W, WA, WA2, R, P, ASP, BS, X, B;
148: PetscScalar scal, beta, rntrn, step;
149: PetscReal q1, q2, xnorm, step1, step2, rnrm = 0.0, btx, xtax;
150: PetscReal ptasp, rtr, wtasp, bstp;
151: PetscReal dzero = 0.0, bsnrm = 0.0;
152: PetscInt maxit;
153: PC pc = ksp->pc;
155: PetscFunctionBegin;
156: PetscCheck(!ksp->transpose_solve, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Currently does not support transpose solve");
158: ksp->its = 0;
159: maxit = ksp->max_it;
160: WA = ksp->work[0];
161: R = ksp->work[1];
162: P = ksp->work[2];
163: ASP = ksp->work[3];
164: BS = ksp->work[4];
165: W = ksp->work[5];
166: WA2 = ksp->work[6];
167: X = ksp->vec_sol;
168: B = ksp->vec_rhs;
170: PetscCheck(pcgP->delta > dzero, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: delta <= 0");
172: /* Initialize variables */
173: PetscCall(VecSet(W, 0.0)); /* W = 0 */
174: PetscCall(VecSet(X, 0.0)); /* X = 0 */
175: PetscCall(PCGetOperators(pc, &Amat, &Pmat));
177: /* Compute: BS = D^{-1} B */
178: PetscCall(PCApplySymmetricLeft(pc, B, BS));
180: if (ksp->normtype != KSP_NORM_NONE) {
181: PetscCall(VecNorm(BS, NORM_2, &bsnrm));
182: KSPCheckNorm(ksp, bsnrm);
183: }
184: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
185: ksp->its = 0;
186: ksp->rnorm = bsnrm;
187: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
188: PetscCall(KSPLogResidualHistory(ksp, bsnrm));
189: PetscCall(KSPMonitor(ksp, 0, bsnrm));
190: PetscCall((*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP));
191: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
193: /* Compute the initial scaled direction and scaled residual */
194: PetscCall(VecCopy(BS, R));
195: PetscCall(VecScale(R, -1.0));
196: PetscCall(VecCopy(R, P));
197: PetscCall(VecDotRealPart(R, R, &rtr));
199: for (PetscInt i = 0; i <= maxit; i++) {
200: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
201: ksp->its++;
202: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
204: /* Compute: asp = D^{-T}*A*D^{-1}*p */
205: PetscCall(PCApplySymmetricRight(pc, P, WA));
206: PetscCall(KSP_MatMult(ksp, Amat, WA, WA2));
207: PetscCall(PCApplySymmetricLeft(pc, WA2, ASP));
209: /* Check for negative curvature */
210: PetscCall(VecDotRealPart(P, ASP, &ptasp));
211: if (ptasp <= dzero) {
212: /* Scaled negative curvature direction: Compute a step so that
213: ||w + step*p|| = delta and QS(w + step*p) is least */
215: if (!i) {
216: PetscCall(VecCopy(P, X));
217: PetscCall(VecNorm(X, NORM_2, &xnorm));
218: KSPCheckNorm(ksp, xnorm);
219: scal = pcgP->delta / xnorm;
220: PetscCall(VecScale(X, scal));
221: } else {
222: /* Compute roots of quadratic */
223: PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
224: PetscCall(VecDotRealPart(W, ASP, &wtasp));
225: PetscCall(VecDotRealPart(BS, P, &bstp));
226: PetscCall(VecCopy(W, X));
227: q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
228: q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
229: if (q1 <= q2) PetscCall(VecAXPY(X, step1, P));
230: else PetscCall(VecAXPY(X, step2, P));
231: }
232: pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
233: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
234: if (!i) {
235: PetscCall(PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta));
236: } else {
237: PetscCall(PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
238: }
240: } else {
241: /* Compute step along p */
242: step = rtr / ptasp;
243: PetscCall(VecCopy(W, X)); /* x = w */
244: PetscCall(VecAXPY(X, step, P)); /* x <- step*p + x */
245: PetscCall(VecNorm(X, NORM_2, &pcgP->ltsnrm));
246: KSPCheckNorm(ksp, pcgP->ltsnrm);
248: if (pcgP->ltsnrm > pcgP->delta) {
249: /* Since the trial iterate is outside the trust region,
250: evaluate a constrained step along p so that
251: ||w + step*p|| = delta
252: The positive step is always better in this case. */
253: if (!i) {
254: scal = pcgP->delta / pcgP->ltsnrm;
255: PetscCall(VecScale(X, scal));
256: } else {
257: /* Compute roots of quadratic */
258: PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
259: PetscCall(VecCopy(W, X));
260: PetscCall(VecAXPY(X, step1, P)); /* x <- step1*p + x */
261: }
262: pcgP->ltsnrm = pcgP->delta;
263: ksp->reason = KSP_CONVERGED_STEP_LENGTH; /* convergence along constrained step */
264: if (!i) {
265: PetscCall(PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta));
266: } else {
267: PetscCall(PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
268: }
270: } else {
271: /* Evaluate the current step */
272: PetscCall(VecCopy(X, W)); /* update interior iterate */
273: PetscCall(VecAXPY(R, -step, ASP)); /* r <- -step*asp + r */
274: if (ksp->normtype != KSP_NORM_NONE) {
275: PetscCall(VecNorm(R, NORM_2, &rnrm));
276: KSPCheckNorm(ksp, rnrm);
277: }
278: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
279: ksp->rnorm = rnrm;
280: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
281: PetscCall(KSPLogResidualHistory(ksp, rnrm));
282: PetscCall(KSPMonitor(ksp, i + 1, rnrm));
283: PetscCall((*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP));
284: if (ksp->reason) { /* convergence for */
285: PetscCall(PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta));
286: }
287: }
288: }
289: if (ksp->reason) break; /* Convergence has been attained */
290: else { /* Compute a new AS-orthogonal direction */ PetscCall(VecDot(R, R, &rntrn));
291: beta = rntrn / rtr;
292: PetscCall(VecAYPX(P, beta, R)); /* p <- r + beta*p */
293: rtr = PetscRealPart(rntrn);
294: }
295: }
296: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
298: /* Unscale x */
299: PetscCall(VecCopy(X, WA2));
300: PetscCall(PCApplySymmetricRight(pc, WA2, X));
302: PetscCall(KSP_MatMult(ksp, Amat, X, WA));
303: PetscCall(VecDotRealPart(B, X, &btx));
304: PetscCall(VecDotRealPart(X, WA, &xtax));
306: pcgP->quadratic = btx + .5 * xtax;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode KSPSetUp_QCG(KSP ksp)
311: {
312: PetscFunctionBegin;
313: /* Get work vectors from user code */
314: PetscCall(KSPSetWorkVecs(ksp, 7));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode KSPDestroy_QCG(KSP ksp)
319: {
320: PetscFunctionBegin;
321: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL));
322: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL));
323: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL));
324: PetscCall(KSPDestroyDefault(ksp));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
329: {
330: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
332: PetscFunctionBegin;
333: cgP->delta = delta;
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
338: {
339: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
341: PetscFunctionBegin;
342: *ltsnrm = cgP->ltsnrm;
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
347: {
348: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
350: PetscFunctionBegin;
351: *quadratic = cgP->quadratic;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems PetscOptionsObject)
356: {
357: PetscReal delta;
358: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
359: PetscBool flg;
361: PetscFunctionBegin;
362: PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
363: PetscCall(PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg));
364: if (flg) PetscCall(KSPQCGSetTrustRegionRadius(ksp, delta));
365: PetscOptionsHeadEnd();
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /*MC
370: KSPQCG - Code to run conjugate gradient method subject to a constraint on the solution norm {cite}`steihaug:83`.
372: Options Database Key:
373: . -ksp_qcg_trustregionradius r - Trust Region Radius
375: Level: developer
377: Notes:
378: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
380: Uses preconditioned conjugate gradient to compute
381: an approximate minimizer of the quadratic function $ q(s) = g^T * s + .5 * s^T * H * s $ subject to the Euclidean norm trust region constraint
382: $ || D * s || \le delta$, where
383: .vb
384: delta is the trust region radius,
385: g is the gradient vector, and
386: H is Hessian matrix,
387: D is a scaling matrix.
388: .ve
390: `KSPConvergedReason` may include
391: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
392: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,
394: Note:
395: Allows symmetric preconditioning with the following scaling matrices:
396: .vb
397: `PCNONE`: D = Identity matrix
398: `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
399: `PCICC`: D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
400: .ve
402: .seealso: [](ch_ksp), `KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`,
403: `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
404: M*/
406: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
407: {
408: KSP_QCG *cgP;
410: PetscFunctionBegin;
411: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3));
412: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1));
413: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
414: PetscCall(PetscNew(&cgP));
416: ksp->data = (void *)cgP;
417: ksp->ops->setup = KSPSetUp_QCG;
418: ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
419: ksp->ops->solve = KSPSolve_QCG;
420: ksp->ops->destroy = KSPDestroy_QCG;
421: ksp->ops->buildsolution = KSPBuildSolutionDefault;
422: ksp->ops->buildresidual = KSPBuildResidualDefault;
423: ksp->ops->view = NULL;
425: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG));
426: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG));
427: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG));
428: cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }