Actual source code: borthog2.c
1: /*
2: Routines used for the orthogonalization of the Hessenberg matrix.
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>
10: /*@C
11: KSPGMRESClassicalGramSchmidtOrthogonalization - This is the basic orthogonalization routine
12: using classical Gram-Schmidt with possible iterative refinement to improve the stability
14: Collective
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
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.
25: Level: intermediate
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 a single step of iterative refinement with `KSPGMRESSetCGSRefinementType()`
32: .seealso: [](ch_ksp), `KSPGMRESCGSRefinementType`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetCGSRefinementType()`,
33: `KSPGMRESGetCGSRefinementType()`, `KSPGMRESGetOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`
34: @*/
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);
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;
48: /* update Hessenberg matrix and do unmodified Gram-Schmidt */
49: hh = HH(0, it);
50: hes = HES(0, it);
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: }
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: }
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: }
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]));
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: }
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: }