1 #include <petsc/private/kspimpl.h>
2
KSPSetUp_CR(KSP ksp)3 static PetscErrorCode KSPSetUp_CR(KSP ksp)
4 {
5 PetscFunctionBegin;
6 PetscCall(KSPSetWorkVecs(ksp, 6));
7 PetscFunctionReturn(PETSC_SUCCESS);
8 }
9
KSPSolve_CR(KSP ksp)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;
18
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];
28
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) */
44
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 }
67
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);
76
77 i = 0;
78 do {
79 PetscCall(KSP_PCApply(ksp, AP, Q)); /* Q <- B* AP */
80
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) */
89
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));
95
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 }
119
120 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
121 ksp->its++;
122 ksp->rnorm = dp;
123 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
124
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;
129
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 }
138
139 /*MC
140 KSPCR - This code implements the (preconditioned) conjugate residuals method {cite}`hs:52`
141
142 Level: beginner
143
144 Notes:
145 The operator and the preconditioner must be symmetric for this method.
146
147 The preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE.
148
149 Support only for left preconditioning.
150
151 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`
152 M*/
KSPCreate_CR(KSP ksp)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));
160
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 }
170