1 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
2 #include <petscblaslapack.h>
3 PETSC_INTERN PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP, PetscReal *, PetscReal *);
4 PETSC_INTERN PetscErrorCode KSPComputeEigenvalues_MINRES(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
5
6 PetscBool QLPcite = PETSC_FALSE;
7 const char QLPCitation[] = "@article{choi2011minres,\n"
8 " title={MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems},\n"
9 " author={Choi, Sou-Cheng T and Paige, Christopher C and Saunders, Michael A},\n"
10 " journal={SIAM Journal on Scientific Computing},\n"
11 " volume={33},\n"
12 " number={4},\n"
13 " pages={1810--1836},\n"
14 " year={2011},\n}\n";
15
16 typedef struct {
17 PetscReal haptol;
18 PetscReal nutol;
19 PetscBool qlp;
20 PetscReal maxxnorm;
21 PetscReal TranCond;
22 PetscBool monitor;
23 PetscViewer viewer;
24 PetscViewerFormat viewer_fmt;
25 // The following arrays are of size ksp->maxit
26 PetscScalar *e, *d;
27 PetscReal *ee, *dd; /* work space for Lanczos algorithm */
28 } KSP_MINRES;
29
KSPSetUp_MINRES(KSP ksp)30 static PetscErrorCode KSPSetUp_MINRES(KSP ksp)
31 {
32 PetscFunctionBegin;
33 PetscCall(KSPSetWorkVecs(ksp, 9));
34 /*
35 If user requested computations of eigenvalues then allocate
36 work space needed
37 */
38 if (ksp->calc_sings) {
39 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
40 PetscInt maxit = ksp->max_it;
41 PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
42 PetscCall(PetscMalloc4(maxit + 1, &minres->e, maxit, &minres->d, maxit, &minres->ee, maxit, &minres->dd));
43
44 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_MINRES;
45 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_MINRES;
46 }
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
50 /* Convenience functions */
51 #define KSPMinresSwap3(V1, V2, V3) \
52 do { \
53 Vec T = V1; \
54 V1 = V2; \
55 V2 = V3; \
56 V3 = T; \
57 } while (0)
58
Norm3(PetscReal a,PetscReal b,PetscReal c)59 static inline PetscReal Norm3(PetscReal a, PetscReal b, PetscReal c)
60 {
61 return PetscSqrtReal(PetscSqr(a) + PetscSqr(b) + PetscSqr(c));
62 }
63
SymOrtho(PetscReal a,PetscReal b,PetscReal * c,PetscReal * s,PetscReal * r)64 static inline void SymOrtho(PetscReal a, PetscReal b, PetscReal *c, PetscReal *s, PetscReal *r)
65 {
66 if (b == 0.0) {
67 if (a == 0.0) *c = 1.0;
68 else *c = PetscCopysignReal(1.0, a);
69 *s = 0.0;
70 *r = PetscAbsReal(a);
71 } else if (a == 0.0) {
72 *c = 0.0;
73 *s = PetscCopysignReal(1.0, b);
74 *r = PetscAbsReal(b);
75 } else if (PetscAbsReal(b) > PetscAbsReal(a)) {
76 PetscReal t = a / b;
77
78 *s = PetscCopysignReal(1.0, b) / PetscSqrtReal(1.0 + t * t);
79 *c = (*s) * t;
80 *r = b / (*s); // computationally better than d = a / c since |c| <= |s|
81 } else {
82 PetscReal t = b / a;
83
84 *c = PetscCopysignReal(1.0, a) / PetscSqrtReal(1.0 + t * t);
85 *s = (*c) * t;
86 *r = a / (*c); // computationally better than d = b / s since |s| <= |c|
87 }
88 }
89
90 /*
91 Code adapted from https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
92 CSP11/Algorithms/MINRESQLP/minresQLP.m
93 */
KSPSolve_MINRES(KSP ksp)94 static PetscErrorCode KSPSolve_MINRES(KSP ksp)
95 {
96 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
97 Mat Amat;
98 Vec X, B, R1, R2, R3, V, W, WL, WL2, XL2, RN;
99 PetscReal alpha, beta, beta1, betan, betal;
100 PetscBool diagonalscale;
101 PetscReal zero = 0.0, dbar, dltan = 0.0, dlta, cs = -1.0, sn = 0.0, epln, eplnn = 0.0, gbar, dlta_QLP;
102 PetscReal gamal3 = 0.0, gamal2 = 0.0, gamal = 0.0, gama = 0.0, gama_tmp;
103 PetscReal taul2 = 0.0, taul = 0.0, tau = 0.0, phi, phi0, phir;
104 PetscReal Axnorm, xnorm, xnorm_tmp, xl2norm = 0.0, pnorm, Anorm = 0.0, gmin = 0.0, gminl = 0.0, gminl2 = 0.0;
105 PetscReal Acond = 1.0, Acondl = 0.0, rnorml, rnorm, rootl, relAresl, relres, relresl, Arnorml, Anorml = 0.0;
106 PetscReal epsx, realmin = PETSC_REAL_MIN, eps = PETSC_MACHINE_EPSILON;
107 PetscReal veplnl2 = 0.0, veplnl = 0.0, vepln = 0.0, etal2 = 0.0, etal = 0.0, eta = 0.0;
108 PetscReal dlta_tmp, sr2 = 0.0, cr2 = -1.0, cr1 = -1.0, sr1 = 0.0;
109 PetscReal ul4 = 0.0, ul3 = 0.0, ul2 = 0.0, ul = 0.0, u = 0.0, ul_QLP = 0.0, u_QLP = 0.0;
110 PetscReal vepln_QLP = 0.0, gamal_QLP = 0.0, gama_QLP = 0.0, gamal_tmp, abs_gama;
111 PetscInt flag = -2, flag0 = -2, QLPiter = 0;
112 PetscInt stored_max_it, eigs;
113 PetscScalar *e = NULL, *d = NULL;
114
115 PetscFunctionBegin;
116 PetscCall(PetscCitationsRegister(QLPCitation, &QLPcite));
117 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
118 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
119
120 eigs = ksp->calc_sings;
121 stored_max_it = ksp->max_it;
122 if (eigs) {
123 e = minres->e;
124 d = minres->d;
125 }
126
127 X = ksp->vec_sol;
128 B = ksp->vec_rhs;
129 R1 = ksp->work[0];
130 R2 = ksp->work[1];
131 R3 = ksp->work[2];
132 V = ksp->work[3];
133 W = ksp->work[4];
134 WL = ksp->work[5];
135 WL2 = ksp->work[6];
136 XL2 = ksp->work[7];
137 RN = ksp->work[8];
138 PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
139
140 ksp->its = 0;
141 ksp->rnorm = 0.0;
142 if (!ksp->guess_zero) {
143 PetscCall(KSP_MatMult(ksp, Amat, X, R2));
144 PetscCall(VecNorm(R2, NORM_2, &Axnorm));
145 PetscCall(VecNorm(X, NORM_2, &xnorm));
146 PetscCall(VecAYPX(R2, -1.0, B));
147 } else {
148 PetscCall(VecCopy(B, R2));
149 Axnorm = 0.0;
150 xnorm = 0.0;
151 }
152 PetscCall(KSP_PCApply(ksp, R2, R3));
153 if (ksp->converged_neg_curve) PetscCall(VecCopy(R3, RN));
154 PetscCall(VecDotRealPart(R3, R2, &beta1));
155 KSPCheckDot(ksp, beta1);
156 if (beta1 < 0.0) {
157 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g", (double)beta1);
158 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
159 PetscFunctionReturn(PETSC_SUCCESS);
160 }
161 beta1 = PetscSqrtReal(beta1);
162
163 rnorm = beta1;
164 if (ksp->normtype == KSP_NORM_PRECONDITIONED) ksp->rnorm = rnorm;
165 else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) PetscCall(VecNorm(R2, NORM_2, &ksp->rnorm));
166 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
167 PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
168 PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
169 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
170
171 relres = rnorm / beta1;
172 betan = beta1;
173 phi0 = beta1;
174 phi = beta1;
175 betan = beta1;
176 beta = 0.0;
177 do {
178 /* Lanczos */
179 ksp->its++;
180 betal = beta;
181 beta = betan;
182 PetscCall(VecAXPBY(V, 1.0 / beta, 0.0, R3));
183 PetscCall(KSP_MatMult(ksp, Amat, V, R3));
184 if (ksp->its > 1) PetscCall(VecAXPY(R3, -beta / betal, R1));
185 PetscCall(VecDotRealPart(R3, V, &alpha));
186 PetscCall(VecAXPY(R3, -alpha / beta, R2));
187 KSPMinresSwap3(R1, R2, R3);
188 if (eigs) {
189 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
190 d[ksp->its - 1] = alpha;
191 e[ksp->its - 1] = beta;
192 }
193
194 PetscCall(KSP_PCApply(ksp, R2, R3));
195 PetscCall(VecDotRealPart(R3, R2, &betan));
196 KSPCheckDot(ksp, betan);
197 if (betan < 0.0) {
198 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite preconditioner %g", (double)betan);
199 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202 betan = PetscSqrtReal(betan);
203
204 pnorm = Norm3(betal, alpha, betan);
205
206 // Apply previous left rotation Q_{k-1}
207 dbar = dltan;
208 epln = eplnn;
209 dlta = cs * dbar + sn * alpha;
210 gbar = sn * dbar - cs * alpha;
211 eplnn = sn * betan;
212 dltan = -cs * betan;
213 dlta_QLP = dlta;
214
215 // Stop if negative curvature is detected and return residual
216 // This is very experimental and maybe changed in the future
217 // based on https://arxiv.org/pdf/2208.07095.pdf
218 if (ksp->converged_neg_curve) {
219 if (cs * gbar >= 0.0) {
220 PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1 %g, gbar %g\n", (double)cs, (double)gbar));
221 ksp->reason = KSP_CONVERGED_NEG_CURVE;
222 PetscCall(VecCopy(RN, X));
223 break;
224 } else {
225 PetscCall(VecAXPBY(RN, -phi * cs, PetscSqr(sn), V));
226 }
227 }
228
229 // Compute the current left plane rotation Q_k
230 gamal3 = gamal2;
231 gamal2 = gamal;
232 gamal = gama;
233 SymOrtho(gbar, betan, &cs, &sn, &gama);
234
235 // Inexactness condition from https://arxiv.org/pdf/2208.07095.pdf
236 rootl = Norm3(gbar, dltan, zero);
237 phir = PetscSqr(phi0 / phi);
238 if (ksp->its > 2 && minres->nutol > 0.0) {
239 PetscReal tmp;
240
241 phir = PetscSqrtReal(phir - 1.0);
242 tmp = rootl / phir;
243 PetscCall(PetscInfo(ksp, "it = %" PetscInt_FMT ": inexact check %g (%g / %g)\n", ksp->its - 2, (double)tmp, (double)rootl, (double)phir));
244 if (tmp < minres->nutol) {
245 ksp->its--;
246 ksp->reason = KSP_CONVERGED_RTOL;
247 break;
248 }
249 }
250
251 gama_tmp = gama;
252 taul2 = taul;
253 taul = tau;
254 tau = cs * phi;
255 Axnorm = Norm3(Axnorm, tau, zero);
256 phi = sn * phi;
257
258 //Apply the previous right plane rotation P{k-2,k}
259 if (ksp->its > 2) {
260 veplnl2 = veplnl;
261 etal2 = etal;
262 etal = eta;
263 dlta_tmp = sr2 * vepln - cr2 * dlta;
264 veplnl = cr2 * vepln + sr2 * dlta;
265 dlta = dlta_tmp;
266 eta = sr2 * gama;
267 gama = -cr2 * gama;
268 }
269
270 // Compute the current right plane rotation P{k-1,k}, P_12, P_23,...
271 if (ksp->its > 1) {
272 SymOrtho(gamal, dlta, &cr1, &sr1, &gamal);
273 vepln = sr1 * gama;
274 gama = -cr1 * gama;
275 }
276
277 // Update xnorm
278 ul4 = ul3;
279 ul3 = ul2;
280 if (ksp->its > 2) ul2 = (taul2 - etal2 * ul4 - veplnl2 * ul3) / gamal2;
281 if (ksp->its > 1) ul = (taul - etal * ul3 - veplnl * ul2) / gamal;
282 xnorm_tmp = Norm3(xl2norm, ul2, ul);
283 if (PetscAbsReal(gama) > realmin && xnorm_tmp < minres->maxxnorm) {
284 u = (tau - eta * ul2 - vepln * ul) / gama;
285 if (Norm3(xnorm_tmp, u, zero) > minres->maxxnorm) {
286 u = 0;
287 flag = 6;
288 }
289 } else {
290 u = 0;
291 flag = 9;
292 }
293 xl2norm = Norm3(xl2norm, ul2, zero);
294 xnorm = Norm3(xl2norm, ul, u);
295
296 // Update w. Update x except if it will become too big
297 //if (Acond < minres->TranCond && flag != flag0 && QLPiter == 0) { // I believe they have a typo in the MATLAB code
298 if ((Acond < minres->TranCond || !minres->qlp) && flag == flag0 && QLPiter == 0) { // MINRES
299 KSPMinresSwap3(WL2, WL, W);
300 PetscCall(VecAXPBY(W, 1.0 / gama_tmp, 0.0, V));
301 if (ksp->its > 1) {
302 Vec T[] = {WL, WL2};
303 PetscScalar alphas[] = {-dlta_QLP / gama_tmp, -epln / gama_tmp};
304 PetscInt nv = (ksp->its == 2 ? 1 : 2);
305
306 PetscCall(VecMAXPY(W, nv, alphas, T));
307 }
308 if (xnorm < minres->maxxnorm) {
309 PetscCall(VecAXPY(X, tau, W));
310 } else {
311 flag = 6;
312 }
313 } else if (minres->qlp) { //MINRES-QLP updates
314 QLPiter = QLPiter + 1;
315 if (QLPiter == 1) {
316 // xl2 = x - wl*ul_QLP - w*u_QLP;
317 PetscScalar maxpys[] = {1.0, -ul_QLP, -u_QLP};
318 Vec maxpyv[] = {X, WL, W};
319
320 PetscCall(VecSet(XL2, 0.0));
321 // construct w_{k-3}, w_{k-2}, w_{k-1}
322 if (ksp->its > 1) {
323 if (ksp->its > 3) { // w_{k-3}
324 //wl2 = gamal3*wl2 + veplnl2*wl + etal*w;
325 PetscCall(VecAXPBYPCZ(WL2, veplnl2, etal, gamal3, WL, W));
326 }
327 if (ksp->its > 2) { // w_{k-2}
328 //wl = gamal_QLP*wl + vepln_QLP*w;
329 PetscCall(VecAXPBY(WL, vepln_QLP, gamal_QLP, W));
330 }
331 // w = gama_QLP*w;
332 PetscCall(VecScale(W, gama_QLP));
333 // xl2 = x - wl*ul_QLP - w*u_QLP;
334 PetscCall(VecMAXPY(XL2, 3, maxpys, maxpyv));
335 }
336 }
337 if (ksp->its == 1) {
338 //wl2 = wl; wl = v*sr1; w = -v*cr1;
339 PetscCall(VecCopy(WL, WL2));
340 PetscCall(VecAXPBY(WL, sr1, 0, V));
341 PetscCall(VecAXPBY(W, -cr1, 0, V));
342 } else if (ksp->its == 2) {
343 //wl2 = wl;
344 //wl = w*cr1 + v*sr1;
345 //w = w*sr1 - v*cr1;
346 PetscCall(VecCopy(WL, WL2));
347 PetscCall(VecAXPBYPCZ(WL, cr1, sr1, 0.0, W, V));
348 PetscCall(VecAXPBY(W, -cr1, sr1, V));
349 } else {
350 //wl2 = wl; wl = w; w = wl2*sr2 - v*cr2;
351 //wl2 = wl2*cr2 + v*sr2; v = wl *cr1 + w*sr1;
352 //w = wl *sr1 - w*cr1; wl = v;
353 PetscCall(VecCopy(WL, WL2));
354 PetscCall(VecCopy(W, WL));
355 PetscCall(VecAXPBYPCZ(W, sr2, -cr2, 0.0, WL2, V));
356 PetscCall(VecAXPBY(WL2, sr2, cr2, V));
357 PetscCall(VecAXPBYPCZ(V, cr1, sr1, 0.0, WL, W));
358 PetscCall(VecAXPBY(W, sr1, -cr1, WL));
359 PetscCall(VecCopy(V, WL));
360 }
361
362 //xl2 = xl2 + wl2*ul2;
363 PetscCall(VecAXPY(XL2, ul2, WL2));
364 //x = xl2 + wl *ul + w*u;
365 PetscCall(VecCopy(XL2, X));
366 PetscCall(VecAXPBYPCZ(X, ul, u, 1.0, WL, W));
367 }
368 // Compute the next right plane rotation P{k-1,k+1}
369 gamal_tmp = gamal;
370 SymOrtho(gamal, eplnn, &cr2, &sr2, &gamal);
371
372 //Store quantities for transferring from MINRES to MINRES-QLP
373 gamal_QLP = gamal_tmp;
374 vepln_QLP = vepln;
375 gama_QLP = gama;
376 ul_QLP = ul;
377 u_QLP = u;
378
379 // Estimate various norms
380 abs_gama = PetscAbsReal(gama);
381 Anorml = Anorm;
382 Anorm = PetscMax(PetscMax(Anorm, pnorm), PetscMax(gamal, abs_gama));
383 if (ksp->its == 1) {
384 gmin = gama;
385 gminl = gmin;
386 } else {
387 gminl2 = gminl;
388 gminl = gmin;
389 gmin = PetscMin(gminl2, PetscMin(gamal, abs_gama));
390 }
391 Acondl = Acond;
392 Acond = Anorm / gmin;
393 rnorml = rnorm;
394 relresl = relres;
395 if (flag != 9) rnorm = phi;
396 relres = rnorm / (Anorm * xnorm + beta1);
397 Arnorml = rnorml * rootl;
398 relAresl = rootl / Anorm;
399
400 // See if any of the stopping criteria are satisfied.
401 epsx = Anorm * xnorm * eps;
402 if (flag == flag0 || flag == 9) {
403 //if (Acond >= Acondlim) flag = 7; // Huge Acond
404 if (epsx >= beta1) flag = 5; // x is an eigenvector
405 if (minres->qlp) { /* We use these indicators only if the QLP variant has been selected */
406 PetscReal t1 = 1.0 + relres;
407 PetscReal t2 = 1.0 + relAresl;
408 if (xnorm >= minres->maxxnorm) flag = 6; // xnorm exceeded its limit
409 if (t2 <= 1) flag = 4; // Accurate LS solution
410 if (t1 <= 1) flag = 3; // Accurate Ax=b solution
411 if (relAresl <= ksp->rtol) flag = 2; // Good enough LS solution
412 if (relres <= ksp->rtol) flag = 1; // Good enough Ax=b solution
413 }
414 }
415
416 if (flag == 2 || flag == 4 || flag == 6 || flag == 7) {
417 Acond = Acondl;
418 rnorm = rnorml;
419 relres = relresl;
420 }
421
422 if (minres->monitor) { /* Mimics MATLAB code with extra flag */
423 PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
424 if (ksp->its == 1) PetscCall(PetscViewerASCIIPrintf(minres->viewer, " flag rnorm Arnorm Compatible LS Anorm Acond xnorm\n"));
425 PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5" PetscInt_FMT " %2" PetscInt_FMT " %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", ksp->its - 1, flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorm));
426 PetscCall(PetscViewerPopFormat(minres->viewer));
427 }
428
429 if (ksp->normtype == KSP_NORM_PRECONDITIONED) ksp->rnorm = rnorm;
430 else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
431 PetscCall(KSP_MatMult(ksp, Amat, X, V));
432 PetscCall(VecAYPX(V, -1.0, B));
433 PetscCall(VecNorm(V, NORM_2, &ksp->rnorm));
434 }
435 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
436 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
437 PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
438 if (!ksp->reason) {
439 switch (flag) {
440 case 1:
441 case 2:
442 case 5: /* XXX */
443 ksp->reason = KSP_CONVERGED_RTOL;
444 break;
445 case 3:
446 case 4:
447 ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
448 break;
449 case 6:
450 ksp->reason = KSP_CONVERGED_STEP_LENGTH;
451 break;
452 default:
453 break;
454 }
455 }
456 if (ksp->reason) break;
457 } while (ksp->its < ksp->max_it);
458
459 if (minres->monitor && flag != 2 && flag != 4 && flag != 6 && flag != 7) {
460 PetscCall(VecNorm(X, NORM_2, &xnorm));
461 PetscCall(KSP_MatMult(ksp, Amat, X, R1));
462 PetscCall(VecAYPX(R1, -1.0, B));
463 PetscCall(VecNorm(R1, NORM_2, &rnorml));
464 PetscCall(KSP_MatMult(ksp, Amat, R1, R2));
465 PetscCall(VecNorm(R2, NORM_2, &Arnorml));
466 relresl = rnorml / (Anorm * xnorm + beta1);
467 relAresl = rnorml > realmin ? Arnorml / (Anorm * rnorml) : 0.0;
468 PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
469 PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5" PetscInt_FMT " %2" PetscInt_FMT " %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", ksp->its, flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorm));
470 PetscCall(PetscViewerPopFormat(minres->viewer));
471 }
472 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
473 PetscFunctionReturn(PETSC_SUCCESS);
474 }
475
476 /* This was the original implementation provided by R. Scheichl */
KSPSolve_MINRES_OLD(KSP ksp)477 static PetscErrorCode KSPSolve_MINRES_OLD(KSP ksp)
478 {
479 PetscInt i;
480 PetscScalar alpha, beta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
481 PetscScalar rho0, rho1, rho2, rho3, dp = 0.0;
482 const PetscScalar none = -1.0;
483 PetscReal np;
484 Vec X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
485 Mat Amat;
486 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
487 PetscBool diagonalscale;
488 PetscInt stored_max_it, eigs;
489 PetscScalar *e = NULL, *d = NULL;
490
491 PetscFunctionBegin;
492 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
493 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
494
495 X = ksp->vec_sol;
496 B = ksp->vec_rhs;
497 R = ksp->work[0];
498 Z = ksp->work[1];
499 U = ksp->work[2];
500 V = ksp->work[3];
501 W = ksp->work[4];
502 UOLD = ksp->work[5];
503 VOLD = ksp->work[6];
504 WOLD = ksp->work[7];
505 WOOLD = ksp->work[8];
506
507 PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
508
509 ksp->its = 0;
510 eigs = ksp->calc_sings;
511 stored_max_it = ksp->max_it;
512 if (eigs) {
513 e = minres->e;
514 d = minres->d;
515 }
516
517 if (!ksp->guess_zero) {
518 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
519 PetscCall(VecAYPX(R, -1.0, B));
520 } else {
521 PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
522 }
523 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
524 PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
525 KSPCheckNorm(ksp, np);
526 PetscCall(VecDot(R, Z, &dp));
527 KSPCheckDot(ksp, dp);
528
529 if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
530 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
531 PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
532 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
533 PetscFunctionReturn(PETSC_SUCCESS);
534 }
535
536 ksp->rnorm = 0.0;
537 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
538 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
539 PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
540 PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
541 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
542
543 dp = PetscAbsScalar(dp);
544 dp = PetscSqrtScalar(dp);
545 beta = dp; /* beta <- sqrt(r'*z) */
546 eta = beta;
547 PetscCall(VecAXPBY(V, 1.0 / beta, 0, R)); /* v <- r / beta */
548 PetscCall(VecAXPBY(U, 1.0 / beta, 0, Z)); /* u <- z / beta */
549
550 i = 0;
551 do {
552 ksp->its = i + 1;
553
554 /* Lanczos */
555
556 PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- A*u */
557 PetscCall(VecDot(U, R, &alpha)); /* alpha <- r'*u */
558 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
559 if (eigs) {
560 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
561 d[i] = alpha;
562 e[i] = beta;
563 }
564
565 if (ksp->its > 1) {
566 Vec T[2];
567 PetscScalar alphas[] = {-alpha, -beta};
568 /* r <- r - alpha v - beta v_old */
569 T[0] = V;
570 T[1] = VOLD;
571 PetscCall(VecMAXPY(R, 2, alphas, T));
572 /* z <- z - alpha u - beta u_old */
573 T[0] = U;
574 T[1] = UOLD;
575 PetscCall(VecMAXPY(Z, 2, alphas, T));
576 } else {
577 PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha v */
578 PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha u */
579 }
580
581 betaold = beta;
582
583 PetscCall(VecDot(R, Z, &dp));
584 KSPCheckDot(ksp, dp);
585 dp = PetscAbsScalar(dp);
586 beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
587
588 /* QR factorisation */
589
590 coold = cold;
591 cold = c;
592 soold = sold;
593 sold = s;
594
595 rho0 = cold * alpha - coold * sold * betaold;
596 rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
597 rho2 = sold * alpha + coold * cold * betaold;
598 rho3 = soold * betaold;
599
600 /* Stop if negative curvature is detected */
601 if (ksp->converged_neg_curve && PetscRealPart(cold * rho0) <= 0.0) {
602 PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1=%g, gbar %g\n", (double)PetscRealPart(cold), -(double)PetscRealPart(rho0)));
603 ksp->reason = KSP_CONVERGED_NEG_CURVE;
604 break;
605 }
606
607 /* Givens rotation */
608
609 c = rho0 / rho1;
610 s = beta / rho1;
611
612 /* Update */
613 /* w_oold <- w_old */
614 /* w_old <- w */
615 KSPMinresSwap3(WOOLD, WOLD, W);
616
617 /* w <- (u - rho2 w_old - rho3 w_oold)/rho1 */
618 PetscCall(VecAXPBY(W, 1.0 / rho1, 0.0, U));
619 if (ksp->its > 1) {
620 Vec T[] = {WOLD, WOOLD};
621 PetscScalar alphas[] = {-rho2 / rho1, -rho3 / rho1};
622 PetscInt nv = (ksp->its == 2 ? 1 : 2);
623
624 PetscCall(VecMAXPY(W, nv, alphas, T));
625 }
626
627 ceta = c * eta;
628 PetscCall(VecAXPY(X, ceta, W)); /* x <- x + c eta w */
629
630 /*
631 when dp is really small we have either convergence or an indefinite operator so compute true
632 residual norm to check for convergence
633 */
634 if (PetscRealPart(dp) < minres->haptol) {
635 PetscCall(PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
636 PetscCall(KSP_MatMult(ksp, Amat, X, VOLD));
637 PetscCall(VecAXPY(VOLD, none, B));
638 PetscCall(VecNorm(VOLD, NORM_2, &np));
639 KSPCheckNorm(ksp, np);
640 } else {
641 /* otherwise compute new residual norm via recurrence relation */
642 np *= PetscAbsScalar(s);
643 }
644
645 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
646 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
647 PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
648 PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
649 if (ksp->reason) break;
650
651 if (PetscRealPart(dp) < minres->haptol) {
652 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
653 PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
654 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
655 break;
656 }
657
658 eta = -s * eta;
659 KSPMinresSwap3(VOLD, V, R);
660 KSPMinresSwap3(UOLD, U, Z);
661 PetscCall(VecScale(V, 1.0 / beta)); /* v <- r / beta */
662 PetscCall(VecScale(U, 1.0 / beta)); /* u <- z / beta */
663
664 i++;
665 } while (i < ksp->max_it);
666 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
667 PetscFunctionReturn(PETSC_SUCCESS);
668 }
669
KSPDestroy_MINRES(KSP ksp)670 static PetscErrorCode KSPDestroy_MINRES(KSP ksp)
671 {
672 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
673
674 PetscFunctionBegin;
675 PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
676 PetscCall(PetscViewerDestroy(&minres->viewer));
677 PetscCall(PetscFree(ksp->data));
678 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", NULL));
679 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", NULL));
680 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", NULL));
681 PetscFunctionReturn(PETSC_SUCCESS);
682 }
683
KSPMINRESSetUseQLP_MINRES(KSP ksp,PetscBool qlp)684 static PetscErrorCode KSPMINRESSetUseQLP_MINRES(KSP ksp, PetscBool qlp)
685 {
686 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
687
688 PetscFunctionBegin;
689 minres->qlp = qlp;
690 PetscFunctionReturn(PETSC_SUCCESS);
691 }
692
KSPMINRESSetRadius_MINRES(KSP ksp,PetscReal radius)693 static PetscErrorCode KSPMINRESSetRadius_MINRES(KSP ksp, PetscReal radius)
694 {
695 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
696
697 PetscFunctionBegin;
698 minres->maxxnorm = radius;
699 PetscFunctionReturn(PETSC_SUCCESS);
700 }
701
KSPMINRESGetUseQLP_MINRES(KSP ksp,PetscBool * qlp)702 static PetscErrorCode KSPMINRESGetUseQLP_MINRES(KSP ksp, PetscBool *qlp)
703 {
704 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
705
706 PetscFunctionBegin;
707 *qlp = minres->qlp;
708 PetscFunctionReturn(PETSC_SUCCESS);
709 }
710
KSPSetFromOptions_MINRES(KSP ksp,PetscOptionItems PetscOptionsObject)711 static PetscErrorCode KSPSetFromOptions_MINRES(KSP ksp, PetscOptionItems PetscOptionsObject)
712 {
713 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
714
715 PetscFunctionBegin;
716 PetscOptionsHeadBegin(PetscOptionsObject, "KSP MINRES options");
717 { /* Allow comparing with the old code (to be removed in a few releases) */
718 PetscBool flg = PETSC_FALSE;
719 PetscCall(PetscOptionsBool("-ksp_minres_old", "Use old implementation (to be removed)", "None", flg, &flg, NULL));
720 if (flg) ksp->ops->solve = KSPSolve_MINRES_OLD;
721 else ksp->ops->solve = KSPSolve_MINRES;
722 }
723 PetscCall(PetscOptionsBool("-ksp_minres_qlp", "Solve with QLP variant", "KSPMINRESSetUseQLP", minres->qlp, &minres->qlp, NULL));
724 PetscCall(PetscOptionsReal("-ksp_minres_radius", "Maximum allowed norm of solution", "KSPMINRESSetRadius", minres->maxxnorm, &minres->maxxnorm, NULL));
725 PetscCall(PetscOptionsReal("-ksp_minres_trancond", "Threshold on condition number to dynamically switch to QLP", "None", minres->TranCond, &minres->TranCond, NULL));
726 PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ksp), PetscOptionsObject->options, PetscOptionsObject->prefix, "-ksp_minres_monitor", &minres->viewer, &minres->viewer_fmt, &minres->monitor));
727 PetscCall(PetscOptionsReal("-ksp_minres_nutol", "Inexactness tolerance", NULL, minres->nutol, &minres->nutol, NULL));
728 PetscOptionsHeadEnd();
729 PetscFunctionReturn(PETSC_SUCCESS);
730 }
731
732 /*@
733 KSPMINRESSetUseQLP - Use the QLP variant of `KSPMINRES`
734
735 Logically Collective
736
737 Input Parameters:
738 + ksp - the iterative context
739 - qlp - a Boolean indicating if the QLP variant should be used
740
741 Level: beginner
742
743 Note:
744 By default, the QLP variant is not used.
745
746 .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESGetUseQLP()`
747 @*/
KSPMINRESSetUseQLP(KSP ksp,PetscBool qlp)748 PetscErrorCode KSPMINRESSetUseQLP(KSP ksp, PetscBool qlp)
749 {
750 PetscFunctionBegin;
751 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
752 PetscValidLogicalCollectiveBool(ksp, qlp, 2);
753 PetscTryMethod(ksp, "KSPMINRESSetUseQLP_C", (KSP, PetscBool), (ksp, qlp));
754 PetscFunctionReturn(PETSC_SUCCESS);
755 }
756
757 /*@
758 KSPMINRESSetRadius - Set the maximum solution norm allowed for use with trust region methods
759
760 Logically Collective
761
762 Input Parameters:
763 + ksp - the iterative context
764 - radius - the value
765
766 Level: beginner
767
768 Options Database Key:
769 . -ksp_minres_radius <real> - maximum allowed solution norm
770
771 Developer Note:
772 Perhaps the KSPXXXSetRadius() should be unified
773
774 .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
775 @*/
KSPMINRESSetRadius(KSP ksp,PetscReal radius)776 PetscErrorCode KSPMINRESSetRadius(KSP ksp, PetscReal radius)
777 {
778 PetscFunctionBegin;
779 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
780 PetscValidLogicalCollectiveReal(ksp, radius, 2);
781 PetscTryMethod(ksp, "KSPMINRESSetRadius_C", (KSP, PetscReal), (ksp, radius));
782 PetscFunctionReturn(PETSC_SUCCESS);
783 }
784
785 /*@
786 KSPMINRESGetUseQLP - Get the flag that indicates if the QLP variant is being used
787
788 Logically Collective
789
790 Input Parameter:
791 . ksp - the iterative context
792
793 Output Parameter:
794 . qlp - a Boolean indicating if the QLP variant is used
795
796 Level: beginner
797
798 .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
799 @*/
KSPMINRESGetUseQLP(KSP ksp,PetscBool * qlp)800 PetscErrorCode KSPMINRESGetUseQLP(KSP ksp, PetscBool *qlp)
801 {
802 PetscFunctionBegin;
803 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
804 PetscAssertPointer(qlp, 2);
805 PetscUseMethod(ksp, "KSPMINRESGetUseQLP_C", (KSP, PetscBool *), (ksp, qlp));
806 PetscFunctionReturn(PETSC_SUCCESS);
807 }
808
809 /*MC
810 KSPMINRES - This code implements the MINRES (Minimum Residual) method and its QLP variant {cite}`paige.saunders:solution`, {cite}`choi2011minres`,
811 {cite}`liu2022newton` for solving linear systems using `KSP`.
812
813 Options Database Keys:
814 + -ksp_minres_qlp <bool> - activates QLP code
815 . -ksp_minres_radius <real> - maximum allowed solution norm
816 . -ksp_minres_trancond <real> - threshold on condition number to dynamically switch to QLP iterations when QLP has been activated
817 . -ksp_minres_monitor - monitors convergence quantities
818 - -ksp_minres_nutol <real> - inexactness tolerance (see https://arxiv.org/pdf/2208.07095.pdf)
819
820 Level: beginner
821
822 Notes:
823 The matrix (operator) and the preconditioner must be symmetric and the preconditioner must also be positive definite for this method.
824
825 `KSPMINRES` is often the best Krylov method for symmetric indefinite matrices.
826
827 Supports only left preconditioning.
828
829 Contributed by:
830 Original MINRES code - Robert Scheichl: maprs@maths.bath.ac.uk
831 QLP variant adapted from: https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
832
833 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`, `KSPMINRESGetUseQLP()`, `KSPMINRESSetUseQLP()`, `KSPMINRESSetRadius()`
834 M*/
KSPCreate_MINRES(KSP ksp)835 PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
836 {
837 KSP_MINRES *minres;
838
839 PetscFunctionBegin;
840 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
841 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
842 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
843 PetscCall(PetscNew(&minres));
844
845 /* this parameter is arbitrary and belongs to the old implementation; but e-50 didn't work for __float128 in one example */
846 #if defined(PETSC_USE_REAL___FLOAT128)
847 minres->haptol = 1.e-100;
848 #elif defined(PETSC_USE_REAL_SINGLE)
849 minres->haptol = 1.e-25;
850 #else
851 minres->haptol = 1.e-50;
852 #endif
853 /* those are set as 1.e7 in the MATLAB code -> use 1.0/sqrt(eps) to support single precision */
854 minres->maxxnorm = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
855 minres->TranCond = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
856
857 ksp->data = (void *)minres;
858
859 ksp->ops->setup = KSPSetUp_MINRES;
860 ksp->ops->solve = KSPSolve_MINRES;
861 ksp->ops->destroy = KSPDestroy_MINRES;
862 ksp->ops->setfromoptions = KSPSetFromOptions_MINRES;
863 ksp->ops->buildsolution = KSPBuildSolutionDefault;
864 ksp->ops->buildresidual = KSPBuildResidualDefault;
865
866 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", KSPMINRESSetRadius_MINRES));
867 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", KSPMINRESSetUseQLP_MINRES));
868 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", KSPMINRESGetUseQLP_MINRES));
869 PetscFunctionReturn(PETSC_SUCCESS);
870 }
871
KSPComputeEigenvalues_MINRES(KSP ksp,PetscInt nmax,PetscReal * r,PetscReal * c,PetscInt * neig)872 PetscErrorCode KSPComputeEigenvalues_MINRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
873 {
874 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
875 PetscScalar *d, *e;
876 PetscReal *ee;
877 PetscInt n = ksp->its;
878 PetscBLASInt bn, lierr = 0, ldz = 1;
879
880 PetscFunctionBegin;
881 PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
882 *neig = n;
883
884 PetscCall(PetscArrayzero(c, nmax));
885 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
886 d = minres->d;
887 e = minres->e;
888 ee = minres->ee;
889
890 /* copy tridiagonal matrix to work space */
891 for (PetscInt j = 0; j < n; j++) {
892 r[j] = PetscRealPart(d[j]);
893 ee[j] = PetscRealPart(e[j + 1]);
894 }
895
896 PetscCall(PetscBLASIntCast(n, &bn));
897 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
898 PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, ee, NULL, &ldz, NULL, &lierr));
899 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
900 PetscCall(PetscFPTrapPop());
901 PetscCall(PetscSortReal(n, r));
902 PetscFunctionReturn(PETSC_SUCCESS);
903 }
904
KSPComputeExtremeSingularValues_MINRES(KSP ksp,PetscReal * emax,PetscReal * emin)905 PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP ksp, PetscReal *emax, PetscReal *emin)
906 {
907 KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
908 PetscScalar *d, *e;
909 PetscReal *dd, *ee;
910 PetscInt n = ksp->its;
911 PetscBLASInt bn, lierr = 0, ldz = 1;
912
913 PetscFunctionBegin;
914 if (!n) {
915 *emax = *emin = 1.0;
916 PetscFunctionReturn(PETSC_SUCCESS);
917 }
918 d = minres->d;
919 e = minres->e;
920 dd = minres->dd;
921 ee = minres->ee;
922
923 /* copy tridiagonal matrix to work space */
924 for (PetscInt j = 0; j < n; j++) {
925 dd[j] = PetscRealPart(d[j]);
926 ee[j] = PetscRealPart(e[j + 1]);
927 }
928
929 PetscCall(PetscBLASIntCast(n, &bn));
930 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
931 PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, ee, NULL, &ldz, NULL, &lierr));
932 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
933 PetscCall(PetscFPTrapPop());
934 for (PetscInt j = 0; j < n; j++) dd[j] = PetscAbsReal(dd[j]);
935 PetscCall(PetscSortReal(n, dd));
936 *emin = dd[0];
937 *emax = dd[n - 1];
938 PetscFunctionReturn(PETSC_SUCCESS);
939 }
940