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