1 #include <../src/ksp/ksp/impls/qcg/qcgimpl.h> /*I "petscksp.h" I*/
2
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 */
KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal delta,PetscReal * step1,PetscReal * step2)20 static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
21 {
22 PetscReal dsq, ptp, pts, rad, sts;
23
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 }
39
40 /*@
41 KSPQCGSetTrustRegionRadius - Sets the radius of the trust region for `KSPQCG`
42
43 Logically Collective
44
45 Input Parameters:
46 + ksp - the iterative context
47 - delta - the trust region radius (Infinity is the default)
48
49 Options Database Key:
50 . -ksp_qcg_trustregionradius <delta> - trust region radius
51
52 Level: advanced
53
54 Developer Note:
55 `KSPMINRESSetRadius()`, for example, does not have TrustRegion in the name
56
57 .seealso: [](ch_ksp), `KSPQCG`, `KSPQCGGetTrialStepNorm()`
58 @*/
KSPQCGSetTrustRegionRadius(KSP ksp,PetscReal delta)59 PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
60 {
61 PetscFunctionBegin;
62 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
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 }
67
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`.
71
72 Not Collective
73
74 Input Parameter:
75 . ksp - the iterative context
76
77 Output Parameter:
78 . tsnorm - the norm
79
80 Level: advanced
81
82 .seealso: [](ch_ksp), `KSPQCG`, `KSPQCGSetTrustRegionRadius()`
83 @*/
KSPQCGGetTrialStepNorm(KSP ksp,PetscReal * tsnorm)84 PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
85 {
86 PetscFunctionBegin;
87 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
88 PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
89 PetscFunctionReturn(PETSC_SUCCESS);
90 }
91
92 /*@
93 KSPQCGGetQuadratic - Gets the value of the quadratic function, evaluated at the new iterate
94
95 Collective
96
97 Input Parameter:
98 . ksp - the iterative context
99
100 Output Parameter:
101 . quadratic - the quadratic function evaluated at the new iterate
102
103 Level: advanced
104
105 Note:
106 The quadratic function is
107
108 $$
109 q(s) = g^T * s + 0.5 * s^T * H * s
110 $$
111
112 which satisfies the Euclidean Norm trust region constraint
113
114 $$
115 || D * s || \le delta,
116 $$
117
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
125
126 .seealso: [](ch_ksp), `KSPQCG`
127 @*/
KSPQCGGetQuadratic(KSP ksp,PetscReal * quadratic)128 PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
129 {
130 PetscFunctionBegin;
131 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
132 PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
KSPSolve_QCG(KSP ksp)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 */
144
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 i, maxit;
153 PC pc = ksp->pc;
154 PetscBool diagonalscale;
155
156 PetscFunctionBegin;
157 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
158 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
159 PetscCheck(!ksp->transpose_solve, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Currently does not support transpose solve");
160
161 ksp->its = 0;
162 maxit = ksp->max_it;
163 WA = ksp->work[0];
164 R = ksp->work[1];
165 P = ksp->work[2];
166 ASP = ksp->work[3];
167 BS = ksp->work[4];
168 W = ksp->work[5];
169 WA2 = ksp->work[6];
170 X = ksp->vec_sol;
171 B = ksp->vec_rhs;
172
173 PetscCheck(pcgP->delta > dzero, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: delta <= 0");
174
175 /* Initialize variables */
176 PetscCall(VecSet(W, 0.0)); /* W = 0 */
177 PetscCall(VecSet(X, 0.0)); /* X = 0 */
178 PetscCall(PCGetOperators(pc, &Amat, &Pmat));
179
180 /* Compute: BS = D^{-1} B */
181 PetscCall(PCApplySymmetricLeft(pc, B, BS));
182
183 if (ksp->normtype != KSP_NORM_NONE) {
184 PetscCall(VecNorm(BS, NORM_2, &bsnrm));
185 KSPCheckNorm(ksp, bsnrm);
186 }
187 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
188 ksp->its = 0;
189 ksp->rnorm = bsnrm;
190 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
191 PetscCall(KSPLogResidualHistory(ksp, bsnrm));
192 PetscCall(KSPMonitor(ksp, 0, bsnrm));
193 PetscCall((*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP));
194 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
195
196 /* Compute the initial scaled direction and scaled residual */
197 PetscCall(VecCopy(BS, R));
198 PetscCall(VecScale(R, -1.0));
199 PetscCall(VecCopy(R, P));
200 PetscCall(VecDotRealPart(R, R, &rtr));
201
202 for (i = 0; i <= maxit; i++) {
203 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
204 ksp->its++;
205 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
206
207 /* Compute: asp = D^{-T}*A*D^{-1}*p */
208 PetscCall(PCApplySymmetricRight(pc, P, WA));
209 PetscCall(KSP_MatMult(ksp, Amat, WA, WA2));
210 PetscCall(PCApplySymmetricLeft(pc, WA2, ASP));
211
212 /* Check for negative curvature */
213 PetscCall(VecDotRealPart(P, ASP, &ptasp));
214 if (ptasp <= dzero) {
215 /* Scaled negative curvature direction: Compute a step so that
216 ||w + step*p|| = delta and QS(w + step*p) is least */
217
218 if (!i) {
219 PetscCall(VecCopy(P, X));
220 PetscCall(VecNorm(X, NORM_2, &xnorm));
221 KSPCheckNorm(ksp, xnorm);
222 scal = pcgP->delta / xnorm;
223 PetscCall(VecScale(X, scal));
224 } else {
225 /* Compute roots of quadratic */
226 PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
227 PetscCall(VecDotRealPart(W, ASP, &wtasp));
228 PetscCall(VecDotRealPart(BS, P, &bstp));
229 PetscCall(VecCopy(W, X));
230 q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
231 q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
232 if (q1 <= q2) {
233 PetscCall(VecAXPY(X, step1, P));
234 } else {
235 PetscCall(VecAXPY(X, step2, P));
236 }
237 }
238 pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
239 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
240 if (!i) {
241 PetscCall(PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta));
242 } else {
243 PetscCall(PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
244 }
245
246 } else {
247 /* Compute step along p */
248 step = rtr / ptasp;
249 PetscCall(VecCopy(W, X)); /* x = w */
250 PetscCall(VecAXPY(X, step, P)); /* x <- step*p + x */
251 PetscCall(VecNorm(X, NORM_2, &pcgP->ltsnrm));
252 KSPCheckNorm(ksp, pcgP->ltsnrm);
253
254 if (pcgP->ltsnrm > pcgP->delta) {
255 /* Since the trial iterate is outside the trust region,
256 evaluate a constrained step along p so that
257 ||w + step*p|| = delta
258 The positive step is always better in this case. */
259 if (!i) {
260 scal = pcgP->delta / pcgP->ltsnrm;
261 PetscCall(VecScale(X, scal));
262 } else {
263 /* Compute roots of quadratic */
264 PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
265 PetscCall(VecCopy(W, X));
266 PetscCall(VecAXPY(X, step1, P)); /* x <- step1*p + x */
267 }
268 pcgP->ltsnrm = pcgP->delta;
269 ksp->reason = KSP_CONVERGED_STEP_LENGTH; /* convergence along constrained step */
270 if (!i) {
271 PetscCall(PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta));
272 } else {
273 PetscCall(PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
274 }
275
276 } else {
277 /* Evaluate the current step */
278 PetscCall(VecCopy(X, W)); /* update interior iterate */
279 PetscCall(VecAXPY(R, -step, ASP)); /* r <- -step*asp + r */
280 if (ksp->normtype != KSP_NORM_NONE) {
281 PetscCall(VecNorm(R, NORM_2, &rnrm));
282 KSPCheckNorm(ksp, rnrm);
283 }
284 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
285 ksp->rnorm = rnrm;
286 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
287 PetscCall(KSPLogResidualHistory(ksp, rnrm));
288 PetscCall(KSPMonitor(ksp, i + 1, rnrm));
289 PetscCall((*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP));
290 if (ksp->reason) { /* convergence for */
291 PetscCall(PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta));
292 }
293 }
294 }
295 if (ksp->reason) break; /* Convergence has been attained */
296 else { /* Compute a new AS-orthogonal direction */ PetscCall(VecDot(R, R, &rntrn));
297 beta = rntrn / rtr;
298 PetscCall(VecAYPX(P, beta, R)); /* p <- r + beta*p */
299 rtr = PetscRealPart(rntrn);
300 }
301 }
302 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
303
304 /* Unscale x */
305 PetscCall(VecCopy(X, WA2));
306 PetscCall(PCApplySymmetricRight(pc, WA2, X));
307
308 PetscCall(KSP_MatMult(ksp, Amat, X, WA));
309 PetscCall(VecDotRealPart(B, X, &btx));
310 PetscCall(VecDotRealPart(X, WA, &xtax));
311
312 pcgP->quadratic = btx + .5 * xtax;
313 PetscFunctionReturn(PETSC_SUCCESS);
314 }
315
KSPSetUp_QCG(KSP ksp)316 static PetscErrorCode KSPSetUp_QCG(KSP ksp)
317 {
318 PetscFunctionBegin;
319 /* Get work vectors from user code */
320 PetscCall(KSPSetWorkVecs(ksp, 7));
321 PetscFunctionReturn(PETSC_SUCCESS);
322 }
323
KSPDestroy_QCG(KSP ksp)324 static PetscErrorCode KSPDestroy_QCG(KSP ksp)
325 {
326 PetscFunctionBegin;
327 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL));
328 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL));
329 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL));
330 PetscCall(KSPDestroyDefault(ksp));
331 PetscFunctionReturn(PETSC_SUCCESS);
332 }
333
KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)334 static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
335 {
336 KSP_QCG *cgP = (KSP_QCG *)ksp->data;
337
338 PetscFunctionBegin;
339 cgP->delta = delta;
340 PetscFunctionReturn(PETSC_SUCCESS);
341 }
342
KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal * ltsnrm)343 static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
344 {
345 KSP_QCG *cgP = (KSP_QCG *)ksp->data;
346
347 PetscFunctionBegin;
348 *ltsnrm = cgP->ltsnrm;
349 PetscFunctionReturn(PETSC_SUCCESS);
350 }
351
KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal * quadratic)352 static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
353 {
354 KSP_QCG *cgP = (KSP_QCG *)ksp->data;
355
356 PetscFunctionBegin;
357 *quadratic = cgP->quadratic;
358 PetscFunctionReturn(PETSC_SUCCESS);
359 }
360
KSPSetFromOptions_QCG(KSP ksp,PetscOptionItems PetscOptionsObject)361 static PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems PetscOptionsObject)
362 {
363 PetscReal delta;
364 KSP_QCG *cgP = (KSP_QCG *)ksp->data;
365 PetscBool flg;
366
367 PetscFunctionBegin;
368 PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
369 PetscCall(PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg));
370 if (flg) PetscCall(KSPQCGSetTrustRegionRadius(ksp, delta));
371 PetscOptionsHeadEnd();
372 PetscFunctionReturn(PETSC_SUCCESS);
373 }
374
375 /*MC
376 KSPQCG - Code to run conjugate gradient method subject to a constraint on the solution norm {cite}`steihaug:83`.
377
378 Options Database Key:
379 . -ksp_qcg_trustregionradius <r> - Trust Region Radius
380
381 Level: developer
382
383 Notes:
384 This is rarely used directly, ir is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
385
386 Uses preconditioned conjugate gradient to compute
387 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
388 $ || D * s || \le delta$, where
389 .vb
390 delta is the trust region radius,
391 g is the gradient vector, and
392 H is Hessian matrix,
393 D is a scaling matrix.
394 .ve
395
396 `KSPConvergedReason` may include
397 + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
398 - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,
399
400 Note:
401 Allows symmetric preconditioning with the following scaling matrices:
402 .vb
403 `PCNONE`: D = Identity matrix
404 `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
405 `PCICC`: D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
406 .ve
407
408 .seealso: [](ch_ksp), `KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`
409 `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
410 M*/
411
KSPCreate_QCG(KSP ksp)412 PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
413 {
414 KSP_QCG *cgP;
415
416 PetscFunctionBegin;
417 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3));
418 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1));
419 PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
420 PetscCall(PetscNew(&cgP));
421
422 ksp->data = (void *)cgP;
423 ksp->ops->setup = KSPSetUp_QCG;
424 ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
425 ksp->ops->solve = KSPSolve_QCG;
426 ksp->ops->destroy = KSPDestroy_QCG;
427 ksp->ops->buildsolution = KSPBuildSolutionDefault;
428 ksp->ops->buildresidual = KSPBuildResidualDefault;
429 ksp->ops->view = NULL;
430
431 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG));
432 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG));
433 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG));
434 cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
435 PetscFunctionReturn(PETSC_SUCCESS);
436 }
437