Actual source code: symmlq.c
2: #include <petsc/private/kspimpl.h>
4: typedef struct {
5: PetscReal haptol;
6: } KSP_SYMMLQ;
8: PetscErrorCode KSPSetUp_SYMMLQ(KSP ksp)
9: {
10: PetscFunctionBegin;
11: PetscCall(KSPSetWorkVecs(ksp, 9));
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
15: PetscErrorCode KSPSolve_SYMMLQ(KSP ksp)
16: {
17: PetscInt i;
18: PetscScalar alpha, beta, ibeta, betaold, beta1, ceta = 0, ceta_oold = 0.0, ceta_old = 0.0, ceta_bar;
19: PetscScalar c = 1.0, cold = 1.0, s = 0.0, sold = 0.0, coold, soold, rho0, rho1, rho2, rho3;
20: PetscScalar dp = 0.0;
21: PetscReal np = 0.0, s_prod;
22: Vec X, B, R, Z, U, V, W, UOLD, VOLD, Wbar;
23: Mat Amat, Pmat;
24: KSP_SYMMLQ *symmlq = (KSP_SYMMLQ *)ksp->data;
25: PetscBool diagonalscale;
27: PetscFunctionBegin;
28: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
29: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
31: X = ksp->vec_sol;
32: B = ksp->vec_rhs;
33: R = ksp->work[0];
34: Z = ksp->work[1];
35: U = ksp->work[2];
36: V = ksp->work[3];
37: W = ksp->work[4];
38: UOLD = ksp->work[5];
39: VOLD = ksp->work[6];
40: Wbar = ksp->work[7];
42: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
44: ksp->its = 0;
46: PetscCall(VecSet(UOLD, 0.0)); /* u_old <- zeros; */
47: PetscCall(VecCopy(UOLD, VOLD)); /* v_old <- u_old; */
48: PetscCall(VecCopy(UOLD, W)); /* w <- u_old; */
49: PetscCall(VecCopy(UOLD, Wbar)); /* w_bar <- u_old; */
50: if (!ksp->guess_zero) {
51: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
52: PetscCall(VecAYPX(R, -1.0, B));
53: } else {
54: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
55: }
57: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
58: PetscCall(VecDot(R, Z, &dp)); /* dp = r'*z; */
59: KSPCheckDot(ksp, dp);
60: if (PetscAbsScalar(dp) < symmlq->haptol) {
61: PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol));
62: ksp->rnorm = 0.0; /* what should we really put here? */
63: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN; /* bugfix proposed by Lourens (lourens.vanzanen@shell.com) */
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: #if !defined(PETSC_USE_COMPLEX)
68: if (dp < 0.0) {
69: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
72: #endif
73: dp = PetscSqrtScalar(dp);
74: beta = dp; /* beta <- sqrt(r'*z) */
75: beta1 = beta;
76: s_prod = PetscAbsScalar(beta1);
78: PetscCall(VecCopy(R, V)); /* v <- r; */
79: PetscCall(VecCopy(Z, U)); /* u <- z; */
80: ibeta = 1.0 / beta;
81: PetscCall(VecScale(V, ibeta)); /* v <- ibeta*v; */
82: PetscCall(VecScale(U, ibeta)); /* u <- ibeta*u; */
83: PetscCall(VecCopy(U, Wbar)); /* w_bar <- u; */
84: if (ksp->normtype != KSP_NORM_NONE) {
85: PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
86: KSPCheckNorm(ksp, np);
87: }
88: PetscCall(KSPLogResidualHistory(ksp, np));
89: PetscCall(KSPMonitor(ksp, 0, np));
90: ksp->rnorm = np;
91: PetscCall((*ksp->converged)(ksp, 0, np, &ksp->reason, ksp->cnvP)); /* test for convergence */
92: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
94: i = 0;
95: ceta = 0.;
96: do {
97: ksp->its = i + 1;
99: /* Update */
100: if (ksp->its > 1) {
101: PetscCall(VecCopy(V, VOLD)); /* v_old <- v; */
102: PetscCall(VecCopy(U, UOLD)); /* u_old <- u; */
104: PetscCall(VecCopy(R, V));
105: PetscCall(VecScale(V, 1.0 / beta)); /* v <- ibeta*r; */
106: PetscCall(VecCopy(Z, U));
107: PetscCall(VecScale(U, 1.0 / beta)); /* u <- ibeta*z; */
109: PetscCall(VecCopy(Wbar, W));
110: PetscCall(VecScale(W, c));
111: PetscCall(VecAXPY(W, s, U)); /* w <- c*w_bar + s*u; (w_k) */
112: PetscCall(VecScale(Wbar, -s));
113: PetscCall(VecAXPY(Wbar, c, U)); /* w_bar <- -s*w_bar + c*u; (w_bar_(k+1)) */
114: PetscCall(VecAXPY(X, ceta, W)); /* x <- x + ceta * w; (xL_k) */
116: ceta_oold = ceta_old;
117: ceta_old = ceta;
118: }
120: /* Lanczos */
121: PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- Amat*u; */
122: PetscCall(VecDot(U, R, &alpha)); /* alpha <- u'*r; */
123: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r; */
125: PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha* v; */
126: PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha* u; */
127: PetscCall(VecAXPY(R, -beta, VOLD)); /* r <- r - beta * v_old; */
128: PetscCall(VecAXPY(Z, -beta, UOLD)); /* z <- z - beta * u_old; */
129: betaold = beta; /* beta_k */
130: PetscCall(VecDot(R, Z, &dp)); /* dp <- r'*z; */
131: KSPCheckDot(ksp, dp);
132: if (PetscAbsScalar(dp) < symmlq->haptol) {
133: PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol));
134: dp = 0.0;
135: }
137: #if !defined(PETSC_USE_COMPLEX)
138: if (dp < 0.0) {
139: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
140: break;
141: }
142: #endif
143: beta = PetscSqrtScalar(dp); /* beta = sqrt(dp); */
145: /* QR factorization */
146: coold = cold;
147: cold = c;
148: soold = sold;
149: sold = s;
150: rho0 = cold * alpha - coold * sold * betaold; /* gamma_bar */
151: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta); /* gamma */
152: rho2 = sold * alpha + coold * cold * betaold; /* delta */
153: rho3 = soold * betaold; /* epsilon */
155: /* Givens rotation: [c -s; s c] (different from the Reference!) */
156: c = rho0 / rho1;
157: s = beta / rho1;
159: if (ksp->its == 1) ceta = beta1 / rho1;
160: else ceta = -(rho2 * ceta_old + rho3 * ceta_oold) / rho1;
162: s_prod = s_prod * PetscAbsScalar(s);
163: if (c == 0.0) np = s_prod * 1.e16;
164: else np = s_prod / PetscAbsScalar(c); /* residual norm for xc_k (CGNORM) */
166: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
167: else ksp->rnorm = 0.0;
168: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
169: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
170: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
171: if (ksp->reason) break;
172: i++;
173: } while (i < ksp->max_it);
175: /* move to the CG point: xc_(k+1) */
176: if (c == 0.0) ceta_bar = ceta * 1.e15;
177: else ceta_bar = ceta / c;
179: PetscCall(VecAXPY(X, ceta_bar, Wbar)); /* x <- x + ceta_bar*w_bar */
181: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*MC
186: KSPSYMMLQ - This code implements the SYMMLQ method.
188: Level: beginner
190: Notes:
191: The operator and the preconditioner must be symmetric for this method.
193: The preconditioner must be POSITIVE-DEFINITE.
195: Supports only left preconditioning.
197: Reference:
198: . * - Paige & Saunders, Solution of sparse indefinite systems of linear equations, SIAM J. Numer. Anal. 12, 1975.
200: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`
201: M*/
202: PETSC_EXTERN PetscErrorCode KSPCreate_SYMMLQ(KSP ksp)
203: {
204: KSP_SYMMLQ *symmlq;
206: PetscFunctionBegin;
207: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
208: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
210: PetscCall(PetscNew(&symmlq));
211: symmlq->haptol = 1.e-18;
212: ksp->data = (void *)symmlq;
214: /*
215: Sets the functions that are associated with this data structure
216: (in C++ this is the same as defining virtual functions)
217: */
218: ksp->ops->setup = KSPSetUp_SYMMLQ;
219: ksp->ops->solve = KSPSolve_SYMMLQ;
220: ksp->ops->destroy = KSPDestroyDefault;
221: ksp->ops->setfromoptions = NULL;
222: ksp->ops->buildsolution = KSPBuildSolutionDefault;
223: ksp->ops->buildresidual = KSPBuildResidualDefault;
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }