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