Actual source code: cr.c
2: #include <petsc/private/kspimpl.h>
4: static PetscErrorCode KSPSetUp_CR(KSP ksp)
5: {
6: KSPSetWorkVecs(ksp, 6);
7: return 0;
8: }
10: static PetscErrorCode KSPSolve_CR(KSP ksp)
11: {
12: PetscInt i = 0;
13: PetscReal dp;
14: PetscScalar ai, bi;
15: PetscScalar apq, btop, bbot;
16: Vec X, B, R, RT, P, AP, ART, Q;
17: Mat Amat, Pmat;
19: X = ksp->vec_sol;
20: B = ksp->vec_rhs;
21: R = ksp->work[0];
22: RT = ksp->work[1];
23: P = ksp->work[2];
24: AP = ksp->work[3];
25: ART = ksp->work[4];
26: Q = ksp->work[5];
28: /* R is the true residual norm, RT is the preconditioned residual norm */
29: PCGetOperators(ksp->pc, &Amat, &Pmat);
30: if (!ksp->guess_zero) {
31: KSP_MatMult(ksp, Amat, X, R); /* R <- A*X */
32: VecAYPX(R, -1.0, B); /* R <- B-R == B-A*X */
33: } else {
34: VecCopy(B, R); /* R <- B (X is 0) */
35: }
36: /* This may be true only on a subset of MPI ranks; setting it here so it will be detected by the first norm computation below */
37: if (ksp->reason == KSP_DIVERGED_PC_FAILED) VecSetInf(R);
38: KSP_PCApply(ksp, R, P); /* P <- B*R */
39: KSP_MatMult(ksp, Amat, P, AP); /* AP <- A*P */
40: VecCopy(P, RT); /* RT <- P */
41: VecCopy(AP, ART); /* ART <- AP */
42: VecDotBegin(RT, ART, &btop); /* (RT,ART) */
44: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
45: VecNormBegin(RT, NORM_2, &dp); /* dp <- RT'*RT */
46: VecDotEnd(RT, ART, &btop); /* (RT,ART) */
47: VecNormEnd(RT, NORM_2, &dp); /* dp <- RT'*RT */
48: KSPCheckNorm(ksp, dp);
49: } else if (ksp->normtype == KSP_NORM_NONE) {
50: dp = 0.0; /* meaningless value that is passed to monitor and convergence test */
51: VecDotEnd(RT, ART, &btop); /* (RT,ART) */
52: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
53: VecNormBegin(R, NORM_2, &dp); /* dp <- R'*R */
54: VecDotEnd(RT, ART, &btop); /* (RT,ART) */
55: VecNormEnd(R, NORM_2, &dp); /* dp <- RT'*RT */
56: KSPCheckNorm(ksp, dp);
57: } else if (ksp->normtype == KSP_NORM_NATURAL) {
58: VecDotEnd(RT, ART, &btop); /* (RT,ART) */
59: dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */
60: } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSPNormType of %d not supported", (int)ksp->normtype);
61: if (PetscAbsScalar(btop) < 0.0) {
62: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
63: PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n");
64: return 0;
65: }
67: ksp->its = 0;
68: KSPMonitor(ksp, 0, dp);
69: PetscObjectSAWsTakeAccess((PetscObject)ksp);
70: ksp->rnorm = dp;
71: PetscObjectSAWsGrantAccess((PetscObject)ksp);
72: KSPLogResidualHistory(ksp, dp);
73: (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP);
74: if (ksp->reason) return 0;
76: i = 0;
77: do {
78: KSP_PCApply(ksp, AP, Q); /* Q <- B* AP */
80: VecDot(AP, Q, &apq);
81: KSPCheckDot(ksp, apq);
82: if (PetscRealPart(apq) <= 0.0) {
83: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
84: PetscInfo(ksp, "KSPSolve_CR:diverging due to indefinite or negative definite PC\n");
85: break;
86: }
87: ai = btop / apq; /* ai = (RT,ART)/(AP,Q) */
89: VecAXPY(X, ai, P); /* X <- X + ai*P */
90: VecAXPY(RT, -ai, Q); /* RT <- RT - ai*Q */
91: KSP_MatMult(ksp, Amat, RT, ART); /* ART <- A*RT */
92: bbot = btop;
93: VecDotBegin(RT, ART, &btop);
95: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
96: VecNormBegin(RT, NORM_2, &dp); /* dp <- || RT || */
97: VecDotEnd(RT, ART, &btop);
98: VecNormEnd(RT, NORM_2, &dp); /* dp <- || RT || */
99: KSPCheckNorm(ksp, dp);
100: } else if (ksp->normtype == KSP_NORM_NATURAL) {
101: VecDotEnd(RT, ART, &btop);
102: dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */
103: } else if (ksp->normtype == KSP_NORM_NONE) {
104: VecDotEnd(RT, ART, &btop);
105: dp = 0.0; /* meaningless value that is passed to monitor and convergence test */
106: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
107: VecAXPY(R, ai, AP); /* R <- R - ai*AP */
108: VecNormBegin(R, NORM_2, &dp); /* dp <- R'*R */
109: VecDotEnd(RT, ART, &btop);
110: VecNormEnd(R, NORM_2, &dp); /* dp <- R'*R */
111: KSPCheckNorm(ksp, dp);
112: } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSPNormType of %d not supported", (int)ksp->normtype);
113: if (PetscAbsScalar(btop) < 0.0) {
114: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
115: PetscInfo(ksp, "diverging due to indefinite or negative definite PC\n");
116: break;
117: }
119: PetscObjectSAWsTakeAccess((PetscObject)ksp);
120: ksp->its++;
121: ksp->rnorm = dp;
122: PetscObjectSAWsGrantAccess((PetscObject)ksp);
124: KSPLogResidualHistory(ksp, dp);
125: KSPMonitor(ksp, i + 1, dp);
126: (*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP);
127: if (ksp->reason) break;
129: bi = btop / bbot;
130: VecAYPX(P, bi, RT); /* P <- RT + Bi P */
131: VecAYPX(AP, bi, ART); /* AP <- ART + Bi AP */
132: i++;
133: } while (i < ksp->max_it);
134: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
135: return 0;
136: }
138: /*MC
139: KSPCR - This code implements the (preconditioned) conjugate residuals method
141: Level: beginner
143: Notes:
144: The operator and the preconditioner must be symmetric for this method.
146: The preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE.
148: Support only for left preconditioning.
150: References:
151: . * - Magnus R. Hestenes and Eduard Stiefel, Methods of Conjugate Gradients for Solving Linear Systems,
152: Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379
154: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`
155: M*/
156: PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp)
157: {
158: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
159: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
160: KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
161: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
163: ksp->ops->setup = KSPSetUp_CR;
164: ksp->ops->solve = KSPSolve_CR;
165: ksp->ops->destroy = KSPDestroyDefault;
166: ksp->ops->buildsolution = KSPBuildSolutionDefault;
167: ksp->ops->buildresidual = KSPBuildResidualDefault;
168: ksp->ops->setfromoptions = NULL;
169: ksp->ops->view = NULL;
170: return 0;
171: }