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;
25: PetscFunctionBegin;
26: X = ksp->vec_sol;
27: B = ksp->vec_rhs;
28: R = ksp->work[0];
29: Z = ksp->work[1];
30: U = ksp->work[2];
31: V = ksp->work[3];
32: W = ksp->work[4];
33: UOLD = ksp->work[5];
34: VOLD = ksp->work[6];
35: Wbar = ksp->work[7];
37: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
39: ksp->its = 0;
41: PetscCall(VecSet(UOLD, 0.0)); /* u_old <- zeros; */
42: PetscCall(VecCopy(UOLD, VOLD)); /* v_old <- u_old; */
43: PetscCall(VecCopy(UOLD, W)); /* w <- u_old; */
44: PetscCall(VecCopy(UOLD, Wbar)); /* w_bar <- u_old; */
45: if (!ksp->guess_zero) {
46: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
47: PetscCall(VecAYPX(R, -1.0, B));
48: } else {
49: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
50: }
52: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
53: PetscCall(VecDot(R, Z, &dp)); /* dp = r'*z; */
54: KSPCheckDot(ksp, dp);
55: if (PetscAbsScalar(dp) < symmlq->haptol) {
56: PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol));
57: ksp->rnorm = 0.0; /* what should we really put here? */
58: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN; /* bugfix proposed by Lourens (lourens.vanzanen@shell.com) */
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: #if !PetscDefined(USE_COMPLEX)
63: if (dp < 0.0) {
64: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
67: #endif
68: dp = PetscSqrtScalar(dp);
69: beta = dp; /* beta <- sqrt(r'*z) */
70: beta1 = beta;
71: s_prod = PetscAbsScalar(beta1);
73: PetscCall(VecCopy(R, V)); /* v <- r; */
74: PetscCall(VecCopy(Z, U)); /* u <- z; */
75: ibeta = 1.0 / beta;
76: PetscCall(VecScale(V, ibeta)); /* v <- ibeta*v; */
77: PetscCall(VecScale(U, ibeta)); /* u <- ibeta*u; */
78: PetscCall(VecCopy(U, Wbar)); /* w_bar <- u; */
79: if (ksp->normtype != KSP_NORM_NONE) {
80: PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
81: KSPCheckNorm(ksp, np);
82: }
83: PetscCall(KSPLogResidualHistory(ksp, np));
84: PetscCall(KSPMonitor(ksp, 0, np));
85: ksp->rnorm = np;
86: PetscCall((*ksp->converged)(ksp, 0, np, &ksp->reason, ksp->cnvP)); /* test for convergence */
87: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
89: i = 0;
90: ceta = 0.;
91: do {
92: ksp->its = i + 1;
94: /* Update */
95: if (ksp->its > 1) {
96: PetscCall(VecCopy(V, VOLD)); /* v_old <- v; */
97: PetscCall(VecCopy(U, UOLD)); /* u_old <- u; */
99: PetscCall(VecCopy(R, V));
100: PetscCall(VecScale(V, 1.0 / beta)); /* v <- ibeta*r; */
101: PetscCall(VecCopy(Z, U));
102: PetscCall(VecScale(U, 1.0 / beta)); /* u <- ibeta*z; */
104: PetscCall(VecCopy(Wbar, W));
105: PetscCall(VecScale(W, c));
106: PetscCall(VecAXPY(W, s, U)); /* w <- c*w_bar + s*u; (w_k) */
107: PetscCall(VecScale(Wbar, -s));
108: PetscCall(VecAXPY(Wbar, c, U)); /* w_bar <- -s*w_bar + c*u; (w_bar_(k+1)) */
109: PetscCall(VecAXPY(X, ceta, W)); /* x <- x + ceta * w; (xL_k) */
111: ceta_oold = ceta_old;
112: ceta_old = ceta;
113: }
115: /* Lanczos */
116: PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- Amat*u; */
117: PetscCall(VecDot(U, R, &alpha)); /* alpha <- u'*r; */
118: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r; */
120: PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha* v; */
121: PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha* u; */
122: PetscCall(VecAXPY(R, -beta, VOLD)); /* r <- r - beta * v_old; */
123: PetscCall(VecAXPY(Z, -beta, UOLD)); /* z <- z - beta * u_old; */
124: betaold = beta; /* beta_k */
125: PetscCall(VecDot(R, Z, &dp)); /* dp <- r'*z; */
126: KSPCheckDot(ksp, dp);
127: if (PetscAbsScalar(dp) < symmlq->haptol) {
128: PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol));
129: dp = 0.0;
130: }
132: #if !PetscDefined(USE_COMPLEX)
133: if (dp < 0.0) {
134: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
135: break;
136: }
137: #endif
138: beta = PetscSqrtScalar(dp); /* beta = sqrt(dp); */
140: /* QR factorization */
141: coold = cold;
142: cold = c;
143: soold = sold;
144: sold = s;
145: rho0 = cold * alpha - coold * sold * betaold; /* gamma_bar */
146: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta); /* gamma */
147: rho2 = sold * alpha + coold * cold * betaold; /* delta */
148: rho3 = soold * betaold; /* epsilon */
150: /* Givens rotation: [c -s; s c] (different from the Reference!) */
151: c = rho0 / rho1;
152: s = beta / rho1;
154: if (ksp->its == 1) ceta = beta1 / rho1;
155: else ceta = -(rho2 * ceta_old + rho3 * ceta_oold) / rho1;
157: s_prod = s_prod * PetscAbsScalar(s);
158: if (c == 0.0) np = s_prod * 1.e16;
159: else np = s_prod / PetscAbsScalar(c); /* residual norm for xc_k (CGNORM) */
161: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
162: else ksp->rnorm = 0.0;
163: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
164: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
165: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
166: if (ksp->reason) break;
167: i++;
168: } while (i < ksp->max_it);
170: /* move to the CG point: xc_(k+1) */
171: if (c == 0.0) ceta_bar = ceta * 1.e15;
172: else ceta_bar = ceta / c;
174: PetscCall(VecAXPY(X, ceta_bar, Wbar)); /* x <- x + ceta_bar*w_bar */
176: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*MC
181: KSPSYMMLQ - Implements the SYMMLQ method {cite}`paige.saunders:solution`.
183: Level: beginner
185: Notes:
186: The operator and the preconditioner must be symmetric for this method.
188: The preconditioner must be POSITIVE-DEFINITE.
190: Supports only left preconditioning.
192: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`
193: M*/
194: PETSC_EXTERN PetscErrorCode KSPCreate_SYMMLQ(KSP ksp)
195: {
196: KSP_SYMMLQ *symmlq;
198: PetscFunctionBegin;
199: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
200: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
202: PetscCall(PetscNew(&symmlq));
203: symmlq->haptol = 1.e-18;
204: ksp->data = (void *)symmlq;
206: /*
207: Sets the functions that are associated with this data structure
208: (in C++ this is the same as defining virtual functions)
209: */
210: ksp->ops->setup = KSPSetUp_SYMMLQ;
211: ksp->ops->solve = KSPSolve_SYMMLQ;
212: ksp->ops->destroy = KSPDestroyDefault;
213: ksp->ops->setfromoptions = NULL;
214: ksp->ops->buildsolution = KSPBuildSolutionDefault;
215: ksp->ops->buildresidual = KSPBuildResidualDefault;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }