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