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;
21: PetscBool diagonalscale;
23: PetscFunctionBegin;
24: /* not sure what residual norm it does use, should use for right preconditioning */
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);
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;
40: /* Compute initial preconditioned residual */
41: PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));
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;
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);
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 */
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));
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;
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;
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;
124: PetscCall(KSPUnwindPreconditioner(ksp, X, T));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*MC
129: KSPCGS - This code implements the CGS (Conjugate Gradient Squared) method {cite}`so:89`.
131: Level: beginner
133: Notes:
134: Does not require a symmetric matrix. Does not apply transpose of the matrix.
136: Supports left and right preconditioning, but not symmetric.
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.
142: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBCGS`, `KSPSetPCSide()`
143: M*/
144: PETSC_EXTERN PetscErrorCode KSPCreate_CGS(KSP ksp)
145: {
146: PetscFunctionBegin;
147: ksp->data = (void *)0;
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));
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: }