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