Actual source code: pipecgrr.c
1: #include <petsc/private/kspimpl.h>
3: /*
4: KSPSetUp_PIPECGRR - Sets up the workspace needed by the PIPECGRR 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_PIPECGRR(KSP ksp)
10: {
11: PetscFunctionBegin;
12: /* get work vectors needed by PIPECGRR */
13: PetscCall(KSPSetWorkVecs(ksp, 9));
14: PetscFunctionReturn(PETSC_SUCCESS);
15: }
17: /*
18: KSPSolve_PIPECGRR - This routine actually applies the pipelined conjugate gradient method with automated residual replacement
19: */
20: static PetscErrorCode KSPSolve_PIPECGRR(KSP ksp)
21: {
22: PetscInt i = 0, replace = 0, nsize;
23: PetscScalar alpha = 0.0, beta = 0.0, gamma = 0.0, gammaold = 0.0, delta = 0.0, alphap = 0.0, betap = 0.0;
24: PetscReal dp = 0.0, nsi = 0.0, sqn = 0.0, Anorm = 0.0, rnp = 0.0, pnp = 0.0, snp = 0.0, unp = 0.0, wnp = 0.0, xnp = 0.0, qnp = 0.0, znp = 0.0, mnz = 5.0, tol = PETSC_SQRT_MACHINE_EPSILON, eps = PETSC_MACHINE_EPSILON;
25: PetscReal ds = 0.0, dz = 0.0, dx = 0.0, dpp = 0.0, dq = 0.0, dm = 0.0, du = 0.0, dw = 0.0, db = 0.0, errr = 0.0, errrprev = 0.0, errs = 0.0, errw = 0.0, errz = 0.0, errncr = 0.0, errncs = 0.0, errncw = 0.0, errncz = 0.0;
26: Vec X, B, Z, P, W, Q, U, M, N, R, S;
27: Mat Amat, Pmat;
29: PetscFunctionBegin;
30: X = ksp->vec_sol;
31: B = ksp->vec_rhs;
32: M = ksp->work[0];
33: Z = ksp->work[1];
34: P = ksp->work[2];
35: N = ksp->work[3];
36: W = ksp->work[4];
37: Q = ksp->work[5];
38: U = ksp->work[6];
39: R = ksp->work[7];
40: S = ksp->work[8];
42: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
44: ksp->its = 0;
45: if (!ksp->guess_zero) {
46: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
47: PetscCall(VecAYPX(R, -1.0, B));
48: } else {
49: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
50: }
52: PetscCall(KSP_PCApply(ksp, R, U)); /* u <- Br */
54: switch (ksp->normtype) {
55: case KSP_NORM_PRECONDITIONED:
56: PetscCall(VecNormBegin(U, NORM_2, &dp)); /* dp <- u'*u = e'*A'*B'*B*A'*e' */
57: PetscCall(VecNormBegin(B, NORM_2, &db));
58: PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)U)));
59: PetscCall(KSP_MatMult(ksp, Amat, U, W)); /* w <- Au */
60: PetscCall(VecNormEnd(U, NORM_2, &dp));
61: PetscCall(VecNormEnd(B, NORM_2, &db));
62: break;
63: case KSP_NORM_UNPRECONDITIONED:
64: PetscCall(VecNormBegin(R, NORM_2, &dp)); /* dp <- r'*r = e'*A'*A*e */
65: PetscCall(VecNormBegin(B, NORM_2, &db));
66: PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R)));
67: PetscCall(KSP_MatMult(ksp, Amat, U, W)); /* w <- Au */
68: PetscCall(VecNormEnd(R, NORM_2, &dp));
69: PetscCall(VecNormEnd(B, NORM_2, &db));
70: break;
71: case KSP_NORM_NATURAL:
72: PetscCall(VecDotBegin(R, U, &gamma)); /* gamma <- u'*r */
73: PetscCall(VecNormBegin(B, NORM_2, &db));
74: PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R)));
75: PetscCall(KSP_MatMult(ksp, Amat, U, W)); /* w <- Au */
76: PetscCall(VecDotEnd(R, U, &gamma));
77: PetscCall(VecNormEnd(B, NORM_2, &db));
78: KSPCheckDot(ksp, gamma);
79: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- r'*u = r'*B*r = e'*A'*B*A*e */
80: break;
81: case KSP_NORM_NONE:
82: PetscCall(KSP_MatMult(ksp, Amat, U, W));
83: dp = 0.0;
84: break;
85: default:
86: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
87: }
88: PetscCall(KSPLogResidualHistory(ksp, dp));
89: PetscCall(KSPMonitor(ksp, 0, dp));
90: ksp->rnorm = dp;
91: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
92: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
94: PetscCall(MatNorm(Amat, NORM_INFINITY, &Anorm));
95: PetscCall(VecGetSize(B, &nsize));
96: nsi = (PetscReal)nsize;
97: sqn = PetscSqrtReal(nsi);
99: do {
100: if (i > 1) {
101: pnp = dpp;
102: snp = ds;
103: qnp = dq;
104: znp = dz;
105: }
106: if (i > 0) {
107: rnp = dp;
108: unp = du;
109: wnp = dw;
110: xnp = dx;
111: alphap = alpha;
112: betap = beta;
113: }
115: if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
116: PetscCall(VecNormBegin(R, NORM_2, &dp));
117: } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
118: PetscCall(VecNormBegin(U, NORM_2, &dp));
119: }
120: if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) PetscCall(VecDotBegin(R, U, &gamma));
121: PetscCall(VecDotBegin(W, U, &delta));
123: if (i > 0) {
124: PetscCall(VecNormBegin(S, NORM_2, &ds));
125: PetscCall(VecNormBegin(Z, NORM_2, &dz));
126: PetscCall(VecNormBegin(P, NORM_2, &dpp));
127: PetscCall(VecNormBegin(Q, NORM_2, &dq));
128: PetscCall(VecNormBegin(M, NORM_2, &dm));
129: }
130: PetscCall(VecNormBegin(X, NORM_2, &dx));
131: PetscCall(VecNormBegin(U, NORM_2, &du));
132: PetscCall(VecNormBegin(W, NORM_2, &dw));
134: PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R)));
135: PetscCall(KSP_PCApply(ksp, W, M)); /* m <- Bw */
136: PetscCall(KSP_MatMult(ksp, Amat, M, N)); /* n <- Am */
138: if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
139: PetscCall(VecNormEnd(R, NORM_2, &dp));
140: } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
141: PetscCall(VecNormEnd(U, NORM_2, &dp));
142: }
143: if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) PetscCall(VecDotEnd(R, U, &gamma));
144: PetscCall(VecDotEnd(W, U, &delta));
146: if (i > 0) {
147: PetscCall(VecNormEnd(S, NORM_2, &ds));
148: PetscCall(VecNormEnd(Z, NORM_2, &dz));
149: PetscCall(VecNormEnd(P, NORM_2, &dpp));
150: PetscCall(VecNormEnd(Q, NORM_2, &dq));
151: PetscCall(VecNormEnd(M, NORM_2, &dm));
152: }
153: PetscCall(VecNormEnd(X, NORM_2, &dx));
154: PetscCall(VecNormEnd(U, NORM_2, &du));
155: PetscCall(VecNormEnd(W, NORM_2, &dw));
157: if (i > 0) {
158: if (ksp->normtype == KSP_NORM_NATURAL) dp = PetscSqrtReal(PetscAbsScalar(gamma));
159: else if (ksp->normtype == KSP_NORM_NONE) dp = 0.0;
161: ksp->rnorm = dp;
162: PetscCall(KSPLogResidualHistory(ksp, dp));
163: PetscCall(KSPMonitor(ksp, i, dp));
164: PetscCall((*ksp->converged)(ksp, i, dp, &ksp->reason, ksp->cnvP));
165: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: if (i == 0) {
169: alpha = gamma / delta;
170: PetscCall(VecCopy(N, Z)); /* z <- n */
171: PetscCall(VecCopy(M, Q)); /* q <- m */
172: PetscCall(VecCopy(U, P)); /* p <- u */
173: PetscCall(VecCopy(W, S)); /* s <- w */
174: } else {
175: beta = gamma / gammaold;
176: alpha = gamma / (delta - beta / alpha * gamma);
177: PetscCall(VecAYPX(Z, beta, N)); /* z <- n + beta * z */
178: PetscCall(VecAYPX(Q, beta, M)); /* q <- m + beta * q */
179: PetscCall(VecAYPX(P, beta, U)); /* p <- u + beta * p */
180: PetscCall(VecAYPX(S, beta, W)); /* s <- w + beta * s */
181: }
182: PetscCall(VecAXPY(X, alpha, P)); /* x <- x + alpha * p */
183: PetscCall(VecAXPY(U, -alpha, Q)); /* u <- u - alpha * q */
184: PetscCall(VecAXPY(W, -alpha, Z)); /* w <- w - alpha * z */
185: PetscCall(VecAXPY(R, -alpha, S)); /* r <- r - alpha * s */
186: gammaold = gamma;
188: if (i > 0) {
189: errncr = PetscSqrtReal(Anorm * xnp + 2.0 * Anorm * PetscAbsScalar(alphap) * dpp + rnp + 2.0 * PetscAbsScalar(alphap) * ds) * eps;
190: errncw = PetscSqrtReal(Anorm * unp + 2.0 * Anorm * PetscAbsScalar(alphap) * dq + wnp + 2.0 * PetscAbsScalar(alphap) * dz) * eps;
191: }
192: if (i > 1) {
193: errncs = PetscSqrtReal(Anorm * unp + 2.0 * Anorm * PetscAbsScalar(betap) * pnp + wnp + 2.0 * PetscAbsScalar(betap) * snp) * eps;
194: errncz = PetscSqrtReal((mnz * sqn + 2) * Anorm * dm + 2.0 * Anorm * PetscAbsScalar(betap) * qnp + 2.0 * PetscAbsScalar(betap) * znp) * eps;
195: }
197: if (i > 0) {
198: if (i == 1) {
199: errr = PetscSqrtReal((mnz * sqn + 1) * Anorm * xnp + db) * eps + PetscSqrtReal(PetscAbsScalar(alphap) * mnz * sqn * Anorm * dpp) * eps + errncr;
200: errs = PetscSqrtReal(mnz * sqn * Anorm * dpp) * eps;
201: errw = PetscSqrtReal(mnz * sqn * Anorm * unp) * eps + PetscSqrtReal(PetscAbsScalar(alphap) * mnz * sqn * Anorm * dq) * eps + errncw;
202: errz = PetscSqrtReal(mnz * sqn * Anorm * dq) * eps;
203: } else if (replace == 1) {
204: errrprev = errr;
205: errr = PetscSqrtReal((mnz * sqn + 1) * Anorm * dx + db) * eps;
206: errs = PetscSqrtReal(mnz * sqn * Anorm * dpp) * eps;
207: errw = PetscSqrtReal(mnz * sqn * Anorm * du) * eps;
208: errz = PetscSqrtReal(mnz * sqn * Anorm * dq) * eps;
209: replace = 0;
210: } else {
211: errrprev = errr;
212: errr = errr + PetscAbsScalar(alphap) * PetscAbsScalar(betap) * errs + PetscAbsScalar(alphap) * errw + errncr + PetscAbsScalar(alphap) * errncs;
213: errs = errw + PetscAbsScalar(betap) * errs + errncs;
214: errw = errw + PetscAbsScalar(alphap) * PetscAbsScalar(betap) * errz + errncw + PetscAbsScalar(alphap) * errncz;
215: errz = PetscAbsScalar(betap) * errz + errncz;
216: }
217: if (i > 1 && errrprev <= (tol * rnp) && errr > (tol * dp)) {
218: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- Ax - b */
219: PetscCall(VecAYPX(R, -1.0, B));
220: PetscCall(KSP_PCApply(ksp, R, U)); /* u <- Br */
221: PetscCall(KSP_MatMult(ksp, Amat, U, W)); /* w <- Au */
222: PetscCall(KSP_MatMult(ksp, Amat, P, S)); /* s <- Ap */
223: PetscCall(KSP_PCApply(ksp, S, Q)); /* q <- Bs */
224: PetscCall(KSP_MatMult(ksp, Amat, Q, Z)); /* z <- Aq */
225: replace = 1;
226: }
227: }
229: i++;
230: ksp->its = i;
232: } while (i <= ksp->max_it);
233: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: /*MC
238: KSPPIPECGRR - Pipelined conjugate gradient method with automated residual replacements {cite}`cools2018analyzing`. [](sec_pipelineksp)
240: Level: intermediate
242: Notes:
243: This method has only a single non-blocking reduction per iteration, compared to 2 blocking for standard `KSPCG`. The
244: non-blocking reduction is overlapped by the matrix-vector product and preconditioner application.
246: `KSPPIPECGRR` improves the robustness of `KSPPIPECG` by adding an automated residual replacement strategy.
247: True residual and other auxiliary variables are computed explicitly in a number of dynamically determined
248: iterations to counteract the accumulation of rounding errors and thus attain a higher maximal final accuracy.
250: See also `KSPPIPECG`, which is identical to `KSPPIPECGRR` without residual replacements.
251: See also `KSPPIPECR`, where the reduction is only overlapped with the matrix-vector product.
253: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for
254: performance of pipelined methods. See [](doc_faq_pipelined)
256: Contributed by:
257: Siegfried Cools, Universiteit Antwerpen, Dept. Mathematics & Computer Science,
258: European FP7 Project on EXascale Algorithms and Advanced Computational Techniques (EXA2CT) / Research Foundation Flanders (FWO)
260: .seealso: [](ch_ksp), [](doc_faq_pipelined), [](sec_pipelineksp), `KSPCreate()`, `KSPSetType()`, `KSPPIPECR`, `KSPGROPPCG`, `KSPPIPECG`, `KSPPGMRES`, `KSPCG`, `KSPPIPEBCGS`, `KSPCGUseSingleReduction()`
261: M*/
262: PETSC_EXTERN PetscErrorCode KSPCreate_PIPECGRR(KSP ksp)
263: {
264: PetscFunctionBegin;
265: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
266: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
267: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
268: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
270: ksp->ops->setup = KSPSetUp_PIPECGRR;
271: ksp->ops->solve = KSPSolve_PIPECGRR;
272: ksp->ops->destroy = KSPDestroyDefault;
273: ksp->ops->view = NULL;
274: ksp->ops->setfromoptions = NULL;
275: ksp->ops->buildsolution = KSPBuildSolutionDefault;
276: ksp->ops->buildresidual = KSPBuildResidualDefault;
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }