1 /*
2 Routines used for the orthogonalization of the Hessenberg matrix.
3
4 Note that for the complex numbers version, the VecDot() and
5 VecMDot() arguments within the code MUST remain in the order
6 given for correct computation of inner products.
7 */
8 #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>
9
10 /*@C
11 KSPGMRESClassicalGramSchmidtOrthogonalization - This is the basic orthogonalization routine
12 using classical Gram-Schmidt with possible iterative refinement to improve the stability
13
14 Collective, No Fortran Support
15
16 Input Parameters:
17 + ksp - `KSP` object, must be associated with `KSPGMRES`, `KSPFGMRES`, or `KSPLGMRES` Krylov method
18 - it - one less than the current GMRES restart iteration, i.e. the size of the Krylov space
19
20 Options Database Keys:
21 + -ksp_gmres_classicalgramschmidt - Activates `KSPGMRESClassicalGramSchmidtOrthogonalization()`
22 - -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is
23 used to increase the stability of the classical Gram-Schmidt orthogonalization.
24
25 Level: intermediate
26
27 Note:
28 Use `KSPGMRESSetCGSRefinementType()` to determine if iterative refinement is to be used.
29 This is much faster than `KSPGMRESModifiedGramSchmidtOrthogonalization()` but has the small possibility of stability issues
30 that can usually be handled by using a single step of iterative refinement with `KSPGMRESSetCGSRefinementType()`
31
32 .seealso: [](ch_ksp), `KSPGMRESCGSRefinementType`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetCGSRefinementType()`,
33 `KSPGMRESGetCGSRefinementType()`, `KSPGMRESGetOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`
34 @*/
KSPGMRESClassicalGramSchmidtOrthogonalization(KSP ksp,PetscInt it)35 PetscErrorCode KSPGMRESClassicalGramSchmidtOrthogonalization(KSP ksp, PetscInt it)
36 {
37 KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
38 PetscInt j;
39 PetscScalar *hh, *hes, *lhh;
40 PetscReal hnrm, wnrm;
41 PetscBool refine = (PetscBool)(gmres->cgstype == KSP_GMRES_CGS_REFINE_ALWAYS);
42
43 PetscFunctionBegin;
44 PetscCall(PetscLogEventBegin(KSP_GMRESOrthogonalization, ksp, 0, 0, 0));
45 if (!gmres->orthogwork) PetscCall(PetscMalloc1(gmres->max_k + 2, &gmres->orthogwork));
46 lhh = gmres->orthogwork;
47
48 /* update Hessenberg matrix and do unmodified Gram-Schmidt */
49 hh = HH(0, it);
50 hes = HES(0, it);
51
52 /* Clear hh and hes since we will accumulate values into them */
53 for (j = 0; j <= it; j++) {
54 hh[j] = 0.0;
55 hes[j] = 0.0;
56 }
57
58 /*
59 This is really a matrix-vector product, with the matrix stored
60 as pointer to rows
61 */
62 PetscCall(VecMDot(VEC_VV(it + 1), it + 1, &(VEC_VV(0)), lhh)); /* <v,vnew> */
63 for (j = 0; j <= it; j++) {
64 KSPCheckDot(ksp, lhh[j]);
65 if (ksp->reason) goto done;
66 lhh[j] = -lhh[j];
67 }
68
69 /*
70 This is really a matrix-vector product:
71 [h[0],h[1],...]*[ v[0]; v[1]; ...] subtracted from v[it+1].
72 */
73 PetscCall(VecMAXPY(VEC_VV(it + 1), it + 1, lhh, &VEC_VV(0)));
74 /* note lhh[j] is -<v,vnew> , hence the subtraction */
75 for (j = 0; j <= it; j++) {
76 hh[j] -= lhh[j]; /* hh += <v,vnew> */
77 hes[j] -= lhh[j]; /* hes += <v,vnew> */
78 }
79
80 /*
81 the second step classical Gram-Schmidt is only necessary
82 when a simple test criteria is not passed
83 */
84 if (gmres->cgstype == KSP_GMRES_CGS_REFINE_IFNEEDED) {
85 hnrm = 0.0;
86 for (j = 0; j <= it; j++) hnrm += PetscRealPart(lhh[j] * PetscConj(lhh[j]));
87
88 hnrm = PetscSqrtReal(hnrm);
89 PetscCall(VecNorm(VEC_VV(it + 1), NORM_2, &wnrm));
90 KSPCheckNorm(ksp, wnrm);
91 if (ksp->reason) goto done;
92 if (wnrm < hnrm) {
93 refine = PETSC_TRUE;
94 PetscCall(PetscInfo(ksp, "Performing iterative refinement wnorm %g hnorm %g\n", (double)wnrm, (double)hnrm));
95 }
96 }
97
98 if (refine) {
99 PetscCall(VecMDot(VEC_VV(it + 1), it + 1, &(VEC_VV(0)), lhh)); /* <v,vnew> */
100 for (j = 0; j <= it; j++) {
101 KSPCheckDot(ksp, lhh[j]);
102 if (ksp->reason) goto done;
103 lhh[j] = -lhh[j];
104 }
105 PetscCall(VecMAXPY(VEC_VV(it + 1), it + 1, lhh, &VEC_VV(0)));
106 /* note lhh[j] is -<v,vnew> , hence the subtraction */
107 for (j = 0; j <= it; j++) {
108 hh[j] -= lhh[j]; /* hh += <v,vnew> */
109 hes[j] -= lhh[j]; /* hes += <v,vnew> */
110 }
111 }
112 done:
113 PetscCall(PetscLogEventEnd(KSP_GMRESOrthogonalization, ksp, 0, 0, 0));
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116