1 /*
2 Note that for the complex numbers version, the VecDot() arguments
3 within the code MUST remain in the order given for correct computation
4 of inner products.
5 */
6 #include <petsc/private/kspimpl.h>
7
KSPSetUp_CGS(KSP ksp)8 static PetscErrorCode KSPSetUp_CGS(KSP ksp)
9 {
10 PetscFunctionBegin;
11 PetscCall(KSPSetWorkVecs(ksp, 7));
12 PetscFunctionReturn(PETSC_SUCCESS);
13 }
14
KSPSolve_CGS(KSP ksp)15 static PetscErrorCode KSPSolve_CGS(KSP ksp)
16 {
17 PetscInt i;
18 PetscScalar rho, rhoold, a, s, b;
19 Vec X, B, V, P, R, RP, T, Q, U, AUQ;
20 PetscReal dp = 0.0;
21 PetscBool diagonalscale;
22
23 PetscFunctionBegin;
24 /* not sure what residual norm it does use, should use for right preconditioning */
25
26 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
27 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
28
29 X = ksp->vec_sol;
30 B = ksp->vec_rhs;
31 R = ksp->work[0];
32 RP = ksp->work[1];
33 V = ksp->work[2];
34 T = ksp->work[3];
35 Q = ksp->work[4];
36 P = ksp->work[5];
37 U = ksp->work[6];
38 AUQ = V;
39
40 /* Compute initial preconditioned residual */
41 PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));
42
43 /* Test for nothing to do */
44 if (ksp->normtype != KSP_NORM_NONE) {
45 PetscCall(VecNorm(R, NORM_2, &dp));
46 KSPCheckNorm(ksp, dp);
47 if (ksp->normtype == KSP_NORM_NATURAL) dp *= dp;
48 } else dp = 0.0;
49
50 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
51 ksp->its = 0;
52 ksp->rnorm = dp;
53 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
54 PetscCall(KSPLogResidualHistory(ksp, dp));
55 PetscCall(KSPMonitor(ksp, 0, dp));
56 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
57 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
58
59 /* Make the initial Rp == R */
60 PetscCall(VecCopy(R, RP));
61 /* added for Fidap */
62 /* Penalize Startup - Isaac Hasbani Trick for CGS
63 Since most initial conditions result in a mostly 0 residual,
64 we change all the 0 values in the vector RP to the maximum.
65 */
66 if (ksp->normtype == KSP_NORM_NATURAL) {
67 PetscReal vr0max;
68 PetscScalar *tmp_RP = NULL;
69 PetscInt numnp = 0, *max_pos = NULL;
70 PetscCall(VecMax(RP, max_pos, &vr0max));
71 PetscCall(VecGetArray(RP, &tmp_RP));
72 PetscCall(VecGetLocalSize(RP, &numnp));
73 for (i = 0; i < numnp; i++) {
74 if (tmp_RP[i] == 0.0) tmp_RP[i] = vr0max;
75 }
76 PetscCall(VecRestoreArray(RP, &tmp_RP));
77 }
78 /* end of addition for Fidap */
79
80 /* Set the initial conditions */
81 PetscCall(VecDot(R, RP, &rhoold)); /* rhoold = (r,rp) */
82 PetscCall(VecCopy(R, U));
83 PetscCall(VecCopy(R, P));
84 PetscCall(KSP_PCApplyBAorAB(ksp, P, V, T));
85
86 i = 0;
87 do {
88 PetscCall(VecDot(V, RP, &s)); /* s <- (v,rp) */
89 KSPCheckDot(ksp, s);
90 a = rhoold / s; /* a <- rho / s */
91 PetscCall(VecWAXPY(Q, -a, V, U)); /* q <- u - a v */
92 PetscCall(VecWAXPY(T, 1.0, U, Q)); /* t <- u + q */
93 PetscCall(VecAXPY(X, a, T)); /* x <- x + a (u + q) */
94 PetscCall(KSP_PCApplyBAorAB(ksp, T, AUQ, U));
95 PetscCall(VecAXPY(R, -a, AUQ)); /* r <- r - a K (u + q) */
96 PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp) */
97 KSPCheckDot(ksp, rho);
98 if (ksp->normtype == KSP_NORM_NATURAL) {
99 dp = PetscAbsScalar(rho);
100 } else if (ksp->normtype != KSP_NORM_NONE) {
101 PetscCall(VecNorm(R, NORM_2, &dp));
102 KSPCheckNorm(ksp, dp);
103 } else dp = 0.0;
104
105 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
106 ksp->its++;
107 ksp->rnorm = dp;
108 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
109 PetscCall(KSPLogResidualHistory(ksp, dp));
110 PetscCall(KSPMonitor(ksp, i + 1, dp));
111 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
112 if (ksp->reason) break;
113
114 b = rho / rhoold; /* b <- rho / rhoold */
115 PetscCall(VecWAXPY(U, b, Q, R)); /* u <- r + b q */
116 PetscCall(VecAXPY(Q, b, P));
117 PetscCall(VecWAXPY(P, b, Q, U)); /* p <- u + b(q + b p) */
118 PetscCall(KSP_PCApplyBAorAB(ksp, P, V, Q)); /* v <- K p */
119 rhoold = rho;
120 i++;
121 } while (i < ksp->max_it);
122 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
123
124 PetscCall(KSPUnwindPreconditioner(ksp, X, T));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 /*MC
129 KSPCGS - This code implements the CGS (Conjugate Gradient Squared) method {cite}`so:89`.
130
131 Level: beginner
132
133 Notes:
134 Does not require a symmetric matrix. Does not apply transpose of the matrix.
135
136 Supports left and right preconditioning, but not symmetric.
137
138 Developer Note:
139 Has this weird support for doing the convergence test with the natural norm, I assume this works only with
140 no preconditioning and symmetric positive definite operator.
141
142 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBCGS`, `KSPSetPCSide()`
143 M*/
KSPCreate_CGS(KSP ksp)144 PETSC_EXTERN PetscErrorCode KSPCreate_CGS(KSP ksp)
145 {
146 PetscFunctionBegin;
147 ksp->data = NULL;
148
149 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
150 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
151 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
152 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_RIGHT, 2));
153 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
154 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
155
156 ksp->ops->setup = KSPSetUp_CGS;
157 ksp->ops->solve = KSPSolve_CGS;
158 ksp->ops->destroy = KSPDestroyDefault;
159 ksp->ops->buildsolution = KSPBuildSolutionDefault;
160 ksp->ops->buildresidual = KSPBuildResidualDefault;
161 ksp->ops->setfromoptions = NULL;
162 ksp->ops->view = NULL;
163 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165