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