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: static PetscErrorCode KSPCGGetObjFcn_CG(KSP ksp, PetscReal *obj)
66: {
67: KSP_CG *cg = (KSP_CG *)ksp->data;
69: PetscFunctionBegin;
70: *obj = cg->obj;
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*
75: KSPSetUp_CG - Sets up the workspace needed by the CG method.
77: This is called once, usually automatically by KSPSolve() or KSPSetUp()
78: but can be called directly by KSPSetUp()
79: */
80: static PetscErrorCode KSPSetUp_CG(KSP ksp)
81: {
82: KSP_CG *cgP = (KSP_CG *)ksp->data;
83: PetscInt maxit = ksp->max_it, nwork = 3;
85: PetscFunctionBegin;
86: /* get work vectors needed by CG */
87: if (cgP->singlereduction) nwork += 2;
88: PetscCall(KSPSetWorkVecs(ksp, nwork));
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));
98: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
99: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
100: }
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
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))
109: /*
110: KSPSolve_CG - This routine actually applies the conjugate gradient method
112: Note : this routine can be replaced with another one (see below) which implements
113: another variant of CG.
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: */
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 testobj;
130: PetscFunctionBegin;
131: cg = (KSP_CG *)ksp->data;
132: eigs = ksp->calc_sings;
133: stored_max_it = ksp->max_it;
134: X = ksp->vec_sol;
135: B = ksp->vec_rhs;
136: R = ksp->work[0];
137: Z = ksp->work[1];
138: P = ksp->work[2];
139: W = Z;
140: r2 = PetscSqr(cg->radius);
142: if (eigs) {
143: e = cg->e;
144: d = cg->d;
145: e[0] = 0.0;
146: }
147: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
149: ksp->its = 0;
150: if (!ksp->guess_zero) {
151: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
153: PetscCall(VecAYPX(R, -1.0, B));
154: if (cg->radius) { /* XXX direction? */
155: PetscCall(VecNorm(X, NORM_2, &norm_d));
156: norm_d *= norm_d;
157: }
158: } else {
159: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
160: norm_d = 0.0;
161: }
162: /* 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 */
163: PetscCall(VecFlag(R, ksp->reason == KSP_DIVERGED_PC_FAILED));
165: switch (ksp->normtype) {
166: case KSP_NORM_PRECONDITIONED:
167: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
168: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z = e'*A'*B'*B*A*e */
169: KSPCheckNorm(ksp, dp);
170: break;
171: case KSP_NORM_UNPRECONDITIONED:
172: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r = e'*A'*A*e */
173: KSPCheckNorm(ksp, dp);
174: break;
175: case KSP_NORM_NATURAL:
176: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
177: PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
178: KSPCheckDot(ksp, beta);
179: dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
180: break;
181: case KSP_NORM_NONE:
182: dp = 0.0;
183: break;
184: default:
185: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
186: }
188: /* Initialize objective function
189: obj = 1/2 x^T A x - x^T b */
190: testobj = (PetscBool)(cg->obj_min < 0.0);
191: PetscCall(VecXDot(R, X, &a));
192: cg->obj = 0.5 * PetscRealPart(a);
193: PetscCall(VecXDot(B, X, &a));
194: cg->obj -= 0.5 * PetscRealPart(a);
196: if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", ksp->its, (double)cg->obj));
197: PetscCall(KSPLogResidualHistory(ksp, dp));
198: PetscCall(KSPMonitor(ksp, ksp->its, dp));
199: ksp->rnorm = dp;
201: PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
203: if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
204: PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
205: ksp->reason = KSP_CONVERGED_ATOL;
206: }
208: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
210: if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
211: if (ksp->normtype != KSP_NORM_NATURAL) {
212: PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
213: KSPCheckDot(ksp, beta);
214: }
216: i = 0;
217: do {
218: ksp->its = i + 1;
219: if (beta == 0.0) {
220: ksp->reason = KSP_CONVERGED_ATOL;
221: PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
222: break;
223: #if !PetscDefined(USE_COMPLEX)
224: } else if (i > 0 && beta * betaold < 0.0) {
225: 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));
226: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
227: PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
228: break;
229: #endif
230: }
231: if (!i) {
232: PetscCall(VecCopy(Z, P)); /* p <- z */
233: if (cg->radius) {
234: PetscCall(VecNorm(P, NORM_2, &norm_p));
235: norm_p *= norm_p;
236: dMp = 0.0;
237: if (!ksp->guess_zero) PetscCall(VecDotRealPart(X, P, &dMp));
238: }
239: b = 0.0;
240: } else {
241: b = beta / betaold;
242: if (eigs) {
243: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
244: e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
245: }
246: PetscCall(VecAYPX(P, b, Z)); /* p <- z + b* p */
247: if (cg->radius) {
248: PetscCall(VecDotRealPart(X, P, &dMp));
249: PetscCall(VecNorm(P, NORM_2, &norm_p));
250: norm_p *= norm_p;
251: }
252: }
253: dpiold = dpi;
254: PetscCall(KSP_MatMult(ksp, Amat, P, W)); /* w <- Ap */
255: PetscCall(VecXDot(P, W, &dpi)); /* dpi <- p'w */
256: KSPCheckDot(ksp, dpi);
257: betaold = beta;
259: if ((dpi == 0.0) || ((i > 0) && ((PetscSign(PetscRealPart(dpi)) * PetscSign(PetscRealPart(dpiold))) < 0.0))) {
260: if (cg->radius) {
261: a = 0.0;
262: if (i == 0) {
263: if (norm_p > 0.0) {
264: a = PetscSqrtReal(r2 / norm_p);
265: } else {
266: PetscCall(VecNorm(R, NORM_2, &dp));
267: a = cg->radius > dp ? 1.0 : cg->radius / dp;
268: }
269: } else if (norm_p > 0.0) {
270: a = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
271: }
272: PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
273: cg->obj += PetscRealPart(a * (0.5 * a * dpi - betaold));
274: }
275: if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " N obj %g\n", i + 1, (double)cg->obj));
276: if (ksp->converged_neg_curve) {
277: PetscCall(PetscInfo(ksp, "converged due to negative curvature: %g\n", (double)(PetscRealPart(dpi))));
278: ksp->reason = KSP_CONVERGED_NEG_CURVE;
279: } else {
280: 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));
281: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
282: PetscCall(PetscInfo(ksp, "diverging due to indefinite matrix\n"));
283: }
284: break;
285: }
286: a = beta / dpi; /* a = beta/p'w */
287: if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
288: if (cg->radius) { /* Steihaugh-Toint */
289: PetscReal norm_dp1 = norm_d + PetscRealPart(a) * (2.0 * dMp + PetscRealPart(a) * norm_p);
290: if (norm_dp1 > r2) {
291: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
292: PetscCall(PetscInfo(ksp, "converged to the trust region radius %g\n", (double)cg->radius));
293: if (norm_p > 0.0) {
294: dp = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
295: PetscCall(VecAXPY(X, dp, P)); /* x <- x + ap */
296: cg->obj += PetscRealPart(dp * (0.5 * dp * dpi - beta));
297: }
298: if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " R obj %g\n", i + 1, (double)cg->obj));
299: break;
300: }
301: }
302: PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
303: PetscCall(VecAXPY(R, -a, W)); /* r <- r - aw */
304: if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
305: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
306: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z */
307: KSPCheckNorm(ksp, dp);
308: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
309: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r */
310: KSPCheckNorm(ksp, dp);
311: } else if (ksp->normtype == KSP_NORM_NATURAL) {
312: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
313: PetscCall(VecXDot(Z, R, &beta)); /* beta <- r'*z */
314: KSPCheckDot(ksp, beta);
315: dp = PetscSqrtReal(PetscAbsScalar(beta));
316: } else {
317: dp = 0.0;
318: }
319: cg->obj -= PetscRealPart(0.5 * a * betaold);
320: if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", i + 1, (double)cg->obj));
322: ksp->rnorm = dp;
323: PetscCall(KSPLogResidualHistory(ksp, dp));
324: PetscCall(KSPMonitor(ksp, i + 1, dp));
325: PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
327: if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
328: PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
329: ksp->reason = KSP_CONVERGED_ATOL;
330: }
332: if (ksp->reason) break;
334: if (cg->radius) {
335: PetscCall(VecNorm(X, NORM_2, &norm_d));
336: norm_d *= norm_d;
337: }
339: if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
340: if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
341: PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
342: KSPCheckDot(ksp, beta);
343: }
345: i++;
346: } while (i < ksp->max_it);
347: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*
352: KSPSolve_CG_SingleReduction
354: This variant of CG is identical in exact arithmetic to the standard algorithm,
355: but is rearranged to use only a single reduction stage per iteration, using additional
356: intermediate vectors.
358: See KSPCGUseSingleReduction_CG()
360: */
361: static PetscErrorCode KSPSolve_CG_SingleReduction(KSP ksp)
362: {
363: PetscInt i, stored_max_it, eigs;
364: PetscScalar dpi = 0.0, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL, delta, dpiold, tmp[2];
365: PetscReal dp = 0.0;
366: Vec X, B, Z, R, P, S, W, tmpvecs[2];
367: KSP_CG *cg;
368: Mat Amat, Pmat;
370: PetscFunctionBegin;
371: PetscCheck(ksp->nwork == 5, PetscObjectComm((PetscObject)ksp), PETSC_ERR_COR, "Unexpected number of work vectors %" PetscInt_FMT " != 5", ksp->nwork);
372: cg = (KSP_CG *)ksp->data;
373: eigs = ksp->calc_sings;
374: stored_max_it = ksp->max_it;
375: X = ksp->vec_sol;
376: B = ksp->vec_rhs;
377: R = ksp->work[0];
378: Z = ksp->work[1];
379: P = ksp->work[2];
380: S = ksp->work[3];
381: W = ksp->work[4];
383: if (eigs) {
384: e = cg->e;
385: d = cg->d;
386: e[0] = 0.0;
387: }
388: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
390: ksp->its = 0;
391: if (!ksp->guess_zero) {
392: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
393: PetscCall(VecAYPX(R, -1.0, B));
394: } else {
395: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
396: }
398: switch (ksp->normtype) {
399: case KSP_NORM_PRECONDITIONED:
400: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
401: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z = e'*A'*B'*B*A'*e' */
402: KSPCheckNorm(ksp, dp);
403: break;
404: case KSP_NORM_UNPRECONDITIONED:
405: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r = e'*A'*A*e */
406: KSPCheckNorm(ksp, dp);
407: break;
408: case KSP_NORM_NATURAL:
409: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
410: PetscCall(KSP_MatMult(ksp, Amat, Z, S));
411: PetscCall(VecXDot(Z, S, &delta)); /* delta <- z'*A*z = r'*B*A*B*r */
412: PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
413: KSPCheckDot(ksp, beta);
414: dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
415: break;
416: case KSP_NORM_NONE:
417: dp = 0.0;
418: break;
419: default:
420: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
421: }
422: PetscCall(KSPLogResidualHistory(ksp, dp));
423: PetscCall(KSPMonitor(ksp, 0, dp));
424: ksp->rnorm = dp;
426: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
427: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
429: if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
430: if (ksp->normtype != KSP_NORM_NATURAL) {
431: PetscCall(KSP_MatMult(ksp, Amat, Z, S));
432: PetscCall(VecXDot(Z, S, &delta)); /* delta <- z'*A*z = r'*B*A*B*r */
433: PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
434: KSPCheckDot(ksp, beta);
435: }
437: i = 0;
438: do {
439: ksp->its = i + 1;
440: if (beta == 0.0) {
441: ksp->reason = KSP_CONVERGED_ATOL;
442: PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
443: break;
444: #if !PetscDefined(USE_COMPLEX)
445: } else if (i > 0 && beta * betaold < 0.0) {
446: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner");
447: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
448: PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
449: break;
450: #endif
451: }
452: if (!i) {
453: PetscCall(VecCopy(Z, P)); /* p <- z */
454: b = 0.0;
455: } else {
456: b = beta / betaold;
457: if (eigs) {
458: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
459: e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
460: }
461: PetscCall(VecAYPX(P, b, Z)); /* p <- z + b* p */
462: }
463: dpiold = dpi;
464: if (!i) {
465: PetscCall(KSP_MatMult(ksp, Amat, P, W)); /* w <- Ap */
466: PetscCall(VecXDot(P, W, &dpi)); /* dpi <- p'w */
467: } else {
468: PetscCall(VecAYPX(W, beta / betaold, S)); /* w <- Ap */
469: dpi = delta - beta * beta * dpiold / (betaold * betaold); /* dpi <- p'w */
470: }
471: betaold = beta;
472: KSPCheckDot(ksp, beta);
474: if ((dpi == 0.0) || ((i > 0) && (PetscRealPart(dpi * dpiold) <= 0.0))) {
475: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix");
476: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
477: PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n"));
478: break;
479: }
480: a = beta / dpi; /* a = beta/p'w */
481: if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
482: PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
483: PetscCall(VecAXPY(R, -a, W)); /* r <- r - aw */
484: if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
485: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
486: PetscCall(KSP_MatMult(ksp, Amat, Z, S));
487: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z */
488: KSPCheckNorm(ksp, dp);
489: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
490: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r */
491: KSPCheckNorm(ksp, dp);
492: } else if (ksp->normtype == KSP_NORM_NATURAL) {
493: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
494: tmpvecs[0] = S;
495: tmpvecs[1] = R;
496: PetscCall(KSP_MatMult(ksp, Amat, Z, S));
497: PetscCall(VecMDot(Z, 2, tmpvecs, tmp)); /* delta <- z'*A*z = r'*B*A*B*r */
498: delta = tmp[0];
499: beta = tmp[1]; /* beta <- z'*r */
500: KSPCheckDot(ksp, beta);
501: dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
502: } else {
503: dp = 0.0;
504: }
505: ksp->rnorm = dp;
506: PetscCall(KSPLogResidualHistory(ksp, dp));
507: PetscCall(KSPMonitor(ksp, i + 1, dp));
508: PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
509: if (ksp->reason) break;
511: if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) {
512: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
513: PetscCall(KSP_MatMult(ksp, Amat, Z, S));
514: }
515: if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
516: tmpvecs[0] = S;
517: tmpvecs[1] = R;
518: PetscCall(VecMDot(Z, 2, tmpvecs, tmp));
519: delta = tmp[0];
520: beta = tmp[1]; /* delta <- z'*A*z = r'*B'*A*B*r */
521: KSPCheckDot(ksp, beta); /* beta <- z'*r */
522: }
524: i++;
525: } while (i < ksp->max_it);
526: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*
531: KSPDestroy_CG - Frees resources allocated in KSPSetup_CG and clears function
532: compositions from KSPCreate_CG. If adding your own KSP implementation,
533: you must be sure to free all allocated resources here to prevent
534: leaks.
535: */
536: PetscErrorCode KSPDestroy_CG(KSP ksp)
537: {
538: KSP_CG *cg = (KSP_CG *)ksp->data;
540: PetscFunctionBegin;
541: PetscCall(PetscFree4(cg->e, cg->d, cg->ee, cg->dd));
542: PetscCall(KSPDestroyDefault(ksp));
543: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", NULL));
544: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
545: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", NULL));
546: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", NULL));
547: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: /*
552: KSPView_CG - Prints information about the current Krylov method being used.
553: If your Krylov method has special options or flags that information
554: should be printed here.
555: */
556: PetscErrorCode KSPView_CG(KSP ksp, PetscViewer viewer)
557: {
558: KSP_CG *cg = (KSP_CG *)ksp->data;
559: PetscBool isascii;
561: PetscFunctionBegin;
562: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
563: if (isascii) {
564: if (PetscDefined(USE_COMPLEX)) PetscCall(PetscViewerASCIIPrintf(viewer, " variant %s\n", KSPCGTypes[cg->type]));
565: if (cg->singlereduction) PetscCall(PetscViewerASCIIPrintf(viewer, " using single-reduction variant\n"));
566: }
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*
571: KSPSetFromOptions_CG - Checks the options database for options related to the
572: conjugate gradient method.
573: */
574: PetscErrorCode KSPSetFromOptions_CG(KSP ksp, PetscOptionItems PetscOptionsObject)
575: {
576: KSP_CG *cg = (KSP_CG *)ksp->data;
577: PetscBool flg, flg2;
579: PetscFunctionBegin;
580: PetscOptionsHeadBegin(PetscOptionsObject, "KSP CG and CGNE options");
581: if (PetscDefined(USE_COMPLEX)) PetscCall(PetscOptionsEnum("-ksp_cg_type", "Matrix is Hermitian or complex symmetric", "KSPCGSetType", KSPCGTypes, (PetscEnum)cg->type, (PetscEnum *)&cg->type, NULL));
582: PetscCall(PetscOptionsBool("-ksp_cg_single_reduction", "Merge inner products into single MPI_Allreduce()", "KSPCGUseSingleReduction", cg->singlereduction, &flg2, &flg));
583: if (flg) PetscCall(KSPCGUseSingleReduction(ksp, flg2));
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: if (cg->singlereduction != flg) ksp->setupstage = KSP_SETUP_NEW;
615: cg->singlereduction = flg;
616: if (cg->singlereduction) {
617: ksp->ops->solve = KSPSolve_CG_SingleReduction;
618: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
619: } else {
620: ksp->ops->solve = KSPSolve_CG;
621: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
622: }
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP ksp, Vec t, Vec v, Vec *V)
627: {
628: PetscFunctionBegin;
629: PetscCall(VecCopy(ksp->work[0], v));
630: *V = v;
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: /*MC
635: KSPCG - The Preconditioned Conjugate Gradient (PCG) iterative method {cite}`hs:52` and {cite}`malek2014preconditioning` for solving linear systems using `KSP`.
637: Options Database Keys:
638: + -ksp_cg_type (hermitian|symmetric) - (for complex matrices only) indicates the matrix is Hermitian or symmetric, see `KSPCGSetType()`
639: - -ksp_cg_single_reduction - performs both inner products needed in the algorithm with a single `MPI_Allreduce()` call, see `KSPCGUseSingleReduction()`
641: Level: beginner
643: Notes:
644: The `KSPCG` method requires both the matrix and preconditioner to be symmetric positive (or negative) (semi) definite.
646: `KSPCG` is the best Krylov method, `KSPType`, when the matrix and preconditioner are symmetric positive definite (SPD).
648: Only left preconditioning is supported with `KSPCG`; there are several ways to motivate preconditioned CG, but they all produce the same algorithm.
649: One can interpret preconditioning $A$ with $B$ to mean any of the following\:
650: .vb
651: (1) Solve a left-preconditioned system $BAx = Bb $, using $ B^{-1}$ to define an inner product in the algorithm.
652: (2) Solve a right-preconditioned system $ABy = b, x = By,$ using $B$ to define an inner product in the algorithm.
653: (3) Solve a symmetrically-preconditioned system, $ E^TAEy = E^Tb, x = Ey, $ where $B = EE^T.$
654: (4) Solve $Ax=b$ with CG, but use the inner product defined by $B$ to define the method.
655: In all cases, the resulting algorithm only requires application of $B$ to vectors, the other inner-product does not appear explicitly in the code
656: .ve
658: For complex numbers there are two different CG methods, one for Hermitian symmetric matrices and one for non-Hermitian symmetric matrices. Use
659: `KSPCGSetType()` to indicate which type you are using.
661: One can use `KSPSetComputeEigenvalues()` and `KSPComputeEigenvalues()` to compute the eigenvalues of the (preconditioned) operator
663: There are two pipelined implementations of CG in PETSc `KSPPIPECG` and `KSPGROPPCG`. These may perform better for very large
664: numbers of MPI processes since they overlap communication and computation so the reduction operations in CG, that is inner products and norms,
665: do not dominate the compute time.
667: Developer Note:
668: KSPSolve_CG() should actually query the matrix to determine if it is Hermitian or symmetric and NOT require the user to
669: indicate it to the `KSP` object.
671: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSetComputeEigenvalues()`, `KSPComputeEigenvalues()`,
672: `KSPCGSetType()`, `KSPCGUseSingleReduction()`, `KSPPIPECG`, `KSPGROPPCG`
673: M*/
675: /*
676: KSPCreate_CG - Creates the data structure for the Krylov method CG and sets the
677: function pointers for all the routines it needs to call (KSPSolve_CG() etc)
679: It must be labeled as PETSC_EXTERN to be dynamically linkable in C++
680: */
681: PETSC_EXTERN PetscErrorCode KSPCreate_CG(KSP ksp)
682: {
683: KSP_CG *cg;
685: PetscFunctionBegin;
686: PetscCall(PetscNew(&cg));
687: cg->type = !PetscDefined(USE_COMPLEX) ? KSP_CG_SYMMETRIC : KSP_CG_HERMITIAN;
688: cg->obj_min = 0.0;
689: ksp->data = (void *)cg;
691: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
692: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
693: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
694: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
696: /*
697: Sets the functions that are associated with this data structure
698: (in C++ this is the same as defining virtual functions)
699: */
700: ksp->ops->setup = KSPSetUp_CG;
701: ksp->ops->solve = KSPSolve_CG;
702: ksp->ops->destroy = KSPDestroy_CG;
703: ksp->ops->view = KSPView_CG;
704: ksp->ops->setfromoptions = KSPSetFromOptions_CG;
705: ksp->ops->buildsolution = KSPBuildSolutionDefault;
706: ksp->ops->buildresidual = KSPBuildResidual_CG;
708: /*
709: Attach the function KSPCGSetType_CG() to this object. The routine
710: KSPCGSetType() checks for this attached function and calls it if it finds
711: it. (Sort of like a dynamic member function that can be added at run time
712: */
713: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CG));
714: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", KSPCGUseSingleReduction_CG));
715: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_CG));
716: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", KSPCGSetObjectiveTarget_CG));
717: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }