1 /*
2 This file implements the conjugate gradient method in PETSc as part of
3 KSP. You can use this as a starting point for implementing your own
4 Krylov method that is not provided with PETSc.
5
6 The following basic routines are required for each Krylov method.
7 KSPCreate_XXX() - Creates the Krylov context
8 KSPSetFromOptions_XXX() - Sets runtime options
9 KSPSolve_XXX() - Runs the Krylov method
10 KSPDestroy_XXX() - Destroys the Krylov context, freeing all
11 memory it needed
12 Here the "_XXX" denotes a particular implementation, in this case
13 we use _CG (e.g. KSPCreate_CG, KSPDestroy_CG). These routines
14 are actually called via the common user interface routines
15 KSPSetType(), KSPSetFromOptions(), KSPSolve(), and KSPDestroy() so the
16 application code interface remains identical for all preconditioners.
17
18 Other basic routines for the KSP objects include
19 KSPSetUp_XXX()
20 KSPView_XXX() - Prints details of solver being used.
21
22 Detailed Notes:
23 By default, this code implements the CG (Conjugate Gradient) method,
24 which is valid for real symmetric (and complex Hermitian) positive
25 definite matrices. Note that for the complex Hermitian case, the
26 VecDot() arguments within the code MUST remain in the order given
27 for correct computation of inner products.
28
29 Reference: Hestenes and Steifel, 1952.
30
31 By switching to the indefinite vector inner product, VecTDot(), the
32 same code is used for the complex symmetric case as well. The user
33 must call KSPCGSetType(ksp,KSP_CG_SYMMETRIC) or use the option
34 -ksp_cg_type symmetric to invoke this variant for the complex case.
35 Note, however, that the complex symmetric code is NOT valid for
36 all such matrices ... and thus we don't recommend using this method.
37 */
38 /*
39 cgimpl.h defines the simple data structured used to store information
40 related to the type of matrix (e.g. complex symmetric) being solved and
41 data used during the optional Lanczos process used to compute eigenvalues
42 */
43 #include <../src/ksp/ksp/impls/cg/cgimpl.h> /*I "petscksp.h" I*/
44 extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP, PetscReal *, PetscReal *);
45 extern PetscErrorCode KSPComputeEigenvalues_CG(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
46
KSPCGSetObjectiveTarget_CG(KSP ksp,PetscReal obj_min)47 static PetscErrorCode KSPCGSetObjectiveTarget_CG(KSP ksp, PetscReal obj_min)
48 {
49 KSP_CG *cg = (KSP_CG *)ksp->data;
50
51 PetscFunctionBegin;
52 cg->obj_min = obj_min;
53 PetscFunctionReturn(PETSC_SUCCESS);
54 }
55
KSPCGSetRadius_CG(KSP ksp,PetscReal radius)56 static PetscErrorCode KSPCGSetRadius_CG(KSP ksp, PetscReal radius)
57 {
58 KSP_CG *cg = (KSP_CG *)ksp->data;
59
60 PetscFunctionBegin;
61 cg->radius = radius;
62 PetscFunctionReturn(PETSC_SUCCESS);
63 }
64
KSPCGGetObjFcn_CG(KSP ksp,PetscReal * obj)65 static PetscErrorCode KSPCGGetObjFcn_CG(KSP ksp, PetscReal *obj)
66 {
67 KSP_CG *cg = (KSP_CG *)ksp->data;
68
69 PetscFunctionBegin;
70 *obj = cg->obj;
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
74 /*
75 KSPSetUp_CG - Sets up the workspace needed by the CG method.
76
77 This is called once, usually automatically by KSPSolve() or KSPSetUp()
78 but can be called directly by KSPSetUp()
79 */
KSPSetUp_CG(KSP ksp)80 static PetscErrorCode KSPSetUp_CG(KSP ksp)
81 {
82 KSP_CG *cgP = (KSP_CG *)ksp->data;
83 PetscInt maxit = ksp->max_it, nwork = 3;
84
85 PetscFunctionBegin;
86 /* get work vectors needed by CG */
87 if (cgP->singlereduction) nwork += 2;
88 PetscCall(KSPSetWorkVecs(ksp, nwork));
89
90 /*
91 If user requested computations of eigenvalues then allocate
92 work space needed
93 */
94 if (ksp->calc_sings) {
95 PetscCall(PetscFree4(cgP->e, cgP->d, cgP->ee, cgP->dd));
96 PetscCall(PetscMalloc4(maxit + 1, &cgP->e, maxit, &cgP->d, maxit, &cgP->ee, maxit, &cgP->dd));
97
98 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
99 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
100 }
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
104 /*
105 A macro used in the following KSPSolve_CG and KSPSolve_CG_SingleReduction routines
106 */
107 #define VecXDot(x, y, a) (cg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))
108
109 /*
110 KSPSolve_CG - This routine actually applies the conjugate gradient method
111
112 Note : this routine can be replaced with another one (see below) which implements
113 another variant of CG.
114
115 Input Parameter:
116 . ksp - the Krylov space object that was set to use conjugate gradient, by, for
117 example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
118 */
KSPSolve_CG(KSP ksp)119 static PetscErrorCode KSPSolve_CG(KSP ksp)
120 {
121 PetscInt i, stored_max_it, eigs;
122 PetscScalar dpi = 0.0, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL, dpiold;
123 PetscReal dp = 0.0;
124 PetscReal r2, norm_p, norm_d, dMp;
125 Vec X, B, Z, R, P, W;
126 KSP_CG *cg;
127 Mat Amat, Pmat;
128 PetscBool diagonalscale, testobj;
129
130 PetscFunctionBegin;
131 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
132 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
133
134 cg = (KSP_CG *)ksp->data;
135 eigs = ksp->calc_sings;
136 stored_max_it = ksp->max_it;
137 X = ksp->vec_sol;
138 B = ksp->vec_rhs;
139 R = ksp->work[0];
140 Z = ksp->work[1];
141 P = ksp->work[2];
142 W = Z;
143 r2 = PetscSqr(cg->radius);
144
145 if (eigs) {
146 e = cg->e;
147 d = cg->d;
148 e[0] = 0.0;
149 }
150 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
151
152 ksp->its = 0;
153 if (!ksp->guess_zero) {
154 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
155
156 PetscCall(VecAYPX(R, -1.0, B));
157 if (cg->radius) { /* XXX direction? */
158 PetscCall(VecNorm(X, NORM_2, &norm_d));
159 norm_d *= norm_d;
160 }
161 } else {
162 PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
163 norm_d = 0.0;
164 }
165 /* 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 */
166 PetscCall(VecFlag(R, ksp->reason == KSP_DIVERGED_PC_FAILED));
167
168 switch (ksp->normtype) {
169 case KSP_NORM_PRECONDITIONED:
170 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
171 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z = e'*A'*B'*B*A*e */
172 KSPCheckNorm(ksp, dp);
173 break;
174 case KSP_NORM_UNPRECONDITIONED:
175 PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r = e'*A'*A*e */
176 KSPCheckNorm(ksp, dp);
177 break;
178 case KSP_NORM_NATURAL:
179 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
180 PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
181 KSPCheckDot(ksp, beta);
182 dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
183 break;
184 case KSP_NORM_NONE:
185 dp = 0.0;
186 break;
187 default:
188 SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
189 }
190
191 /* Initialize objective function
192 obj = 1/2 x^T A x - x^T b */
193 testobj = (PetscBool)(cg->obj_min < 0.0);
194 PetscCall(VecXDot(R, X, &a));
195 cg->obj = 0.5 * PetscRealPart(a);
196 PetscCall(VecXDot(B, X, &a));
197 cg->obj -= 0.5 * PetscRealPart(a);
198
199 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", ksp->its, (double)cg->obj));
200 PetscCall(KSPLogResidualHistory(ksp, dp));
201 PetscCall(KSPMonitor(ksp, ksp->its, dp));
202 ksp->rnorm = dp;
203
204 PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
205
206 if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
207 PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
208 ksp->reason = KSP_CONVERGED_ATOL;
209 }
210
211 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
212
213 if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
214 if (ksp->normtype != KSP_NORM_NATURAL) {
215 PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
216 KSPCheckDot(ksp, beta);
217 }
218
219 i = 0;
220 do {
221 ksp->its = i + 1;
222 if (beta == 0.0) {
223 ksp->reason = KSP_CONVERGED_ATOL;
224 PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
225 break;
226 #if !defined(PETSC_USE_COMPLEX)
227 } else if ((i > 0) && (beta * betaold < 0.0)) {
228 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner, beta %g, betaold %g", (double)PetscRealPart(beta), (double)PetscRealPart(betaold));
229 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
230 PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
231 break;
232 #endif
233 }
234 if (!i) {
235 PetscCall(VecCopy(Z, P)); /* p <- z */
236 if (cg->radius) {
237 PetscCall(VecNorm(P, NORM_2, &norm_p));
238 norm_p *= norm_p;
239 dMp = 0.0;
240 if (!ksp->guess_zero) PetscCall(VecDotRealPart(X, P, &dMp));
241 }
242 b = 0.0;
243 } else {
244 b = beta / betaold;
245 if (eigs) {
246 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
247 e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
248 }
249 PetscCall(VecAYPX(P, b, Z)); /* p <- z + b* p */
250 if (cg->radius) {
251 PetscCall(VecDotRealPart(X, P, &dMp));
252 PetscCall(VecNorm(P, NORM_2, &norm_p));
253 norm_p *= norm_p;
254 }
255 }
256 dpiold = dpi;
257 PetscCall(KSP_MatMult(ksp, Amat, P, W)); /* w <- Ap */
258 PetscCall(VecXDot(P, W, &dpi)); /* dpi <- p'w */
259 KSPCheckDot(ksp, dpi);
260 betaold = beta;
261
262 if ((dpi == 0.0) || ((i > 0) && ((PetscSign(PetscRealPart(dpi)) * PetscSign(PetscRealPart(dpiold))) < 0.0))) {
263 if (cg->radius) {
264 a = 0.0;
265 if (i == 0) {
266 if (norm_p > 0.0) {
267 a = PetscSqrtReal(r2 / norm_p);
268 } else {
269 PetscCall(VecNorm(R, NORM_2, &dp));
270 a = cg->radius > dp ? 1.0 : cg->radius / dp;
271 }
272 } else if (norm_p > 0.0) {
273 a = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
274 }
275 PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
276 cg->obj += PetscRealPart(a * (0.5 * a * dpi - betaold));
277 }
278 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " N obj %g\n", i + 1, (double)cg->obj));
279 if (ksp->converged_neg_curve) {
280 PetscCall(PetscInfo(ksp, "converged due to negative curvature: %g\n", (double)(PetscRealPart(dpi))));
281 ksp->reason = KSP_CONVERGED_NEG_CURVE;
282 } else {
283 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix, dpi %g, dpiold %g", (double)PetscRealPart(dpi), (double)PetscRealPart(dpiold));
284 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
285 PetscCall(PetscInfo(ksp, "diverging due to indefinite matrix\n"));
286 }
287 break;
288 }
289 a = beta / dpi; /* a = beta/p'w */
290 if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
291 if (cg->radius) { /* Steihaugh-Toint */
292 PetscReal norm_dp1 = norm_d + PetscRealPart(a) * (2.0 * dMp + PetscRealPart(a) * norm_p);
293 if (norm_dp1 > r2) {
294 ksp->reason = KSP_CONVERGED_STEP_LENGTH;
295 PetscCall(PetscInfo(ksp, "converged to the trust region radius %g\n", (double)cg->radius));
296 if (norm_p > 0.0) {
297 dp = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
298 PetscCall(VecAXPY(X, dp, P)); /* x <- x + ap */
299 cg->obj += PetscRealPart(dp * (0.5 * dp * dpi - beta));
300 }
301 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " R obj %g\n", i + 1, (double)cg->obj));
302 break;
303 }
304 }
305 PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
306 PetscCall(VecAXPY(R, -a, W)); /* r <- r - aw */
307 if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
308 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
309 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z */
310 KSPCheckNorm(ksp, dp);
311 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
312 PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r */
313 KSPCheckNorm(ksp, dp);
314 } else if (ksp->normtype == KSP_NORM_NATURAL) {
315 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
316 PetscCall(VecXDot(Z, R, &beta)); /* beta <- r'*z */
317 KSPCheckDot(ksp, beta);
318 dp = PetscSqrtReal(PetscAbsScalar(beta));
319 } else {
320 dp = 0.0;
321 }
322 cg->obj -= PetscRealPart(0.5 * a * betaold);
323 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", i + 1, (double)cg->obj));
324
325 ksp->rnorm = dp;
326 PetscCall(KSPLogResidualHistory(ksp, dp));
327 PetscCall(KSPMonitor(ksp, i + 1, dp));
328 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
329
330 if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
331 PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
332 ksp->reason = KSP_CONVERGED_ATOL;
333 }
334
335 if (ksp->reason) break;
336
337 if (cg->radius) {
338 PetscCall(VecNorm(X, NORM_2, &norm_d));
339 norm_d *= norm_d;
340 }
341
342 if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
343 if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
344 PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
345 KSPCheckDot(ksp, beta);
346 }
347
348 i++;
349 } while (i < ksp->max_it);
350 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
351 PetscFunctionReturn(PETSC_SUCCESS);
352 }
353
354 /*
355 KSPSolve_CG_SingleReduction
356
357 This variant of CG is identical in exact arithmetic to the standard algorithm,
358 but is rearranged to use only a single reduction stage per iteration, using additional
359 intermediate vectors.
360
361 See KSPCGUseSingleReduction_CG()
362
363 */
KSPSolve_CG_SingleReduction(KSP ksp)364 static PetscErrorCode KSPSolve_CG_SingleReduction(KSP ksp)
365 {
366 PetscInt i, stored_max_it, eigs;
367 PetscScalar dpi = 0.0, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL, delta, dpiold, tmp[2];
368 PetscReal dp = 0.0;
369 Vec X, B, Z, R, P, S, W, tmpvecs[2];
370 KSP_CG *cg;
371 Mat Amat, Pmat;
372 PetscBool diagonalscale;
373
374 PetscFunctionBegin;
375 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
376 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
377
378 cg = (KSP_CG *)ksp->data;
379 eigs = ksp->calc_sings;
380 stored_max_it = ksp->max_it;
381 X = ksp->vec_sol;
382 B = ksp->vec_rhs;
383 R = ksp->work[0];
384 Z = ksp->work[1];
385 P = ksp->work[2];
386 S = ksp->work[3];
387 W = ksp->work[4];
388
389 if (eigs) {
390 e = cg->e;
391 d = cg->d;
392 e[0] = 0.0;
393 }
394 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
395
396 ksp->its = 0;
397 if (!ksp->guess_zero) {
398 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
399 PetscCall(VecAYPX(R, -1.0, B));
400 } else {
401 PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
402 }
403
404 switch (ksp->normtype) {
405 case KSP_NORM_PRECONDITIONED:
406 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
407 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z = e'*A'*B'*B*A'*e' */
408 KSPCheckNorm(ksp, dp);
409 break;
410 case KSP_NORM_UNPRECONDITIONED:
411 PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r = e'*A'*A*e */
412 KSPCheckNorm(ksp, dp);
413 break;
414 case KSP_NORM_NATURAL:
415 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
416 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
417 PetscCall(VecXDot(Z, S, &delta)); /* delta <- z'*A*z = r'*B*A*B*r */
418 PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
419 KSPCheckDot(ksp, beta);
420 dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
421 break;
422 case KSP_NORM_NONE:
423 dp = 0.0;
424 break;
425 default:
426 SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
427 }
428 PetscCall(KSPLogResidualHistory(ksp, dp));
429 PetscCall(KSPMonitor(ksp, 0, dp));
430 ksp->rnorm = dp;
431
432 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
433 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
434
435 if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
436 if (ksp->normtype != KSP_NORM_NATURAL) {
437 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
438 PetscCall(VecXDot(Z, S, &delta)); /* delta <- z'*A*z = r'*B*A*B*r */
439 PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
440 KSPCheckDot(ksp, beta);
441 }
442
443 i = 0;
444 do {
445 ksp->its = i + 1;
446 if (beta == 0.0) {
447 ksp->reason = KSP_CONVERGED_ATOL;
448 PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
449 break;
450 #if !defined(PETSC_USE_COMPLEX)
451 } else if ((i > 0) && (beta * betaold < 0.0)) {
452 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner");
453 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
454 PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
455 break;
456 #endif
457 }
458 if (!i) {
459 PetscCall(VecCopy(Z, P)); /* p <- z */
460 b = 0.0;
461 } else {
462 b = beta / betaold;
463 if (eigs) {
464 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
465 e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
466 }
467 PetscCall(VecAYPX(P, b, Z)); /* p <- z + b* p */
468 }
469 dpiold = dpi;
470 if (!i) {
471 PetscCall(KSP_MatMult(ksp, Amat, P, W)); /* w <- Ap */
472 PetscCall(VecXDot(P, W, &dpi)); /* dpi <- p'w */
473 } else {
474 PetscCall(VecAYPX(W, beta / betaold, S)); /* w <- Ap */
475 dpi = delta - beta * beta * dpiold / (betaold * betaold); /* dpi <- p'w */
476 }
477 betaold = beta;
478 KSPCheckDot(ksp, beta);
479
480 if ((dpi == 0.0) || ((i > 0) && (PetscRealPart(dpi * dpiold) <= 0.0))) {
481 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix");
482 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
483 PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n"));
484 break;
485 }
486 a = beta / dpi; /* a = beta/p'w */
487 if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
488 PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
489 PetscCall(VecAXPY(R, -a, W)); /* r <- r - aw */
490 if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
491 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
492 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
493 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z */
494 KSPCheckNorm(ksp, dp);
495 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
496 PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r */
497 KSPCheckNorm(ksp, dp);
498 } else if (ksp->normtype == KSP_NORM_NATURAL) {
499 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
500 tmpvecs[0] = S;
501 tmpvecs[1] = R;
502 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
503 PetscCall(VecMDot(Z, 2, tmpvecs, tmp)); /* delta <- z'*A*z = r'*B*A*B*r */
504 delta = tmp[0];
505 beta = tmp[1]; /* beta <- z'*r */
506 KSPCheckDot(ksp, beta);
507 dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
508 } else {
509 dp = 0.0;
510 }
511 ksp->rnorm = dp;
512 PetscCall(KSPLogResidualHistory(ksp, dp));
513 PetscCall(KSPMonitor(ksp, i + 1, dp));
514 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
515 if (ksp->reason) break;
516
517 if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) {
518 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
519 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
520 }
521 if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
522 tmpvecs[0] = S;
523 tmpvecs[1] = R;
524 PetscCall(VecMDot(Z, 2, tmpvecs, tmp));
525 delta = tmp[0];
526 beta = tmp[1]; /* delta <- z'*A*z = r'*B'*A*B*r */
527 KSPCheckDot(ksp, beta); /* beta <- z'*r */
528 }
529
530 i++;
531 } while (i < ksp->max_it);
532 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
533 PetscFunctionReturn(PETSC_SUCCESS);
534 }
535
536 /*
537 KSPDestroy_CG - Frees resources allocated in KSPSetup_CG and clears function
538 compositions from KSPCreate_CG. If adding your own KSP implementation,
539 you must be sure to free all allocated resources here to prevent
540 leaks.
541 */
KSPDestroy_CG(KSP ksp)542 PetscErrorCode KSPDestroy_CG(KSP ksp)
543 {
544 KSP_CG *cg = (KSP_CG *)ksp->data;
545
546 PetscFunctionBegin;
547 PetscCall(PetscFree4(cg->e, cg->d, cg->ee, cg->dd));
548 PetscCall(KSPDestroyDefault(ksp));
549 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", NULL));
550 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
551 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", NULL));
552 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", NULL));
553 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
554 PetscFunctionReturn(PETSC_SUCCESS);
555 }
556
557 /*
558 KSPView_CG - Prints information about the current Krylov method being used.
559 If your Krylov method has special options or flags that information
560 should be printed here.
561 */
KSPView_CG(KSP ksp,PetscViewer viewer)562 PetscErrorCode KSPView_CG(KSP ksp, PetscViewer viewer)
563 {
564 KSP_CG *cg = (KSP_CG *)ksp->data;
565 PetscBool isascii;
566
567 PetscFunctionBegin;
568 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
569 if (isascii) {
570 #if defined(PETSC_USE_COMPLEX)
571 PetscCall(PetscViewerASCIIPrintf(viewer, " variant %s\n", KSPCGTypes[cg->type]));
572 #endif
573 if (cg->singlereduction) PetscCall(PetscViewerASCIIPrintf(viewer, " using single-reduction variant\n"));
574 }
575 PetscFunctionReturn(PETSC_SUCCESS);
576 }
577
578 /*
579 KSPSetFromOptions_CG - Checks the options database for options related to the
580 conjugate gradient method.
581 */
KSPSetFromOptions_CG(KSP ksp,PetscOptionItems PetscOptionsObject)582 PetscErrorCode KSPSetFromOptions_CG(KSP ksp, PetscOptionItems PetscOptionsObject)
583 {
584 KSP_CG *cg = (KSP_CG *)ksp->data;
585 PetscBool flg;
586
587 PetscFunctionBegin;
588 PetscOptionsHeadBegin(PetscOptionsObject, "KSP CG and CGNE options");
589 #if defined(PETSC_USE_COMPLEX)
590 PetscCall(PetscOptionsEnum("-ksp_cg_type", "Matrix is Hermitian or complex symmetric", "KSPCGSetType", KSPCGTypes, (PetscEnum)cg->type, (PetscEnum *)&cg->type, NULL));
591 #endif
592 PetscCall(PetscOptionsBool("-ksp_cg_single_reduction", "Merge inner products into single MPI_Allreduce()", "KSPCGUseSingleReduction", cg->singlereduction, &cg->singlereduction, &flg));
593 if (flg) PetscCall(KSPCGUseSingleReduction(ksp, cg->singlereduction));
594 PetscOptionsHeadEnd();
595 PetscFunctionReturn(PETSC_SUCCESS);
596 }
597
598 /*
599 KSPCGSetType_CG - This is an option that is SPECIFIC to this particular Krylov method.
600 This routine is registered below in KSPCreate_CG() and called from the
601 routine KSPCGSetType() (see the file cgtype.c).
602 */
KSPCGSetType_CG(KSP ksp,KSPCGType type)603 PetscErrorCode KSPCGSetType_CG(KSP ksp, KSPCGType type)
604 {
605 KSP_CG *cg = (KSP_CG *)ksp->data;
606
607 PetscFunctionBegin;
608 cg->type = type;
609 PetscFunctionReturn(PETSC_SUCCESS);
610 }
611
612 /*
613 KSPCGUseSingleReduction_CG
614
615 This routine sets a flag to use a variant of CG. Note that (in somewhat
616 atypical fashion) it also swaps out the routine called when KSPSolve()
617 is invoked.
618 */
KSPCGUseSingleReduction_CG(KSP ksp,PetscBool flg)619 static PetscErrorCode KSPCGUseSingleReduction_CG(KSP ksp, PetscBool flg)
620 {
621 KSP_CG *cg = (KSP_CG *)ksp->data;
622
623 PetscFunctionBegin;
624 cg->singlereduction = flg;
625 if (cg->singlereduction) {
626 ksp->ops->solve = KSPSolve_CG_SingleReduction;
627 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
628 } else {
629 ksp->ops->solve = KSPSolve_CG;
630 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
631 }
632 PetscFunctionReturn(PETSC_SUCCESS);
633 }
634
KSPBuildResidual_CG(KSP ksp,Vec t,Vec v,Vec * V)635 PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP ksp, Vec t, Vec v, Vec *V)
636 {
637 PetscFunctionBegin;
638 PetscCall(VecCopy(ksp->work[0], v));
639 *V = v;
640 PetscFunctionReturn(PETSC_SUCCESS);
641 }
642
643 /*MC
644 KSPCG - The Preconditioned Conjugate Gradient (PCG) iterative method {cite}`hs:52` and {cite}`malek2014preconditioning` for solving linear systems using `KSP`.
645
646 Options Database Keys:
647 + -ksp_cg_type Hermitian - (for complex matrices only) indicates the matrix is Hermitian, see `KSPCGSetType()`
648 . -ksp_cg_type symmetric - (for complex matrices only) indicates the matrix is symmetric
649 - -ksp_cg_single_reduction - performs both inner products needed in the algorithm with a single `MPI_Allreduce()` call, see `KSPCGUseSingleReduction()`
650
651 Level: beginner
652
653 Notes:
654 The `KSPCG` method requires both the matrix and preconditioner to be symmetric positive (or negative) (semi) definite.
655
656 `KSPCG` is the best Krylov method, `KSPType`, when the matrix and preconditioner are symmetric positive definite (SPD).
657
658 Only left preconditioning is supported with `KSPCG`; there are several ways to motivate preconditioned CG, but they all produce the same algorithm.
659 One can interpret preconditioning $A$ with $B$ to mean any of the following\:
660 .vb
661 (1) Solve a left-preconditioned system $BAx = Bb $, using $ B^{-1}$ to define an inner product in the algorithm.
662 (2) Solve a right-preconditioned system $ABy = b, x = By,$ using $B$ to define an inner product in the algorithm.
663 (3) Solve a symmetrically-preconditioned system, $ E^TAEy = E^Tb, x = Ey, $ where $B = EE^T.$
664 (4) Solve $Ax=b$ with CG, but use the inner product defined by $B$ to define the method.
665 In all cases, the resulting algorithm only requires application of $B$ to vectors, the other inner-product does not appear explicitly in the code
666 .ve
667
668 For complex numbers there are two different CG methods, one for Hermitian symmetric matrices and one for non-Hermitian symmetric matrices. Use
669 `KSPCGSetType()` to indicate which type you are using.
670
671 One can use `KSPSetComputeEigenvalues()` and `KSPComputeEigenvalues()` to compute the eigenvalues of the (preconditioned) operator
672
673 There are two pipelined implementations of CG in PETSc `KSPPIPECG` and `KSPGROPPCG`. These may perform better for very large
674 numbers of MPI processes since they overlap communication and computation so the reduction operations in CG, that is inner products and norms,
675 do not dominate the compute time.
676
677 Developer Note:
678 KSPSolve_CG() should actually query the matrix to determine if it is Hermitian or symmetric and NOT require the user to
679 indicate it to the `KSP` object.
680
681 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSetComputeEigenvalues()`, `KSPComputeEigenvalues()`
682 `KSPCGSetType()`, `KSPCGUseSingleReduction()`, `KSPPIPECG`, `KSPGROPPCG`
683 M*/
684
685 /*
686 KSPCreate_CG - Creates the data structure for the Krylov method CG and sets the
687 function pointers for all the routines it needs to call (KSPSolve_CG() etc)
688
689 It must be labeled as PETSC_EXTERN to be dynamically linkable in C++
690 */
KSPCreate_CG(KSP ksp)691 PETSC_EXTERN PetscErrorCode KSPCreate_CG(KSP ksp)
692 {
693 KSP_CG *cg;
694
695 PetscFunctionBegin;
696 PetscCall(PetscNew(&cg));
697 cg->type = !PetscDefined(USE_COMPLEX) ? KSP_CG_SYMMETRIC : KSP_CG_HERMITIAN;
698 cg->obj_min = 0.0;
699 ksp->data = (void *)cg;
700
701 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
702 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
703 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
704 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
705
706 /*
707 Sets the functions that are associated with this data structure
708 (in C++ this is the same as defining virtual functions)
709 */
710 ksp->ops->setup = KSPSetUp_CG;
711 ksp->ops->solve = KSPSolve_CG;
712 ksp->ops->destroy = KSPDestroy_CG;
713 ksp->ops->view = KSPView_CG;
714 ksp->ops->setfromoptions = KSPSetFromOptions_CG;
715 ksp->ops->buildsolution = KSPBuildSolutionDefault;
716 ksp->ops->buildresidual = KSPBuildResidual_CG;
717
718 /*
719 Attach the function KSPCGSetType_CG() to this object. The routine
720 KSPCGSetType() checks for this attached function and calls it if it finds
721 it. (Sort of like a dynamic member function that can be added at run time
722 */
723 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CG));
724 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", KSPCGUseSingleReduction_CG));
725 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_CG));
726 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", KSPCGSetObjectiveTarget_CG));
727 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
728 PetscFunctionReturn(PETSC_SUCCESS);
729 }
730