Actual source code: groppcg.c
1: #include <petsc/private/kspimpl.h>
3: /*
4: KSPSetUp_GROPPCG - Sets up the workspace needed by the GROPPCG method.
6: This is called once, usually automatically by KSPSolve() or KSPSetUp()
7: but can be called directly by KSPSetUp()
8: */
9: static PetscErrorCode KSPSetUp_GROPPCG(KSP ksp)
10: {
11: KSPSetWorkVecs(ksp, 6);
12: return 0;
13: }
15: /*
16: KSPSolve_GROPPCG
18: Input Parameter:
19: . ksp - the Krylov space object that was set to use conjugate gradient, by, for
20: example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
21: */
22: static PetscErrorCode KSPSolve_GROPPCG(KSP ksp)
23: {
24: PetscInt i;
25: PetscScalar alpha, beta = 0.0, gamma, gammaNew, t;
26: PetscReal dp = 0.0;
27: Vec x, b, r, p, s, S, z, Z;
28: Mat Amat, Pmat;
29: PetscBool diagonalscale;
31: PCGetDiagonalScale(ksp->pc, &diagonalscale);
34: x = ksp->vec_sol;
35: b = ksp->vec_rhs;
36: r = ksp->work[0];
37: p = ksp->work[1];
38: s = ksp->work[2];
39: S = ksp->work[3];
40: z = ksp->work[4];
41: Z = ksp->work[5];
43: PCGetOperators(ksp->pc, &Amat, &Pmat);
45: ksp->its = 0;
46: if (!ksp->guess_zero) {
47: KSP_MatMult(ksp, Amat, x, r); /* r <- b - Ax */
48: VecAYPX(r, -1.0, b);
49: } else {
50: VecCopy(b, r); /* r <- b (x is 0) */
51: }
53: KSP_PCApply(ksp, r, z); /* z <- Br */
54: VecCopy(z, p); /* p <- z */
55: VecDotBegin(r, z, &gamma); /* gamma <- z'*r */
56: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
57: KSP_MatMult(ksp, Amat, p, s); /* s <- Ap */
58: VecDotEnd(r, z, &gamma); /* gamma <- z'*r */
60: switch (ksp->normtype) {
61: case KSP_NORM_PRECONDITIONED:
62: /* This could be merged with the computation of gamma above */
63: VecNorm(z, NORM_2, &dp); /* dp <- z'*z = e'*A'*B'*B*A'*e' */
64: break;
65: case KSP_NORM_UNPRECONDITIONED:
66: /* This could be merged with the computation of gamma above */
67: VecNorm(r, NORM_2, &dp); /* dp <- r'*r = e'*A'*A*e */
68: break;
69: case KSP_NORM_NATURAL:
70: KSPCheckDot(ksp, gamma);
71: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
72: break;
73: case KSP_NORM_NONE:
74: dp = 0.0;
75: break;
76: default:
77: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
78: }
79: KSPLogResidualHistory(ksp, dp);
80: KSPMonitor(ksp, 0, dp);
81: ksp->rnorm = dp;
82: (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP); /* test for convergence */
83: if (ksp->reason) return 0;
85: i = 0;
86: do {
87: ksp->its = i + 1;
88: i++;
90: VecDotBegin(p, s, &t);
91: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));
93: KSP_PCApply(ksp, s, S); /* S <- Bs */
95: VecDotEnd(p, s, &t);
97: alpha = gamma / t;
98: VecAXPY(x, alpha, p); /* x <- x + alpha * p */
99: VecAXPY(r, -alpha, s); /* r <- r - alpha * s */
100: VecAXPY(z, -alpha, S); /* z <- z - alpha * S */
102: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
103: VecNormBegin(r, NORM_2, &dp);
104: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
105: VecNormBegin(z, NORM_2, &dp);
106: }
107: VecDotBegin(r, z, &gammaNew);
108: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
110: KSP_MatMult(ksp, Amat, z, Z); /* Z <- Az */
112: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
113: VecNormEnd(r, NORM_2, &dp);
114: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
115: VecNormEnd(z, NORM_2, &dp);
116: }
117: VecDotEnd(r, z, &gammaNew);
119: if (ksp->normtype == KSP_NORM_NATURAL) {
120: KSPCheckDot(ksp, gammaNew);
121: dp = PetscSqrtReal(PetscAbsScalar(gammaNew)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
122: } else if (ksp->normtype == KSP_NORM_NONE) {
123: dp = 0.0;
124: }
125: ksp->rnorm = dp;
126: KSPLogResidualHistory(ksp, dp);
127: KSPMonitor(ksp, i, dp);
128: (*ksp->converged)(ksp, i, dp, &ksp->reason, ksp->cnvP);
129: if (ksp->reason) return 0;
131: beta = gammaNew / gamma;
132: gamma = gammaNew;
133: VecAYPX(p, beta, z); /* p <- z + beta * p */
134: VecAYPX(s, beta, Z); /* s <- Z + beta * s */
136: } while (i < ksp->max_it);
138: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
139: return 0;
140: }
142: PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP, Vec, Vec, Vec *);
144: /*MC
145: KSPGROPPCG - A pipelined conjugate gradient method developed by Bill Gropp. [](sec_pipelineksp)
147: Level: intermediate
149: Notes:
150: This method has two reductions, one of which is overlapped with the matrix-vector product and one of which is
151: overlapped with the preconditioner.
153: See also `KSPPIPECG`, which has only a single reduction that overlaps both the matrix-vector product and the preconditioner.
155: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
156: See [](doc_faq_pipelined)
158: Contributed by:
159: Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders
161: Reference:
162: http://www.cs.uiuc.edu/~wgropp/bib/talks/tdata/2012/icerm.pdf
164: .seealso: [](chapter_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPPIPECG2()`, `KSPSetType()`, `KSPPIPECG`, `KSPPIPECR`, `KSPPGMRES`, `KSPCG`, `KSPCGUseSingleReduction()`
165: M*/
167: PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
168: {
169: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
170: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
171: KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
172: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
174: ksp->ops->setup = KSPSetUp_GROPPCG;
175: ksp->ops->solve = KSPSolve_GROPPCG;
176: ksp->ops->destroy = KSPDestroyDefault;
177: ksp->ops->view = NULL;
178: ksp->ops->setfromoptions = NULL;
179: ksp->ops->buildsolution = KSPBuildSolutionDefault;
180: ksp->ops->buildresidual = KSPBuildResidual_CG;
181: return 0;
182: }