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