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