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