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;
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_PRECONDITIONED) ksp->rnorm = rnorm;
165: else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
166: PetscCall(VecNorm(R2, NORM_2, &ksp->rnorm));
167: }
168: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
169: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
170: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
171: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
173: relres = rnorm / beta1;
174: betan = beta1;
175: phi0 = beta1;
176: phi = beta1;
177: betan = beta1;
178: beta = 0.0;
179: do {
180: /* Lanczos */
181: ksp->its++;
182: betal = beta;
183: beta = betan;
184: PetscCall(VecAXPBY(V, 1.0 / beta, 0.0, R3));
185: PetscCall(KSP_MatMult(ksp, Amat, V, R3));
186: if (ksp->its > 1) PetscCall(VecAXPY(R3, -beta / betal, R1));
187: PetscCall(VecDotRealPart(R3, V, &alpha));
188: PetscCall(VecAXPY(R3, -alpha / beta, R2));
189: KSPMinresSwap3(R1, R2, R3);
190: if (eigs) {
191: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
192: d[ksp->its - 1] = alpha;
193: e[ksp->its - 1] = beta;
194: }
196: PetscCall(KSP_PCApply(ksp, R2, R3));
197: PetscCall(VecDotRealPart(R3, R2, &betan));
198: KSPCheckDot(ksp, betan);
199: if (betan < 0.0) {
200: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite preconditioner %g", (double)betan);
201: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
204: betan = PetscSqrtReal(betan);
206: pnorm = Norm3(betal, alpha, betan);
208: // Apply previous left rotation Q_{k-1}
209: dbar = dltan;
210: epln = eplnn;
211: dlta = cs * dbar + sn * alpha;
212: gbar = sn * dbar - cs * alpha;
213: eplnn = sn * betan;
214: dltan = -cs * betan;
215: dlta_QLP = dlta;
217: // Stop if negative curvature is detected and return residual
218: // This is very experimental and maybe changed in the future
219: // based on https://arxiv.org/pdf/2208.07095.pdf
220: if (ksp->converged_neg_curve) {
221: if (cs * gbar >= 0.0) {
222: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1 %g, gbar %g\n", (double)cs, (double)gbar));
223: ksp->reason = KSP_CONVERGED_NEG_CURVE;
224: PetscCall(VecCopy(RN, X));
225: break;
226: } else {
227: PetscCall(VecAXPBY(RN, -phi * cs, PetscSqr(sn), V));
228: }
229: }
231: // Compute the current left plane rotation Q_k
232: gamal3 = gamal2;
233: gamal2 = gamal;
234: gamal = gama;
235: SymOrtho(gbar, betan, &cs, &sn, &gama);
237: // Inexactness condition from https://arxiv.org/pdf/2208.07095.pdf
238: rootl = Norm3(gbar, dltan, zero);
239: phir = PetscSqr(phi0 / phi);
240: if (ksp->its > 2 && minres->nutol > 0.0) {
241: PetscReal tmp;
243: phir = PetscSqrtReal(phir - 1.0);
244: tmp = rootl / phir;
245: PetscCall(PetscInfo(ksp, "it = %" PetscInt_FMT ": inexact check %g (%g / %g)\n", ksp->its - 2, (double)tmp, (double)rootl, (double)phir));
246: if (tmp < minres->nutol) {
247: ksp->its--;
248: ksp->reason = KSP_CONVERGED_RTOL;
249: break;
250: }
251: }
253: gama_tmp = gama;
254: taul2 = taul;
255: taul = tau;
256: tau = cs * phi;
257: Axnorm = Norm3(Axnorm, tau, zero);
258: phi = sn * phi;
260: //Apply the previous right plane rotation P{k-2,k}
261: if (ksp->its > 2) {
262: veplnl2 = veplnl;
263: etal2 = etal;
264: etal = eta;
265: dlta_tmp = sr2 * vepln - cr2 * dlta;
266: veplnl = cr2 * vepln + sr2 * dlta;
267: dlta = dlta_tmp;
268: eta = sr2 * gama;
269: gama = -cr2 * gama;
270: }
272: // Compute the current right plane rotation P{k-1,k}, P_12, P_23,...
273: if (ksp->its > 1) {
274: SymOrtho(gamal, dlta, &cr1, &sr1, &gamal);
275: vepln = sr1 * gama;
276: gama = -cr1 * gama;
277: }
279: // Update xnorm
280: ul4 = ul3;
281: ul3 = ul2;
282: if (ksp->its > 2) ul2 = (taul2 - etal2 * ul4 - veplnl2 * ul3) / gamal2;
283: if (ksp->its > 1) ul = (taul - etal * ul3 - veplnl * ul2) / gamal;
284: xnorm_tmp = Norm3(xl2norm, ul2, ul);
285: if (PetscAbsReal(gama) > realmin && xnorm_tmp < minres->maxxnorm) {
286: u = (tau - eta * ul2 - vepln * ul) / gama;
287: if (Norm3(xnorm_tmp, u, zero) > minres->maxxnorm) {
288: u = 0;
289: flag = 6;
290: }
291: } else {
292: u = 0;
293: flag = 9;
294: }
295: xl2norm = Norm3(xl2norm, ul2, zero);
296: xnorm = Norm3(xl2norm, ul, u);
298: // Update w. Update x except if it will become too big
299: //if (Acond < minres->TranCond && flag != flag0 && QLPiter == 0) { // I believe they have a typo in the MATLAB code
300: if ((Acond < minres->TranCond || !minres->qlp) && flag == flag0 && QLPiter == 0) { // MINRES
301: KSPMinresSwap3(WL2, WL, W);
302: PetscCall(VecAXPBY(W, 1.0 / gama_tmp, 0.0, V));
303: if (ksp->its > 1) {
304: Vec T[] = {WL, WL2};
305: PetscScalar alphas[] = {-dlta_QLP / gama_tmp, -epln / gama_tmp};
306: PetscInt nv = (ksp->its == 2 ? 1 : 2);
308: PetscCall(VecMAXPY(W, nv, alphas, T));
309: }
310: if (xnorm < minres->maxxnorm) {
311: PetscCall(VecAXPY(X, tau, W));
312: } else {
313: flag = 6;
314: }
315: } else if (minres->qlp) { //MINRES-QLP updates
316: QLPiter = QLPiter + 1;
317: if (QLPiter == 1) {
318: // xl2 = x - wl*ul_QLP - w*u_QLP;
319: PetscScalar maxpys[] = {1.0, -ul_QLP, -u_QLP};
320: Vec maxpyv[] = {X, WL, W};
322: PetscCall(VecSet(XL2, 0.0));
323: // construct w_{k-3}, w_{k-2}, w_{k-1}
324: if (ksp->its > 1) {
325: if (ksp->its > 3) { // w_{k-3}
326: //wl2 = gamal3*wl2 + veplnl2*wl + etal*w;
327: PetscCall(VecAXPBYPCZ(WL2, veplnl2, etal, gamal3, WL, W));
328: }
329: if (ksp->its > 2) { // w_{k-2}
330: //wl = gamal_QLP*wl + vepln_QLP*w;
331: PetscCall(VecAXPBY(WL, vepln_QLP, gamal_QLP, W));
332: }
333: // w = gama_QLP*w;
334: PetscCall(VecScale(W, gama_QLP));
335: // xl2 = x - wl*ul_QLP - w*u_QLP;
336: PetscCall(VecMAXPY(XL2, 3, maxpys, maxpyv));
337: }
338: }
339: if (ksp->its == 1) {
340: //wl2 = wl; wl = v*sr1; w = -v*cr1;
341: PetscCall(VecCopy(WL, WL2));
342: PetscCall(VecAXPBY(WL, sr1, 0, V));
343: PetscCall(VecAXPBY(W, -cr1, 0, V));
344: } else if (ksp->its == 2) {
345: //wl2 = wl;
346: //wl = w*cr1 + v*sr1;
347: //w = w*sr1 - v*cr1;
348: PetscCall(VecCopy(WL, WL2));
349: PetscCall(VecAXPBYPCZ(WL, cr1, sr1, 0.0, W, V));
350: PetscCall(VecAXPBY(W, -cr1, sr1, V));
351: } else {
352: //wl2 = wl; wl = w; w = wl2*sr2 - v*cr2;
353: //wl2 = wl2*cr2 + v*sr2; v = wl *cr1 + w*sr1;
354: //w = wl *sr1 - w*cr1; wl = v;
355: PetscCall(VecCopy(WL, WL2));
356: PetscCall(VecCopy(W, WL));
357: PetscCall(VecAXPBYPCZ(W, sr2, -cr2, 0.0, WL2, V));
358: PetscCall(VecAXPBY(WL2, sr2, cr2, V));
359: PetscCall(VecAXPBYPCZ(V, cr1, sr1, 0.0, WL, W));
360: PetscCall(VecAXPBY(W, sr1, -cr1, WL));
361: PetscCall(VecCopy(V, WL));
362: }
364: //xl2 = xl2 + wl2*ul2;
365: PetscCall(VecAXPY(XL2, ul2, WL2));
366: //x = xl2 + wl *ul + w*u;
367: PetscCall(VecCopy(XL2, X));
368: PetscCall(VecAXPBYPCZ(X, ul, u, 1.0, WL, W));
369: }
370: // Compute the next right plane rotation P{k-1,k+1}
371: gamal_tmp = gamal;
372: SymOrtho(gamal, eplnn, &cr2, &sr2, &gamal);
374: //Store quantities for transferring from MINRES to MINRES-QLP
375: gamal_QLP = gamal_tmp;
376: vepln_QLP = vepln;
377: gama_QLP = gama;
378: ul_QLP = ul;
379: u_QLP = u;
381: // Estimate various norms
382: abs_gama = PetscAbsReal(gama);
383: Anorml = Anorm;
384: Anorm = PetscMax(PetscMax(Anorm, pnorm), PetscMax(gamal, abs_gama));
385: if (ksp->its == 1) {
386: gmin = gama;
387: gminl = gmin;
388: } else {
389: gminl2 = gminl;
390: gminl = gmin;
391: gmin = PetscMin(gminl2, PetscMin(gamal, abs_gama));
392: }
393: Acondl = Acond;
394: Acond = Anorm / gmin;
395: rnorml = rnorm;
396: relresl = relres;
397: if (flag != 9) rnorm = phi;
398: relres = rnorm / (Anorm * xnorm + beta1);
399: Arnorml = rnorml * rootl;
400: relAresl = rootl / Anorm;
402: // See if any of the stopping criteria are satisfied.
403: epsx = Anorm * xnorm * eps;
404: if (flag == flag0 || flag == 9) {
405: //if (Acond >= Acondlim) flag = 7; // Huge Acond
406: if (epsx >= beta1) flag = 5; // x is an eigenvector
407: if (minres->qlp) { /* We use these indicators only if the QLP variant has been selected */
408: PetscReal t1 = 1.0 + relres;
409: PetscReal t2 = 1.0 + relAresl;
410: if (xnorm >= minres->maxxnorm) flag = 6; // xnorm exceeded its limit
411: if (t2 <= 1) flag = 4; // Accurate LS solution
412: if (t1 <= 1) flag = 3; // Accurate Ax=b solution
413: if (relAresl <= ksp->rtol) flag = 2; // Good enough LS solution
414: if (relres <= ksp->rtol) flag = 1; // Good enough Ax=b solution
415: }
416: }
418: if (flag == 2 || flag == 4 || flag == 6 || flag == 7) {
419: Acond = Acondl;
420: rnorm = rnorml;
421: relres = relresl;
422: }
424: if (minres->monitor) { /* Mimics MATLAB code with extra flag */
425: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
426: if (ksp->its == 1) PetscCall(PetscViewerASCIIPrintf(minres->viewer, " flag rnorm Arnorm Compatible LS Anorm Acond xnorm\n"));
427: 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));
428: PetscCall(PetscViewerPopFormat(minres->viewer));
429: }
431: if (ksp->normtype == KSP_NORM_PRECONDITIONED) ksp->rnorm = rnorm;
432: else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
433: PetscCall(KSP_MatMult(ksp, Amat, X, V));
434: PetscCall(VecAYPX(V, -1.0, B));
435: PetscCall(VecNorm(V, NORM_2, &ksp->rnorm));
436: }
437: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
438: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
439: PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
440: if (!ksp->reason) {
441: switch (flag) {
442: case 1:
443: case 2:
444: case 5: /* XXX */
445: ksp->reason = KSP_CONVERGED_RTOL;
446: break;
447: case 3:
448: case 4:
449: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
450: break;
451: case 6:
452: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
453: break;
454: default:
455: break;
456: }
457: }
458: if (ksp->reason) break;
459: } while (ksp->its < ksp->max_it);
461: if (minres->monitor && flag != 2 && flag != 4 && flag != 6 && flag != 7) {
462: PetscCall(VecNorm(X, NORM_2, &xnorm));
463: PetscCall(KSP_MatMult(ksp, Amat, X, R1));
464: PetscCall(VecAYPX(R1, -1.0, B));
465: PetscCall(VecNorm(R1, NORM_2, &rnorml));
466: PetscCall(KSP_MatMult(ksp, Amat, R1, R2));
467: PetscCall(VecNorm(R2, NORM_2, &Arnorml));
468: relresl = rnorml / (Anorm * xnorm + beta1);
469: relAresl = rnorml > realmin ? Arnorml / (Anorm * rnorml) : 0.0;
470: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
471: 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));
472: PetscCall(PetscViewerPopFormat(minres->viewer));
473: }
474: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: /* This was the original implementation provided by R. Scheichl */
479: static PetscErrorCode KSPSolve_MINRES_OLD(KSP ksp)
480: {
481: PetscInt i;
482: PetscScalar alpha, beta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
483: PetscScalar rho0, rho1, rho2, rho3, dp = 0.0;
484: const PetscScalar none = -1.0;
485: PetscReal np;
486: Vec X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
487: Mat Amat;
488: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
489: PetscBool diagonalscale;
490: PetscInt stored_max_it, eigs;
491: PetscScalar *e = NULL, *d = NULL;
493: PetscFunctionBegin;
494: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
495: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
497: X = ksp->vec_sol;
498: B = ksp->vec_rhs;
499: R = ksp->work[0];
500: Z = ksp->work[1];
501: U = ksp->work[2];
502: V = ksp->work[3];
503: W = ksp->work[4];
504: UOLD = ksp->work[5];
505: VOLD = ksp->work[6];
506: WOLD = ksp->work[7];
507: WOOLD = ksp->work[8];
509: PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
511: ksp->its = 0;
512: eigs = ksp->calc_sings;
513: stored_max_it = ksp->max_it;
514: if (eigs) {
515: e = minres->e;
516: d = minres->d;
517: }
519: if (!ksp->guess_zero) {
520: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
521: PetscCall(VecAYPX(R, -1.0, B));
522: } else {
523: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
524: }
525: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
526: PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
527: KSPCheckNorm(ksp, np);
528: PetscCall(VecDot(R, Z, &dp));
529: KSPCheckDot(ksp, dp);
531: if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
532: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
533: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
534: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: ksp->rnorm = 0.0;
539: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
540: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
541: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
542: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
543: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
545: dp = PetscAbsScalar(dp);
546: dp = PetscSqrtScalar(dp);
547: beta = dp; /* beta <- sqrt(r'*z) */
548: eta = beta;
549: PetscCall(VecAXPBY(V, 1.0 / beta, 0, R)); /* v <- r / beta */
550: PetscCall(VecAXPBY(U, 1.0 / beta, 0, Z)); /* u <- z / beta */
552: i = 0;
553: do {
554: ksp->its = i + 1;
556: /* Lanczos */
558: PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- A*u */
559: PetscCall(VecDot(U, R, &alpha)); /* alpha <- r'*u */
560: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
561: if (eigs) {
562: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
563: d[i] = alpha;
564: e[i] = beta;
565: }
567: if (ksp->its > 1) {
568: Vec T[2];
569: PetscScalar alphas[] = {-alpha, -beta};
570: /* r <- r - alpha v - beta v_old */
571: T[0] = V;
572: T[1] = VOLD;
573: PetscCall(VecMAXPY(R, 2, alphas, T));
574: /* z <- z - alpha u - beta u_old */
575: T[0] = U;
576: T[1] = UOLD;
577: PetscCall(VecMAXPY(Z, 2, alphas, T));
578: } else {
579: PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha v */
580: PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha u */
581: }
583: betaold = beta;
585: PetscCall(VecDot(R, Z, &dp));
586: KSPCheckDot(ksp, dp);
587: dp = PetscAbsScalar(dp);
588: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
590: /* QR factorisation */
592: coold = cold;
593: cold = c;
594: soold = sold;
595: sold = s;
597: rho0 = cold * alpha - coold * sold * betaold;
598: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
599: rho2 = sold * alpha + coold * cold * betaold;
600: rho3 = soold * betaold;
602: /* Stop if negative curvature is detected */
603: if (ksp->converged_neg_curve && PetscRealPart(cold * rho0) <= 0.0) {
604: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1=%g, gbar %g\n", (double)PetscRealPart(cold), -(double)PetscRealPart(rho0)));
605: ksp->reason = KSP_CONVERGED_NEG_CURVE;
606: break;
607: }
609: /* Givens rotation */
611: c = rho0 / rho1;
612: s = beta / rho1;
614: /* Update */
615: /* w_oold <- w_old */
616: /* w_old <- w */
617: KSPMinresSwap3(WOOLD, WOLD, W);
619: /* w <- (u - rho2 w_old - rho3 w_oold)/rho1 */
620: PetscCall(VecAXPBY(W, 1.0 / rho1, 0.0, U));
621: if (ksp->its > 1) {
622: Vec T[] = {WOLD, WOOLD};
623: PetscScalar alphas[] = {-rho2 / rho1, -rho3 / rho1};
624: PetscInt nv = (ksp->its == 2 ? 1 : 2);
626: PetscCall(VecMAXPY(W, nv, alphas, T));
627: }
629: ceta = c * eta;
630: PetscCall(VecAXPY(X, ceta, W)); /* x <- x + c eta w */
632: /*
633: when dp is really small we have either convergence or an indefinite operator so compute true
634: residual norm to check for convergence
635: */
636: if (PetscRealPart(dp) < minres->haptol) {
637: PetscCall(PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
638: PetscCall(KSP_MatMult(ksp, Amat, X, VOLD));
639: PetscCall(VecAXPY(VOLD, none, B));
640: PetscCall(VecNorm(VOLD, NORM_2, &np));
641: KSPCheckNorm(ksp, np);
642: } else {
643: /* otherwise compute new residual norm via recurrence relation */
644: np *= PetscAbsScalar(s);
645: }
647: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
648: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
649: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
650: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
651: if (ksp->reason) break;
653: if (PetscRealPart(dp) < minres->haptol) {
654: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
655: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
656: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
657: break;
658: }
660: eta = -s * eta;
661: KSPMinresSwap3(VOLD, V, R);
662: KSPMinresSwap3(UOLD, U, Z);
663: PetscCall(VecScale(V, 1.0 / beta)); /* v <- r / beta */
664: PetscCall(VecScale(U, 1.0 / beta)); /* u <- z / beta */
666: i++;
667: } while (i < ksp->max_it);
668: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode KSPDestroy_MINRES(KSP ksp)
673: {
674: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
676: PetscFunctionBegin;
677: PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
678: PetscCall(PetscViewerDestroy(&minres->viewer));
679: PetscCall(PetscFree(ksp->data));
680: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", NULL));
681: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", NULL));
682: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", NULL));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: static PetscErrorCode KSPMINRESSetUseQLP_MINRES(KSP ksp, PetscBool qlp)
687: {
688: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
690: PetscFunctionBegin;
691: minres->qlp = qlp;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: static PetscErrorCode KSPMINRESSetRadius_MINRES(KSP ksp, PetscReal radius)
696: {
697: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
699: PetscFunctionBegin;
700: minres->maxxnorm = radius;
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: static PetscErrorCode KSPMINRESGetUseQLP_MINRES(KSP ksp, PetscBool *qlp)
705: {
706: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
708: PetscFunctionBegin;
709: *qlp = minres->qlp;
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode KSPSetFromOptions_MINRES(KSP ksp, PetscOptionItems PetscOptionsObject)
714: {
715: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
717: PetscFunctionBegin;
718: PetscOptionsHeadBegin(PetscOptionsObject, "KSP MINRES options");
719: { /* Allow comparing with the old code (to be removed in a few releases) */
720: PetscBool flg = PETSC_FALSE;
721: PetscCall(PetscOptionsBool("-ksp_minres_old", "Use old implementation (to be removed)", "None", flg, &flg, NULL));
722: if (flg) ksp->ops->solve = KSPSolve_MINRES_OLD;
723: else ksp->ops->solve = KSPSolve_MINRES;
724: }
725: PetscCall(PetscOptionsBool("-ksp_minres_qlp", "Solve with QLP variant", "KSPMINRESSetUseQLP", minres->qlp, &minres->qlp, NULL));
726: PetscCall(PetscOptionsReal("-ksp_minres_radius", "Maximum allowed norm of solution", "KSPMINRESSetRadius", minres->maxxnorm, &minres->maxxnorm, NULL));
727: PetscCall(PetscOptionsReal("-ksp_minres_trancond", "Threshold on condition number to dynamically switch to QLP", "None", minres->TranCond, &minres->TranCond, NULL));
728: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ksp), PetscOptionsObject->options, PetscOptionsObject->prefix, "-ksp_minres_monitor", &minres->viewer, &minres->viewer_fmt, &minres->monitor));
729: PetscCall(PetscOptionsReal("-ksp_minres_nutol", "Inexactness tolerance", NULL, minres->nutol, &minres->nutol, NULL));
730: PetscOptionsHeadEnd();
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /*@
735: KSPMINRESSetUseQLP - Use the QLP variant of `KSPMINRES`
737: Logically Collective
739: Input Parameters:
740: + ksp - the iterative context
741: - qlp - a Boolean indicating if the QLP variant should be used
743: Level: beginner
745: Note:
746: By default, the QLP variant is not used.
748: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESGetUseQLP()`
749: @*/
750: PetscErrorCode KSPMINRESSetUseQLP(KSP ksp, PetscBool qlp)
751: {
752: PetscFunctionBegin;
755: PetscTryMethod(ksp, "KSPMINRESSetUseQLP_C", (KSP, PetscBool), (ksp, qlp));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: /*@
760: KSPMINRESSetRadius - Set the maximum solution norm allowed for use with trust region methods
762: Logically Collective
764: Input Parameters:
765: + ksp - the iterative context
766: - radius - the value
768: Level: beginner
770: Options Database Key:
771: . -ksp_minres_radius <real> - maximum allowed solution norm
773: Developer Note:
774: Perhaps the KSPXXXSetRadius() should be unified
776: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
777: @*/
778: PetscErrorCode KSPMINRESSetRadius(KSP ksp, PetscReal radius)
779: {
780: PetscFunctionBegin;
783: PetscTryMethod(ksp, "KSPMINRESSetRadius_C", (KSP, PetscReal), (ksp, radius));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: /*@
788: KSPMINRESGetUseQLP - Get the flag that indicates if the QLP variant is being used
790: Logically Collective
792: Input Parameter:
793: . ksp - the iterative context
795: Output Parameter:
796: . qlp - a Boolean indicating if the QLP variant is used
798: Level: beginner
800: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
801: @*/
802: PetscErrorCode KSPMINRESGetUseQLP(KSP ksp, PetscBool *qlp)
803: {
804: PetscFunctionBegin;
806: PetscAssertPointer(qlp, 2);
807: PetscUseMethod(ksp, "KSPMINRESGetUseQLP_C", (KSP, PetscBool *), (ksp, qlp));
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: /*MC
812: KSPMINRES - This code implements the MINRES (Minimum Residual) method and its QLP variant {cite}`paige.saunders:solution`, {cite}`choi2011minres`,
813: {cite}`liu2022newton` for solving linear systems using `KSP`.
815: Options Database Keys:
816: + -ksp_minres_qlp <bool> - activates QLP code
817: . -ksp_minres_radius <real> - maximum allowed solution norm
818: . -ksp_minres_trancond <real> - threshold on condition number to dynamically switch to QLP iterations when QLP has been activated
819: . -ksp_minres_monitor - monitors convergence quantities
820: - -ksp_minres_nutol <real> - inexactness tolerance (see https://arxiv.org/pdf/2208.07095.pdf)
822: Level: beginner
824: Notes:
825: The matrix (operator) and the preconditioner must be symmetric and the preconditioner must also be positive definite for this method.
827: `KSPMINRES` is often the best Krylov method for symmetric indefinite matrices.
829: Supports only left preconditioning.
831: Contributed by:
832: Original MINRES code - Robert Scheichl: maprs@maths.bath.ac.uk
833: QLP variant adapted from: https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
835: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`, `KSPMINRESGetUseQLP()`, `KSPMINRESSetUseQLP()`, `KSPMINRESSetRadius()`
836: M*/
837: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
838: {
839: KSP_MINRES *minres;
841: PetscFunctionBegin;
842: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
843: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
844: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
845: PetscCall(PetscNew(&minres));
847: /* this parameter is arbitrary and belongs to the old implementation; but e-50 didn't work for __float128 in one example */
848: #if defined(PETSC_USE_REAL___FLOAT128)
849: minres->haptol = 1.e-100;
850: #elif defined(PETSC_USE_REAL_SINGLE)
851: minres->haptol = 1.e-25;
852: #else
853: minres->haptol = 1.e-50;
854: #endif
855: /* those are set as 1.e7 in the MATLAB code -> use 1.0/sqrt(eps) to support single precision */
856: minres->maxxnorm = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
857: minres->TranCond = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
859: ksp->data = (void *)minres;
861: ksp->ops->setup = KSPSetUp_MINRES;
862: ksp->ops->solve = KSPSolve_MINRES;
863: ksp->ops->destroy = KSPDestroy_MINRES;
864: ksp->ops->setfromoptions = KSPSetFromOptions_MINRES;
865: ksp->ops->buildsolution = KSPBuildSolutionDefault;
866: ksp->ops->buildresidual = KSPBuildResidualDefault;
868: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", KSPMINRESSetRadius_MINRES));
869: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", KSPMINRESSetUseQLP_MINRES));
870: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", KSPMINRESGetUseQLP_MINRES));
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: PetscErrorCode KSPComputeEigenvalues_MINRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
875: {
876: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
877: PetscScalar *d, *e;
878: PetscReal *ee;
879: PetscInt n = ksp->its;
880: PetscBLASInt bn, lierr = 0, ldz = 1;
882: PetscFunctionBegin;
883: PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
884: *neig = n;
886: PetscCall(PetscArrayzero(c, nmax));
887: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
888: d = minres->d;
889: e = minres->e;
890: ee = minres->ee;
892: /* copy tridiagonal matrix to work space */
893: for (PetscInt j = 0; j < n; j++) {
894: r[j] = PetscRealPart(d[j]);
895: ee[j] = PetscRealPart(e[j + 1]);
896: }
898: PetscCall(PetscBLASIntCast(n, &bn));
899: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
900: PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, ee, NULL, &ldz, NULL, &lierr));
901: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
902: PetscCall(PetscFPTrapPop());
903: PetscCall(PetscSortReal(n, r));
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP ksp, PetscReal *emax, PetscReal *emin)
908: {
909: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
910: PetscScalar *d, *e;
911: PetscReal *dd, *ee;
912: PetscInt n = ksp->its;
913: PetscBLASInt bn, lierr = 0, ldz = 1;
915: PetscFunctionBegin;
916: if (!n) {
917: *emax = *emin = 1.0;
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
920: d = minres->d;
921: e = minres->e;
922: dd = minres->dd;
923: ee = minres->ee;
925: /* copy tridiagonal matrix to work space */
926: for (PetscInt j = 0; j < n; j++) {
927: dd[j] = PetscRealPart(d[j]);
928: ee[j] = PetscRealPart(e[j + 1]);
929: }
931: PetscCall(PetscBLASIntCast(n, &bn));
932: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
933: PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, ee, NULL, &ldz, NULL, &lierr));
934: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
935: PetscCall(PetscFPTrapPop());
936: for (PetscInt j = 0; j < n; j++) dd[j] = PetscAbsReal(dd[j]);
937: PetscCall(PetscSortReal(n, dd));
938: *emin = dd[0];
939: *emax = dd[n - 1];
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }