1 #include <petsc/private/kspimpl.h>
2
3 typedef struct {
4 PetscReal haptol;
5 } KSP_SYMMLQ;
6
KSPSetUp_SYMMLQ(KSP ksp)7 static PetscErrorCode KSPSetUp_SYMMLQ(KSP ksp)
8 {
9 PetscFunctionBegin;
10 PetscCall(KSPSetWorkVecs(ksp, 9));
11 PetscFunctionReturn(PETSC_SUCCESS);
12 }
13
KSPSolve_SYMMLQ(KSP ksp)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;
25
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);
29
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];
40
41 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
42
43 ksp->its = 0;
44
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 }
55
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 }
65
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);
76
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);
92
93 i = 0;
94 ceta = 0.;
95 do {
96 ksp->its = i + 1;
97
98 /* Update */
99 if (ksp->its > 1) {
100 PetscCall(VecCopy(V, VOLD)); /* v_old <- v; */
101 PetscCall(VecCopy(U, UOLD)); /* u_old <- u; */
102
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; */
107
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) */
114
115 ceta_oold = ceta_old;
116 ceta_old = ceta;
117 }
118
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; */
123
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 }
135
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); */
143
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 */
153
154 /* Givens rotation: [c -s; s c] (different from the Reference!) */
155 c = rho0 / rho1;
156 s = beta / rho1;
157
158 if (ksp->its == 1) ceta = beta1 / rho1;
159 else ceta = -(rho2 * ceta_old + rho3 * ceta_oold) / rho1;
160
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) */
164
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);
173
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;
177
178 PetscCall(VecAXPY(X, ceta_bar, Wbar)); /* x <- x + ceta_bar*w_bar */
179
180 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
184 /*MC
185 KSPSYMMLQ - Implements the SYMMLQ method {cite}`paige.saunders:solution`.
186
187 Level: beginner
188
189 Notes:
190 The operator and the preconditioner must be symmetric for this method.
191
192 The preconditioner must be POSITIVE-DEFINITE.
193
194 Supports only left preconditioning.
195
196 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`
197 M*/
KSPCreate_SYMMLQ(KSP ksp)198 PETSC_EXTERN PetscErrorCode KSPCreate_SYMMLQ(KSP ksp)
199 {
200 KSP_SYMMLQ *symmlq;
201
202 PetscFunctionBegin;
203 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
204 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
205
206 PetscCall(PetscNew(&symmlq));
207 symmlq->haptol = 1.e-18;
208 ksp->data = (void *)symmlq;
209
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 }
222