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