Actual source code: minres.c
1: #include <petsc/private/kspimpl.h>
2: #include <petscblaslapack.h>
3: PETSC_INTERN PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP, PetscReal *, PetscReal *);
4: PETSC_INTERN PetscErrorCode KSPComputeEigenvalues_MINRES(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
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";
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;
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));
44: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_MINRES;
45: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_MINRES;
46: }
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
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)
59: static inline PetscReal Norm3(PetscReal a, PetscReal b, PetscReal c)
60: {
61: return PetscSqrtReal(PetscSqr(a) + PetscSqr(b) + PetscSqr(c));
62: }
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;
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;
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: }
90: /*
91: Code adapted from https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
92: CSP11/Algorithms/MINRESQLP/minresQLP.m
93: */
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, xnorml = 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;
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);
120: eigs = ksp->calc_sings;
121: stored_max_it = ksp->max_it;
122: if (eigs) {
123: e = minres->e;
124: d = minres->d;
125: }
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));
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);
163: rnorm = beta1;
164: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
165: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
166: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
167: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
168: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
170: relres = rnorm / beta1;
171: betan = beta1;
172: phi0 = beta1;
173: phi = beta1;
174: betan = beta1;
175: beta = 0.0;
176: do {
177: /* Lanczos */
178: ksp->its++;
179: betal = beta;
180: beta = betan;
181: PetscCall(VecAXPBY(V, 1.0 / beta, 0.0, R3));
182: PetscCall(KSP_MatMult(ksp, Amat, V, R3));
183: if (ksp->its > 1) PetscCall(VecAXPY(R3, -beta / betal, R1));
184: PetscCall(VecDotRealPart(R3, V, &alpha));
185: PetscCall(VecAXPY(R3, -alpha / beta, R2));
186: KSPMinresSwap3(R1, R2, R3);
187: if (eigs) {
188: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
189: d[ksp->its - 1] = alpha;
190: e[ksp->its - 1] = beta;
191: }
193: PetscCall(KSP_PCApply(ksp, R2, R3));
194: PetscCall(VecDotRealPart(R3, R2, &betan));
195: KSPCheckDot(ksp, betan);
196: if (betan < 0.0) {
197: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite preconditioner %g", (double)betan);
198: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
201: betan = PetscSqrtReal(betan);
203: pnorm = Norm3(betal, alpha, betan);
205: // Apply previous left rotation Q_{k-1}
206: dbar = dltan;
207: epln = eplnn;
208: dlta = cs * dbar + sn * alpha;
209: gbar = sn * dbar - cs * alpha;
210: eplnn = sn * betan;
211: dltan = -cs * betan;
212: dlta_QLP = dlta;
214: // Stop if negative curvature is detected and return residual
215: // This is very experimental and maybe changed in the future
216: // based on https://arxiv.org/pdf/2208.07095.pdf
217: if (ksp->converged_neg_curve) {
218: if (cs * gbar >= 0.0) {
219: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1 %g, gbar %g\n", (double)cs, (double)gbar));
220: ksp->reason = KSP_CONVERGED_NEG_CURVE;
221: PetscCall(VecCopy(RN, X));
222: break;
223: } else {
224: PetscCall(VecAXPBY(RN, -phi * cs, PetscSqr(sn), V));
225: }
226: }
228: // Compute the current left plane rotation Q_k
229: gamal3 = gamal2;
230: gamal2 = gamal;
231: gamal = gama;
232: SymOrtho(gbar, betan, &cs, &sn, &gama);
234: // Inexactness condition from https://arxiv.org/pdf/2208.07095.pdf
235: rootl = Norm3(gbar, dltan, zero);
236: phir = PetscSqr(phi0 / phi);
237: if (ksp->its > 2 && minres->nutol > 0.0) {
238: PetscReal tmp;
240: phir = PetscSqrtReal(phir - 1.0);
241: tmp = rootl / phir;
242: PetscCall(PetscInfo(ksp, "it = %" PetscInt_FMT ": inexact check %g (%g / %g)\n", ksp->its - 2, (double)tmp, (double)rootl, (double)phir));
243: if (tmp < minres->nutol) {
244: ksp->its--;
245: ksp->reason = KSP_CONVERGED_RTOL;
246: break;
247: }
248: }
250: gama_tmp = gama;
251: taul2 = taul;
252: taul = tau;
253: tau = cs * phi;
254: Axnorm = Norm3(Axnorm, tau, zero);
255: phi = sn * phi;
257: //Apply the previous right plane rotation P{k-2,k}
258: if (ksp->its > 2) {
259: veplnl2 = veplnl;
260: etal2 = etal;
261: etal = eta;
262: dlta_tmp = sr2 * vepln - cr2 * dlta;
263: veplnl = cr2 * vepln + sr2 * dlta;
264: dlta = dlta_tmp;
265: eta = sr2 * gama;
266: gama = -cr2 * gama;
267: }
269: // Compute the current right plane rotation P{k-1,k}, P_12, P_23,...
270: if (ksp->its > 1) {
271: SymOrtho(gamal, dlta, &cr1, &sr1, &gamal);
272: vepln = sr1 * gama;
273: gama = -cr1 * gama;
274: }
276: // Update xnorm
277: xnorml = 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);
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);
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};
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: }
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);
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;
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;
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: }
416: if (flag == 2 || flag == 4 || flag == 6 || flag == 7) {
417: Acond = Acondl;
418: rnorm = rnorml;
419: relres = relresl;
420: }
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 %5d %2d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", (int)ksp->its - 1, (int)flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorml));
426: PetscCall(PetscViewerPopFormat(minres->viewer));
427: }
429: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
430: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
431: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
432: PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
433: if (!ksp->reason) {
434: switch (flag) {
435: case 1:
436: case 2:
437: case 5: /* XXX */
438: ksp->reason = KSP_CONVERGED_RTOL;
439: break;
440: case 3:
441: case 4:
442: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
443: break;
444: case 6:
445: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
446: break;
447: default:
448: break;
449: }
450: }
451: if (ksp->reason) break;
452: } while (ksp->its < ksp->max_it);
454: if (minres->monitor && flag != 2 && flag != 4 && flag != 6 && flag != 7) {
455: PetscCall(VecNorm(X, NORM_2, &xnorm));
456: PetscCall(KSP_MatMult(ksp, Amat, X, R1));
457: PetscCall(VecAYPX(R1, -1.0, B));
458: PetscCall(VecNorm(R1, NORM_2, &rnorml));
459: PetscCall(KSP_MatMult(ksp, Amat, R1, R2));
460: PetscCall(VecNorm(R2, NORM_2, &Arnorml));
461: relresl = rnorml / (Anorm * xnorm + beta1);
462: relAresl = rnorml > realmin ? Arnorml / (Anorm * rnorml) : 0.0;
463: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
464: PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5d %2d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", (int)ksp->its, (int)flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorml));
465: PetscCall(PetscViewerPopFormat(minres->viewer));
466: }
467: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /* This was the original implementation provided by R. Scheichl */
472: static PetscErrorCode KSPSolve_MINRES_OLD(KSP ksp)
473: {
474: PetscInt i;
475: PetscScalar alpha, beta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
476: PetscScalar rho0, rho1, rho2, rho3, dp = 0.0;
477: const PetscScalar none = -1.0;
478: PetscReal np;
479: Vec X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
480: Mat Amat;
481: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
482: PetscBool diagonalscale;
483: PetscInt stored_max_it, eigs;
484: PetscScalar *e = NULL, *d = NULL;
486: PetscFunctionBegin;
487: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
488: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
490: X = ksp->vec_sol;
491: B = ksp->vec_rhs;
492: R = ksp->work[0];
493: Z = ksp->work[1];
494: U = ksp->work[2];
495: V = ksp->work[3];
496: W = ksp->work[4];
497: UOLD = ksp->work[5];
498: VOLD = ksp->work[6];
499: WOLD = ksp->work[7];
500: WOOLD = ksp->work[8];
502: PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
504: ksp->its = 0;
505: eigs = ksp->calc_sings;
506: stored_max_it = ksp->max_it;
507: if (eigs) {
508: e = minres->e;
509: d = minres->d;
510: }
512: if (!ksp->guess_zero) {
513: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
514: PetscCall(VecAYPX(R, -1.0, B));
515: } else {
516: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
517: }
518: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
519: PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
520: KSPCheckNorm(ksp, np);
521: PetscCall(VecDot(R, Z, &dp));
522: KSPCheckDot(ksp, dp);
524: if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
525: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
526: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
527: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: ksp->rnorm = 0.0;
532: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
533: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
534: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
535: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
536: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
538: dp = PetscAbsScalar(dp);
539: dp = PetscSqrtScalar(dp);
540: beta = dp; /* beta <- sqrt(r'*z) */
541: eta = beta;
542: PetscCall(VecAXPBY(V, 1.0 / beta, 0, R)); /* v <- r / beta */
543: PetscCall(VecAXPBY(U, 1.0 / beta, 0, Z)); /* u <- z / beta */
545: i = 0;
546: do {
547: ksp->its = i + 1;
549: /* Lanczos */
551: PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- A*u */
552: PetscCall(VecDot(U, R, &alpha)); /* alpha <- r'*u */
553: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
554: if (eigs) {
555: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
556: d[i] = alpha;
557: e[i] = beta;
558: }
560: if (ksp->its > 1) {
561: Vec T[2];
562: PetscScalar alphas[] = {-alpha, -beta};
563: /* r <- r - alpha v - beta v_old */
564: T[0] = V;
565: T[1] = VOLD;
566: PetscCall(VecMAXPY(R, 2, alphas, T));
567: /* z <- z - alpha u - beta u_old */
568: T[0] = U;
569: T[1] = UOLD;
570: PetscCall(VecMAXPY(Z, 2, alphas, T));
571: } else {
572: PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha v */
573: PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha u */
574: }
576: betaold = beta;
578: PetscCall(VecDot(R, Z, &dp));
579: KSPCheckDot(ksp, dp);
580: dp = PetscAbsScalar(dp);
581: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
583: /* QR factorisation */
585: coold = cold;
586: cold = c;
587: soold = sold;
588: sold = s;
590: rho0 = cold * alpha - coold * sold * betaold;
591: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
592: rho2 = sold * alpha + coold * cold * betaold;
593: rho3 = soold * betaold;
595: /* Stop if negative curvature is detected */
596: if (ksp->converged_neg_curve && PetscRealPart(cold * rho0) <= 0.0) {
597: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1=%g, gbar %g\n", (double)PetscRealPart(cold), -(double)PetscRealPart(rho0)));
598: ksp->reason = KSP_CONVERGED_NEG_CURVE;
599: break;
600: }
602: /* Givens rotation */
604: c = rho0 / rho1;
605: s = beta / rho1;
607: /* Update */
608: /* w_oold <- w_old */
609: /* w_old <- w */
610: KSPMinresSwap3(WOOLD, WOLD, W);
612: /* w <- (u - rho2 w_old - rho3 w_oold)/rho1 */
613: PetscCall(VecAXPBY(W, 1.0 / rho1, 0.0, U));
614: if (ksp->its > 1) {
615: Vec T[] = {WOLD, WOOLD};
616: PetscScalar alphas[] = {-rho2 / rho1, -rho3 / rho1};
617: PetscInt nv = (ksp->its == 2 ? 1 : 2);
619: PetscCall(VecMAXPY(W, nv, alphas, T));
620: }
622: ceta = c * eta;
623: PetscCall(VecAXPY(X, ceta, W)); /* x <- x + c eta w */
625: /*
626: when dp is really small we have either convergence or an indefinite operator so compute true
627: residual norm to check for convergence
628: */
629: if (PetscRealPart(dp) < minres->haptol) {
630: PetscCall(PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
631: PetscCall(KSP_MatMult(ksp, Amat, X, VOLD));
632: PetscCall(VecAXPY(VOLD, none, B));
633: PetscCall(VecNorm(VOLD, NORM_2, &np));
634: KSPCheckNorm(ksp, np);
635: } else {
636: /* otherwise compute new residual norm via recurrence relation */
637: np *= PetscAbsScalar(s);
638: }
640: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
641: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
642: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
643: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
644: if (ksp->reason) break;
646: if (PetscRealPart(dp) < minres->haptol) {
647: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
648: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
649: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
650: break;
651: }
653: eta = -s * eta;
654: KSPMinresSwap3(VOLD, V, R);
655: KSPMinresSwap3(UOLD, U, Z);
656: PetscCall(VecScale(V, 1.0 / beta)); /* v <- r / beta */
657: PetscCall(VecScale(U, 1.0 / beta)); /* u <- z / beta */
659: i++;
660: } while (i < ksp->max_it);
661: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: static PetscErrorCode KSPDestroy_MINRES(KSP ksp)
666: {
667: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
669: PetscFunctionBegin;
670: PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
671: PetscCall(PetscOptionsRestoreViewer(&minres->viewer));
672: PetscCall(PetscFree(ksp->data));
673: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", NULL));
674: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", NULL));
675: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", NULL));
676: PetscFunctionReturn(PETSC_SUCCESS);
677: }
679: static PetscErrorCode KSPMINRESSetUseQLP_MINRES(KSP ksp, PetscBool qlp)
680: {
681: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
683: PetscFunctionBegin;
684: minres->qlp = qlp;
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: static PetscErrorCode KSPMINRESSetRadius_MINRES(KSP ksp, PetscReal radius)
689: {
690: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
692: PetscFunctionBegin;
693: minres->maxxnorm = radius;
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: static PetscErrorCode KSPMINRESGetUseQLP_MINRES(KSP ksp, PetscBool *qlp)
698: {
699: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
701: PetscFunctionBegin;
702: *qlp = minres->qlp;
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: static PetscErrorCode KSPSetFromOptions_MINRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
707: {
708: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
710: PetscFunctionBegin;
711: PetscOptionsHeadBegin(PetscOptionsObject, "KSP MINRES options");
712: { /* Allow comparing with the old code (to be removed in a few releases) */
713: PetscBool flg = PETSC_FALSE;
714: PetscCall(PetscOptionsBool("-ksp_minres_old", "Use old implementation (to be removed)", "None", flg, &flg, NULL));
715: if (flg) ksp->ops->solve = KSPSolve_MINRES_OLD;
716: else ksp->ops->solve = KSPSolve_MINRES;
717: }
718: PetscCall(PetscOptionsBool("-ksp_minres_qlp", "Solve with QLP variant", "KSPMINRESSetUseQLP", minres->qlp, &minres->qlp, NULL));
719: PetscCall(PetscOptionsReal("-ksp_minres_radius", "Maximum allowed norm of solution", "KSPMINRESSetRadius", minres->maxxnorm, &minres->maxxnorm, NULL));
720: PetscCall(PetscOptionsReal("-ksp_minres_trancond", "Threshold on condition number to dynamically switch to QLP", "None", minres->TranCond, &minres->TranCond, NULL));
721: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp), PetscOptionsObject->options, PetscOptionsObject->prefix, "-ksp_minres_monitor", &minres->viewer, &minres->viewer_fmt, &minres->monitor));
722: PetscCall(PetscOptionsReal("-ksp_minres_nutol", "Inexactness tolerance", NULL, minres->nutol, &minres->nutol, NULL));
723: PetscOptionsHeadEnd();
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: /*@
728: KSPMINRESSetUseQLP - Use the QLP variant of `KSPMINRES`
730: Logically Collective
732: Input Parameters:
733: + ksp - the iterative context
734: - qlp - a Boolean indicating if the QLP variant should be used
736: Level: beginner
738: Note:
739: By default, the QLP variant is not used.
741: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESGetUseQLP()`
742: @*/
743: PetscErrorCode KSPMINRESSetUseQLP(KSP ksp, PetscBool qlp)
744: {
745: PetscFunctionBegin;
748: PetscTryMethod(ksp, "KSPMINRESSetUseQLP_C", (KSP, PetscBool), (ksp, qlp));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: /*@
753: KSPMINRESSetRadius - Set the maximum solution norm allowed for use with trust region methods
755: Logically Collective
757: Input Parameters:
758: + ksp - the iterative context
759: - radius - the value
761: Level: beginner
763: Options Database Key:
764: . -ksp_minres_radius <real> - maximum allowed solution norm
766: Developer Note:
767: Perhaps the KSPXXXSetRadius() should be unified
769: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
770: @*/
771: PetscErrorCode KSPMINRESSetRadius(KSP ksp, PetscReal radius)
772: {
773: PetscFunctionBegin;
776: PetscTryMethod(ksp, "KSPMINRESSetRadius_C", (KSP, PetscReal), (ksp, radius));
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: /*@
781: KSPMINRESGetUseQLP - Get the flag that indicates if the QLP variant is being used
783: Logically Collective
785: Input Parameter:
786: . ksp - the iterative context
788: Output Parameter:
789: . qlp - a Boolean indicating if the QLP variant is used
791: Level: beginner
793: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
794: @*/
795: PetscErrorCode KSPMINRESGetUseQLP(KSP ksp, PetscBool *qlp)
796: {
797: PetscFunctionBegin;
799: PetscAssertPointer(qlp, 2);
800: PetscUseMethod(ksp, "KSPMINRESGetUseQLP_C", (KSP, PetscBool *), (ksp, qlp));
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*MC
805: KSPMINRES - This code implements the MINRES (Minimum Residual) method and its QLP variant {cite}`paige.saunders:solution`, {cite}`choi2011minres`,
806: {cite}`liu2022newton`.
808: Options Database Keys:
809: + -ksp_minres_qlp <bool> - activates QLP code
810: . -ksp_minres_radius <real> - maximum allowed solution norm
811: . -ksp_minres_trancond <real> - threshold on condition number to dynamically switch to QLP iterations when QLP has been activated
812: . -ksp_minres_monitor - monitors convergence quantities
813: - -ksp_minres_nutol <real> - inexactness tolerance (see https://arxiv.org/pdf/2208.07095.pdf)
815: Level: beginner
817: Notes:
818: The operator and the preconditioner must be symmetric and the preconditioner must be positive definite for this method.
820: Supports only left preconditioning.
822: Contributed by:
823: Original MINRES code - Robert Scheichl: maprs@maths.bath.ac.uk
824: QLP variant adapted from: https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
826: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`, `KSPMINRESGetUseQLP()`, `KSPMINRESSetUseQLP()`, `KSPMINRESSetRadius()`
827: M*/
828: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
829: {
830: KSP_MINRES *minres;
832: PetscFunctionBegin;
833: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
834: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
835: PetscCall(PetscNew(&minres));
837: /* this parameter is arbitrary and belongs to the old implementation; but e-50 didn't work for __float128 in one example */
838: #if defined(PETSC_USE_REAL___FLOAT128)
839: minres->haptol = 1.e-100;
840: #elif defined(PETSC_USE_REAL_SINGLE)
841: minres->haptol = 1.e-25;
842: #else
843: minres->haptol = 1.e-50;
844: #endif
845: /* those are set as 1.e7 in the MATLAB code -> use 1.0/sqrt(eps) to support single precision */
846: minres->maxxnorm = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
847: minres->TranCond = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
849: ksp->data = (void *)minres;
851: ksp->ops->setup = KSPSetUp_MINRES;
852: ksp->ops->solve = KSPSolve_MINRES;
853: ksp->ops->destroy = KSPDestroy_MINRES;
854: ksp->ops->setfromoptions = KSPSetFromOptions_MINRES;
855: ksp->ops->buildsolution = KSPBuildSolutionDefault;
856: ksp->ops->buildresidual = KSPBuildResidualDefault;
858: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", KSPMINRESSetRadius_MINRES));
859: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", KSPMINRESSetUseQLP_MINRES));
860: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", KSPMINRESGetUseQLP_MINRES));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: PetscErrorCode KSPComputeEigenvalues_MINRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
865: {
866: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
867: PetscScalar *d, *e;
868: PetscReal *ee;
869: PetscInt n = ksp->its;
870: PetscBLASInt bn, lierr = 0, ldz = 1;
872: PetscFunctionBegin;
873: PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
874: *neig = n;
876: PetscCall(PetscArrayzero(c, nmax));
877: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
878: d = minres->d;
879: e = minres->e;
880: ee = minres->ee;
882: /* copy tridiagonal matrix to work space */
883: for (PetscInt j = 0; j < n; j++) {
884: r[j] = PetscRealPart(d[j]);
885: ee[j] = PetscRealPart(e[j + 1]);
886: }
888: PetscCall(PetscBLASIntCast(n, &bn));
889: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
890: PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, ee, NULL, &ldz, NULL, &lierr));
891: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
892: PetscCall(PetscFPTrapPop());
893: PetscCall(PetscSortReal(n, r));
894: PetscFunctionReturn(PETSC_SUCCESS);
895: }
897: PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP ksp, PetscReal *emax, PetscReal *emin)
898: {
899: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
900: PetscScalar *d, *e;
901: PetscReal *dd, *ee;
902: PetscInt n = ksp->its;
903: PetscBLASInt bn, lierr = 0, ldz = 1;
905: PetscFunctionBegin;
906: if (!n) {
907: *emax = *emin = 1.0;
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
910: d = minres->d;
911: e = minres->e;
912: dd = minres->dd;
913: ee = minres->ee;
915: /* copy tridiagonal matrix to work space */
916: for (PetscInt j = 0; j < n; j++) {
917: dd[j] = PetscRealPart(d[j]);
918: ee[j] = PetscRealPart(e[j + 1]);
919: }
921: PetscCall(PetscBLASIntCast(n, &bn));
922: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
923: PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, ee, NULL, &ldz, NULL, &lierr));
924: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
925: PetscCall(PetscFPTrapPop());
926: for (PetscInt j = 0; j < n; j++) dd[j] = PetscAbsReal(dd[j]);
927: PetscCall(PetscSortReal(n, dd));
928: *emin = dd[0];
929: *emax = dd[n - 1];
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }