Actual source code: orthog.c
1: /*
2: The KSP orthogonalization routines, used in GMRES and other solvers.
3: */
4: #include <petsc/private/kspimpl.h>
6: /*@
7: KSPOrthogonalizationSet - Sets the orthogonalization routine used by `KSPGMRES` and other solvers.
9: Logically Collective
11: Input Parameters:
12: + ksp - the Krylov space solver context
13: - orthog - orthogonalization function; see `KSPOrthogonalizationFn` for the calling sequence
15: Options Database Key:
16: . -ksp_orthogonalization (cgs|mgs) - choose between classical (default) or modified Gram-Schmidt for orthogonalization
18: Level: intermediate
20: Notes:
21: This is used by solvers that explicitly orthogonalize a set of vectors, such as the `KSPGMRES` variants
22: and `KSPIDR`; it has no effect on methods such as `KSPBCGS` that never call it.
24: Two orthogonalization routines are predefined, `KSPOrthogonalizationModifiedGramSchmidt()` and the default
25: `KSPOrthogonalizationClassicalGramSchmidt()`.
27: Use `KSPOrthogonalizationSetCGSRefinementType()` to determine if iterative refinement is used to increase stability.
29: .seealso: [](ch_ksp), `KSPOrthogonalizationFn`, `KSPOrthogonalizationGet()`, `KSPOrthogonalizationSetCGSRefinementType()`, `KSPOrthogonalizationClassicalGramSchmidt()`, `KSPOrthogonalizationModifiedGramSchmidt()`
30: @*/
31: PetscErrorCode KSPOrthogonalizationSet(KSP ksp, KSPOrthogonalizationFn *orthog)
32: {
33: PetscFunctionBegin;
36: ksp->orthog = orthog;
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: /*@
41: KSPOrthogonalizationGet - Gets the orthogonalization routine used by `KSPGMRES` and other solvers.
43: Not Collective
45: Input Parameter:
46: . ksp - the Krylov space solver context
48: Output Parameter:
49: . orthog - orthogonalization function; see `KSPOrthogonalizationFn` for the calling sequence
51: Level: intermediate
53: Notes:
54: Two orthogonalization routines are predefined, `KSPOrthogonalizationModifiedGramSchmidt()` and the default
55: `KSPOrthogonalizationClassicalGramSchmidt()`.
57: Use `KSPOrthogonalizationSetCGSRefinementType()` to determine if iterative refinement is used to increase stability.
59: .seealso: [](ch_ksp), `KSPOrthogonalizationFn`, `KSPOrthogonalizationSetCGSRefinementType()`, `KSPOrthogonalizationClassicalGramSchmidt()`, `KSPOrthogonalizationModifiedGramSchmidt()`
60: @*/
61: PetscErrorCode KSPOrthogonalizationGet(KSP ksp, KSPOrthogonalizationFn **orthog)
62: {
63: PetscFunctionBegin;
65: PetscAssertPointer(orthog, 2);
66: *orthog = ksp->orthog;
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: /*@
71: KSPOrthogonalizationSetCGSRefinementType - Sets the type of iterative refinement to use in the classical Gram-Schmidt
72: orthogonalization used by `KSPGMRES` and other solvers.
74: Logically Collective
76: Input Parameters:
77: + ksp - the Krylov space solver context
78: - type - the type of refinement
80: Options Database Key:
81: . -ksp_orthogonalization_cgs_refinement_type (refine_never|refine_ifneeded|refine_always) - refinement type
83: Level: intermediate
85: Notes:
86: This option applies only if the orthogonalization method is classical Gram-Schmidt (CGS), see `KSPOrthogonalizationSet()`.
88: The default refinement type is `KSP_ORTHOGONALIZATION_CGS_REFINE_NEVER`.
90: For a very small set of problems, not using refinement, that is `KSP_ORTHOGONALIZATION_CGS_REFINE_NEVER`, may be unstable, thus causing `KSPSolve()`
91: to not converge.
93: .seealso: [](ch_ksp), `KSPOrthogonalizationSet()`, `KSPOrthogonalizationCGSRefinementType`, `KSPOrthogonalizationClassicalGramSchmidt()`, `KSPOrthogonalizationGetCGSRefinementType()`,
94: `KSPOrthogonalizationGet()`
95: @*/
96: PetscErrorCode KSPOrthogonalizationSetCGSRefinementType(KSP ksp, KSPOrthogonalizationCGSRefinementType type)
97: {
98: PetscFunctionBegin;
101: ksp->cgstype = type;
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*@
106: KSPOrthogonalizationGetCGSRefinementType - Gets the type of iterative refinement to use in the classical Gram-Schmidt
107: orthogonalization used by `KSPGMRES` and other solvers.
109: Not Collective
111: Input Parameter:
112: . ksp - the Krylov space solver context
114: Output Parameter:
115: . type - the type of refinement
117: Level: intermediate
119: .seealso: [](ch_ksp), `KSPOrthogonalizationSet()`, `KSPOrthogonalizationCGSRefinementType`, `KSPOrthogonalizationClassicalGramSchmidt()`, `KSPOrthogonalizationSetCGSRefinementType()`,
120: `KSPOrthogonalizationGet()`
121: @*/
122: PetscErrorCode KSPOrthogonalizationGetCGSRefinementType(KSP ksp, KSPOrthogonalizationCGSRefinementType *type)
123: {
124: PetscFunctionBegin;
126: PetscAssertPointer(type, 2);
127: *type = ksp->cgstype;
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: KSPOrthogonalizationModifiedGramSchmidt - This is the basic orthogonalization routine
133: using modified Gram-Schmidt.
135: Collective
137: Input Parameters:
138: + ksp - the Krylov space solver context
139: . V - array of previously computed orthonormal vectors
140: . n - number of vectors
141: - x - vector to be orthogonalized, modified on output (may be `NULL`)
143: Output Parameter:
144: . h - computed orthogonalization coefficients
146: Options Database Key:
147: . -ksp_orthogonalization mgs - choose modified Gram-Schmidt (MGS) for orthogonalization
149: Level: intermediate
151: Notes:
152: If no `x` is given, then the vector to be orthogonalized is assumed to be located at `V[n]`.
153: The input vectors `V` must be orthogonal and with unit two-norm. The output vector `x` is
154: not normalized.
156: In general this is much slower than `KSPOrthogonalizationClassicalGramSchmidt()` but has better stability properties.
158: .seealso: [](ch_ksp), `KSPOrthogonalizationSet()`, `KSPOrthogonalizationClassicalGramSchmidt()`, `KSPOrthogonalizationGet()`
159: @*/
160: PetscErrorCode KSPOrthogonalizationModifiedGramSchmidt(KSP ksp, Vec V[], PetscInt n, Vec x, PetscScalar h[])
161: {
162: PetscInt j;
163: PetscScalar *hh = h;
164: Vec z = x;
166: PetscFunctionBegin;
168: PetscAssertPointer(V, 2);
171: PetscAssertPointer(h, 5);
172: PetscCall(PetscLogEventBegin(KSP_Orthogonalization, ksp, 0, 0, 0));
173: if (!z) z = V[n];
174: for (j = 0; j < n; j++) {
175: /* (z, v(j)) */
176: PetscCall(VecDot(z, V[j], hh));
177: KSPCheckDot(ksp, *hh);
178: if (ksp->reason) break;
179: /* z <- z - hh[j] v(j) */
180: PetscCall(VecAXPY(z, -(*hh++), V[j]));
181: }
182: PetscCall(PetscLogEventEnd(KSP_Orthogonalization, ksp, 0, 0, 0));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: KSPOrthogonalizationClassicalGramSchmidt - This is the basic orthogonalization routine
188: using classical Gram-Schmidt with possible iterative refinement to improve the stability.
190: Collective
192: Input Parameters:
193: + ksp - the Krylov space solver context
194: . V - array of previously computed orthonormal vectors
195: . n - number of vectors
196: - x - vector to be orthogonalized, modified on output (may be `NULL`)
198: Output Parameter:
199: . h - computed orthogonalization coefficients
201: Options Database Keys:
202: + -ksp_orthogonalization cgs - choose classical Gram-Schmidt (CGS) for orthogonalization
203: - -ksp_orthogonalization_cgs_refinement_type (refine_never|refine_ifneeded|refine_always) - determine if iterative refinement is used to increase the stability of the
204: classical Gram-Schmidt orthogonalization
206: Level: intermediate
208: Notes:
209: If no `x` is given, then the vector to be orthogonalized is assumed to be located at `V[n]`.
210: The input vectors `V` must be orthogonal and with unit two-norm. The output vector `x` is
211: not normalized.
213: Use `KSPOrthogonalizationSetCGSRefinementType()` to determine if iterative refinement is to be used.
214: This is much faster than `KSPOrthogonalizationModifiedGramSchmidt()` but has the small possibility of stability issues
215: that can usually be handled by using a single step of iterative refinement with `KSPOrthogonalizationSetCGSRefinementType()`.
217: .seealso: [](ch_ksp), `KSPOrthogonalizationCGSRefinementType`, `KSPOrthogonalizationSet()`, `KSPOrthogonalizationSetCGSRefinementType()`,
218: `KSPOrthogonalizationGetCGSRefinementType()`, `KSPOrthogonalizationGet()`, `KSPOrthogonalizationModifiedGramSchmidt()`
219: @*/
220: PetscErrorCode KSPOrthogonalizationClassicalGramSchmidt(KSP ksp, Vec V[], PetscInt n, Vec x, PetscScalar h[])
221: {
222: PetscInt j;
223: PetscScalar *hh = h, *lhh;
224: Vec z = x;
225: PetscReal hnrm, wnrm;
226: PetscBool refine = (PetscBool)(ksp->cgstype == KSP_ORTHOGONALIZATION_CGS_REFINE_ALWAYS);
228: PetscFunctionBegin;
230: PetscAssertPointer(V, 2);
233: PetscAssertPointer(h, 5);
234: PetscCall(PetscLogEventBegin(KSP_Orthogonalization, ksp, 0, 0, 0));
235: if (!z) z = V[n];
236: if (ksp->lorthogwork < n) {
237: PetscCall(PetscFree(ksp->orthogwork));
238: ksp->lorthogwork = PetscMax(30, PetscMax(2 * ksp->lorthogwork, n));
239: PetscCall(PetscMalloc1(ksp->lorthogwork, &ksp->orthogwork));
240: }
241: lhh = ksp->orthogwork;
243: /* Clear hh since we will accumulate values into them */
244: for (j = 0; j < n; j++) hh[j] = 0.0;
246: /*
247: This is really a matrix-vector product, with the matrix stored
248: as pointer to rows
249: */
250: PetscCall(VecMDot(z, n, V, lhh)); /* <v,z> */
251: for (j = 0; j < n; j++) {
252: KSPCheckDot(ksp, lhh[j]);
253: if (ksp->reason) goto done;
254: lhh[j] = -lhh[j];
255: }
257: /*
258: This is really a matrix-vector product:
259: [h[0],h[1],...]*[ v[0]; v[1]; ...] subtracted from z.
260: */
261: PetscCall(VecMAXPY(z, n, lhh, V));
262: /* note lhh[j] is -<v,z> , hence the subtraction */
263: for (j = 0; j < n; j++) {
264: hh[j] -= lhh[j]; /* hh += <v,z> */
265: }
267: /*
268: the second step classical Gram-Schmidt is only necessary
269: when a simple test criteria is not passed
270: */
271: if (ksp->cgstype == KSP_ORTHOGONALIZATION_CGS_REFINE_IFNEEDED) {
272: hnrm = 0.0;
273: for (j = 0; j < n; j++) hnrm += PetscRealPart(lhh[j] * PetscConj(lhh[j]));
275: hnrm = PetscSqrtReal(hnrm);
276: PetscCall(VecNorm(z, NORM_2, &wnrm));
277: KSPCheckNorm(ksp, wnrm);
278: if (ksp->reason) goto done;
279: if (wnrm < hnrm) {
280: refine = PETSC_TRUE;
281: PetscCall(PetscInfo(ksp, "Performing iterative refinement wnorm %g hnorm %g\n", (double)wnrm, (double)hnrm));
282: }
283: }
285: if (refine) {
286: PetscCall(VecMDot(z, n, V, lhh)); /* <v,z> */
287: for (j = 0; j < n; j++) {
288: KSPCheckDot(ksp, lhh[j]);
289: if (ksp->reason) goto done;
290: lhh[j] = -lhh[j];
291: }
292: PetscCall(VecMAXPY(z, n, lhh, V));
293: /* note lhh[j] is -<v,z> , hence the subtraction */
294: for (j = 0; j < n; j++) {
295: hh[j] -= lhh[j]; /* hh += <v,z> */
296: }
297: }
298: done:
299: PetscCall(PetscLogEventEnd(KSP_Orthogonalization, ksp, 0, 0, 0));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }