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: .seealso: `KSPQCG`
55: @*/
56: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
57: {
58: PetscFunctionBegin;
60: PetscCheck(delta >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
61: PetscTryMethod(ksp, "KSPQCGSetTrustRegionRadius_C", (KSP, PetscReal), (ksp, delta));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*@
66: KSPQCGGetTrialStepNorm - Gets the norm of a trial step vector in `KSPQCG`. The WCG step may be
67: constrained, so this is not necessarily the length of the ultimate step taken in `KSPQCG`.
69: Not Collective
71: Input Parameter:
72: . ksp - the iterative context
74: Output Parameter:
75: . tsnorm - the norm
77: Level: advanced
79: .seealso: `KSPQCG`
80: @*/
81: PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
82: {
83: PetscFunctionBegin;
85: PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*@
90: KSPQCGGetQuadratic - Gets the value of the quadratic function, evaluated at the new iterate
92: Collective
94: Input Parameter:
95: . ksp - the iterative context
97: Output Parameter:
98: . quadratic - the quadratic function evaluated at the new iterate
100: Level: advanced
102: Note:
103: $ q(s) = g^T * s + 0.5 * s^T * H * s $ which satisfies the Euclidean Norm trust region constraint
104: .vb
105: || D * s || <= delta,
106: .ve
107: where
108: .vb
109: delta is the trust region radius,
110: g is the gradient vector, and
111: H is Hessian matrix,
112: D is a scaling matrix.
113: .ve
115: .seealso: [](ch_ksp), `KSPQCG`
116: @*/
117: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
118: {
119: PetscFunctionBegin;
121: PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode KSPSolve_QCG(KSP ksp)
126: {
127: /*
128: Correspondence with documentation above:
129: B = g = gradient,
130: X = s = step
131: Note: This is not coded correctly for complex arithmetic!
132: */
134: KSP_QCG *pcgP = (KSP_QCG *)ksp->data;
135: Mat Amat, Pmat;
136: Vec W, WA, WA2, R, P, ASP, BS, X, B;
137: PetscScalar scal, beta, rntrn, step;
138: PetscReal q1, q2, xnorm, step1, step2, rnrm = 0.0, btx, xtax;
139: PetscReal ptasp, rtr, wtasp, bstp;
140: PetscReal dzero = 0.0, bsnrm = 0.0;
141: PetscInt i, maxit;
142: PC pc = ksp->pc;
143: PetscBool diagonalscale;
145: PetscFunctionBegin;
146: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
147: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
148: PetscCheck(!ksp->transpose_solve, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Currently does not support transpose solve");
150: ksp->its = 0;
151: maxit = ksp->max_it;
152: WA = ksp->work[0];
153: R = ksp->work[1];
154: P = ksp->work[2];
155: ASP = ksp->work[3];
156: BS = ksp->work[4];
157: W = ksp->work[5];
158: WA2 = ksp->work[6];
159: X = ksp->vec_sol;
160: B = ksp->vec_rhs;
162: PetscCheck(pcgP->delta > dzero, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: delta <= 0");
164: /* Initialize variables */
165: PetscCall(VecSet(W, 0.0)); /* W = 0 */
166: PetscCall(VecSet(X, 0.0)); /* X = 0 */
167: PetscCall(PCGetOperators(pc, &Amat, &Pmat));
169: /* Compute: BS = D^{-1} B */
170: PetscCall(PCApplySymmetricLeft(pc, B, BS));
172: if (ksp->normtype != KSP_NORM_NONE) {
173: PetscCall(VecNorm(BS, NORM_2, &bsnrm));
174: KSPCheckNorm(ksp, bsnrm);
175: }
176: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
177: ksp->its = 0;
178: ksp->rnorm = bsnrm;
179: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
180: PetscCall(KSPLogResidualHistory(ksp, bsnrm));
181: PetscCall(KSPMonitor(ksp, 0, bsnrm));
182: PetscCall((*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP));
183: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
185: /* Compute the initial scaled direction and scaled residual */
186: PetscCall(VecCopy(BS, R));
187: PetscCall(VecScale(R, -1.0));
188: PetscCall(VecCopy(R, P));
189: PetscCall(VecDotRealPart(R, R, &rtr));
191: for (i = 0; i <= maxit; i++) {
192: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
193: ksp->its++;
194: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
196: /* Compute: asp = D^{-T}*A*D^{-1}*p */
197: PetscCall(PCApplySymmetricRight(pc, P, WA));
198: PetscCall(KSP_MatMult(ksp, Amat, WA, WA2));
199: PetscCall(PCApplySymmetricLeft(pc, WA2, ASP));
201: /* Check for negative curvature */
202: PetscCall(VecDotRealPart(P, ASP, &ptasp));
203: if (ptasp <= dzero) {
204: /* Scaled negative curvature direction: Compute a step so that
205: ||w + step*p|| = delta and QS(w + step*p) is least */
207: if (!i) {
208: PetscCall(VecCopy(P, X));
209: PetscCall(VecNorm(X, NORM_2, &xnorm));
210: KSPCheckNorm(ksp, xnorm);
211: scal = pcgP->delta / xnorm;
212: PetscCall(VecScale(X, scal));
213: } else {
214: /* Compute roots of quadratic */
215: PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
216: PetscCall(VecDotRealPart(W, ASP, &wtasp));
217: PetscCall(VecDotRealPart(BS, P, &bstp));
218: PetscCall(VecCopy(W, X));
219: q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
220: q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
221: if (q1 <= q2) {
222: PetscCall(VecAXPY(X, step1, P));
223: } else {
224: PetscCall(VecAXPY(X, step2, P));
225: }
226: }
227: pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
228: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
229: if (!i) {
230: PetscCall(PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta));
231: } else {
232: PetscCall(PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
233: }
235: } else {
236: /* Compute step along p */
237: step = rtr / ptasp;
238: PetscCall(VecCopy(W, X)); /* x = w */
239: PetscCall(VecAXPY(X, step, P)); /* x <- step*p + x */
240: PetscCall(VecNorm(X, NORM_2, &pcgP->ltsnrm));
241: KSPCheckNorm(ksp, pcgP->ltsnrm);
243: if (pcgP->ltsnrm > pcgP->delta) {
244: /* Since the trial iterate is outside the trust region,
245: evaluate a constrained step along p so that
246: ||w + step*p|| = delta
247: The positive step is always better in this case. */
248: if (!i) {
249: scal = pcgP->delta / pcgP->ltsnrm;
250: PetscCall(VecScale(X, scal));
251: } else {
252: /* Compute roots of quadratic */
253: PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
254: PetscCall(VecCopy(W, X));
255: PetscCall(VecAXPY(X, step1, P)); /* x <- step1*p + x */
256: }
257: pcgP->ltsnrm = pcgP->delta;
258: ksp->reason = KSP_CONVERGED_STEP_LENGTH; /* convergence along constrained step */
259: if (!i) {
260: PetscCall(PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta));
261: } else {
262: PetscCall(PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
263: }
265: } else {
266: /* Evaluate the current step */
267: PetscCall(VecCopy(X, W)); /* update interior iterate */
268: PetscCall(VecAXPY(R, -step, ASP)); /* r <- -step*asp + r */
269: if (ksp->normtype != KSP_NORM_NONE) {
270: PetscCall(VecNorm(R, NORM_2, &rnrm));
271: KSPCheckNorm(ksp, rnrm);
272: }
273: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
274: ksp->rnorm = rnrm;
275: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
276: PetscCall(KSPLogResidualHistory(ksp, rnrm));
277: PetscCall(KSPMonitor(ksp, i + 1, rnrm));
278: PetscCall((*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP));
279: if (ksp->reason) { /* convergence for */
280: PetscCall(PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta));
281: }
282: }
283: }
284: if (ksp->reason) break; /* Convergence has been attained */
285: else { /* Compute a new AS-orthogonal direction */ PetscCall(VecDot(R, R, &rntrn));
286: beta = rntrn / rtr;
287: PetscCall(VecAYPX(P, beta, R)); /* p <- r + beta*p */
288: rtr = PetscRealPart(rntrn);
289: }
290: }
291: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
293: /* Unscale x */
294: PetscCall(VecCopy(X, WA2));
295: PetscCall(PCApplySymmetricRight(pc, WA2, X));
297: PetscCall(KSP_MatMult(ksp, Amat, X, WA));
298: PetscCall(VecDotRealPart(B, X, &btx));
299: PetscCall(VecDotRealPart(X, WA, &xtax));
301: pcgP->quadratic = btx + .5 * xtax;
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode KSPSetUp_QCG(KSP ksp)
306: {
307: PetscFunctionBegin;
308: /* Get work vectors from user code */
309: PetscCall(KSPSetWorkVecs(ksp, 7));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode KSPDestroy_QCG(KSP ksp)
314: {
315: PetscFunctionBegin;
316: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL));
317: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL));
318: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL));
319: PetscCall(KSPDestroyDefault(ksp));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
324: {
325: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
327: PetscFunctionBegin;
328: cgP->delta = delta;
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
333: {
334: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
336: PetscFunctionBegin;
337: *ltsnrm = cgP->ltsnrm;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
342: {
343: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
345: PetscFunctionBegin;
346: *quadratic = cgP->quadratic;
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: static PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
351: {
352: PetscReal delta;
353: KSP_QCG *cgP = (KSP_QCG *)ksp->data;
354: PetscBool flg;
356: PetscFunctionBegin;
357: PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
358: PetscCall(PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg));
359: if (flg) PetscCall(KSPQCGSetTrustRegionRadius(ksp, delta));
360: PetscOptionsHeadEnd();
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: /*MC
365: KSPQCG - Code to run conjugate gradient method subject to a constraint on the solution norm.
367: Options Database Key:
368: . -ksp_qcg_trustregionradius <r> - Trust Region Radius
370: Level: developer
372: Notes:
373: This is rarely used directly, ir is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
375: Uses preconditioned conjugate gradient to compute
376: 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
377: .vb
378: || D * s || <= delta,
379: .ve
380: where
381: .vb
382: delta is the trust region radius,
383: g is the gradient vector, and
384: H is Hessian matrix,
385: D is a scaling matrix.
386: .ve
388: `KSPConvergedReason` may be
389: .vb
390: KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
391: KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step,
392: .ve
393: and other `KSP` converged/diverged reasons
395: Notes:
396: Currently we allow symmetric preconditioning with the following scaling matrices:
397: .vb
398: `PCNONE`: D = Identity matrix
399: `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
400: `PCICC`: D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
401: .ve
403: References:
404: . * - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
405: SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).
407: .seealso: [](ch_ksp), 'KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`
408: `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
409: M*/
411: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
412: {
413: KSP_QCG *cgP;
415: PetscFunctionBegin;
416: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3));
417: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1));
418: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
419: PetscCall(PetscNew(&cgP));
421: ksp->data = (void *)cgP;
422: ksp->ops->setup = KSPSetUp_QCG;
423: ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
424: ksp->ops->solve = KSPSolve_QCG;
425: ksp->ops->destroy = KSPDestroy_QCG;
426: ksp->ops->buildsolution = KSPBuildSolutionDefault;
427: ksp->ops->buildresidual = KSPBuildResidualDefault;
428: ksp->ops->view = NULL;
430: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG));
431: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG));
432: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG));
433: cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }