Actual source code: minres.c
2: #include <petsc/private/kspimpl.h>
4: typedef struct {
5: PetscReal haptol;
6: } KSP_MINRES;
8: static PetscErrorCode KSPSetUp_MINRES(KSP ksp)
9: {
10: KSPSetWorkVecs(ksp, 9);
11: return 0;
12: }
14: static PetscErrorCode KSPSolve_MINRES(KSP ksp)
15: {
16: PetscInt i;
17: PetscScalar alpha, beta, ibeta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
18: PetscScalar rho0, rho1, irho1, rho2, rho3, dp = 0.0;
19: const PetscScalar none = -1.0;
20: PetscReal np;
21: Vec X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
22: Mat Amat, Pmat;
23: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
24: PetscBool diagonalscale;
26: PCGetDiagonalScale(ksp->pc, &diagonalscale);
29: X = ksp->vec_sol;
30: B = ksp->vec_rhs;
31: R = ksp->work[0];
32: Z = ksp->work[1];
33: U = ksp->work[2];
34: V = ksp->work[3];
35: W = ksp->work[4];
36: UOLD = ksp->work[5];
37: VOLD = ksp->work[6];
38: WOLD = ksp->work[7];
39: WOOLD = ksp->work[8];
41: PCGetOperators(ksp->pc, &Amat, &Pmat);
43: ksp->its = 0;
45: VecSet(UOLD, 0.0); /* u_old <- 0 */
46: VecSet(VOLD, 0.0); /* v_old <- 0 */
47: VecSet(W, 0.0); /* w <- 0 */
48: VecSet(WOLD, 0.0); /* w_old <- 0 */
50: if (!ksp->guess_zero) {
51: KSP_MatMult(ksp, Amat, X, R); /* r <- b - A*x */
52: VecAYPX(R, -1.0, B);
53: } else {
54: VecCopy(B, R); /* r <- b (x is 0) */
55: }
56: KSP_PCApply(ksp, R, Z); /* z <- B*r */
57: VecNorm(Z, NORM_2, &np); /* np <- ||z|| */
58: KSPCheckNorm(ksp, np);
59: VecDot(R, Z, &dp);
60: KSPCheckDot(ksp, dp);
62: if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
64: PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol);
65: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
66: return 0;
67: }
69: ksp->rnorm = 0.0;
70: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
71: KSPLogResidualHistory(ksp, ksp->rnorm);
72: KSPMonitor(ksp, 0, ksp->rnorm);
73: (*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP); /* test for convergence */
74: if (ksp->reason) return 0;
76: dp = PetscAbsScalar(dp);
77: dp = PetscSqrtScalar(dp);
78: beta = dp; /* beta <- sqrt(r'*z */
79: eta = beta;
81: VecCopy(R, V);
82: VecCopy(Z, U);
83: ibeta = 1.0 / beta;
84: VecScale(V, ibeta); /* v <- r / beta */
85: VecScale(U, ibeta); /* u <- z / beta */
87: i = 0;
88: do {
89: ksp->its = i + 1;
91: /* Lanczos */
93: KSP_MatMult(ksp, Amat, U, R); /* r <- A*u */
94: VecDot(U, R, &alpha); /* alpha <- r'*u */
95: KSP_PCApply(ksp, R, Z); /* z <- B*r */
97: VecAXPY(R, -alpha, V); /* r <- r - alpha v */
98: VecAXPY(Z, -alpha, U); /* z <- z - alpha u */
99: VecAXPY(R, -beta, VOLD); /* r <- r - beta v_old */
100: VecAXPY(Z, -beta, UOLD); /* z <- z - beta u_old */
102: betaold = beta;
104: VecDot(R, Z, &dp);
105: KSPCheckDot(ksp, dp);
106: dp = PetscAbsScalar(dp);
107: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
109: /* QR factorisation */
111: coold = cold;
112: cold = c;
113: soold = sold;
114: sold = s;
116: rho0 = cold * alpha - coold * sold * betaold;
117: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
118: rho2 = sold * alpha + coold * cold * betaold;
119: rho3 = soold * betaold;
121: /* Givens rotation */
123: c = rho0 / rho1;
124: s = beta / rho1;
126: /* Update */
128: VecCopy(WOLD, WOOLD); /* w_oold <- w_old */
129: VecCopy(W, WOLD); /* w_old <- w */
131: VecCopy(U, W); /* w <- u */
132: VecAXPY(W, -rho2, WOLD); /* w <- w - rho2 w_old */
133: VecAXPY(W, -rho3, WOOLD); /* w <- w - rho3 w_oold */
134: irho1 = 1.0 / rho1;
135: VecScale(W, irho1); /* w <- w / rho1 */
137: ceta = c * eta;
138: VecAXPY(X, ceta, W); /* x <- x + c eta w */
140: /*
141: when dp is really small we have either convergence or an indefinite operator so compute true
142: residual norm to check for convergence
143: */
144: if (PetscRealPart(dp) < minres->haptol) {
145: PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol);
146: KSP_MatMult(ksp, Amat, X, VOLD);
147: VecAXPY(VOLD, none, B);
148: VecNorm(VOLD, NORM_2, &np);
149: KSPCheckNorm(ksp, np);
150: } else {
151: /* otherwise compute new residual norm via recurrence relation */
152: np *= PetscAbsScalar(s);
153: }
155: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
156: KSPLogResidualHistory(ksp, ksp->rnorm);
157: KSPMonitor(ksp, i + 1, ksp->rnorm);
158: (*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP); /* test for convergence */
159: if (ksp->reason) break;
161: if (PetscRealPart(dp) < minres->haptol) {
163: PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol);
164: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
165: break;
166: }
168: eta = -s * eta;
169: VecCopy(V, VOLD);
170: VecCopy(U, UOLD);
171: VecCopy(R, V);
172: VecCopy(Z, U);
173: ibeta = 1.0 / beta;
174: VecScale(V, ibeta); /* v <- r / beta */
175: VecScale(U, ibeta); /* u <- z / beta */
177: i++;
178: } while (i < ksp->max_it);
179: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
180: return 0;
181: }
183: /*MC
184: KSPMINRES - This code implements the MINRES (Minimum Residual) method.
186: Level: beginner
188: Notes:
189: The operator and the preconditioner must be symmetric and the preconditioner must be positive definite for this method.
191: Supports only left preconditioning.
193: Reference:
194: . * - Paige & Saunders, Solution of sparse indefinite systems of linear equations, SIAM J. Numer. Anal. 12, 1975.
196: Contributed by:
197: Robert Scheichl: maprs@maths.bath.ac.uk
199: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`
200: M*/
201: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
202: {
203: KSP_MINRES *minres;
205: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
206: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
207: PetscNew(&minres);
209: /* this parameter is arbitrary; but e-50 didn't work for __float128 in one example */
210: #if defined(PETSC_USE_REAL___FLOAT128)
211: minres->haptol = 1.e-100;
212: #elif defined(PETSC_USE_REAL_SINGLE)
213: minres->haptol = 1.e-25;
214: #else
215: minres->haptol = 1.e-50;
216: #endif
217: ksp->data = (void *)minres;
219: /*
220: Sets the functions that are associated with this data structure
221: (in C++ this is the same as defining virtual functions)
222: */
223: ksp->ops->setup = KSPSetUp_MINRES;
224: ksp->ops->solve = KSPSolve_MINRES;
225: ksp->ops->destroy = KSPDestroyDefault;
226: ksp->ops->setfromoptions = NULL;
227: ksp->ops->buildsolution = KSPBuildSolutionDefault;
228: ksp->ops->buildresidual = KSPBuildResidualDefault;
229: return 0;
230: }