Actual source code: pipecgrr.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/kspimpl.h>
4: /*
5: KSPSetUp_PIPECGRR - Sets up the workspace needed by the PIPECGRR method.
7: This is called once, usually automatically by KSPSolve() or KSPSetUp()
8: but can be called directly by KSPSetUp()
9: */
12: PetscErrorCode KSPSetUp_PIPECGRR(KSP ksp)
13: {
17: /* get work vectors needed by PIPECGRR */
18: KSPSetWorkVecs(ksp,9);
19: return(0);
20: }
22: /*
23: KSPSolve_PIPECGRR - This routine actually applies the pipelined conjugate gradient method with automated residual replacement
25: Input Parameter:
26: . ksp - the Krylov space object that was set to use conjugate gradient, by, for
27: example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
28: */
31: PetscErrorCode KSPSolve_PIPECGRR(KSP ksp)
32: {
34: PetscInt i,replace,totreplaces;
35: PetscScalar alpha = 0.0,beta = 0.0,gamma = 0.0,gammaold = 0.0,delta = 0.0;
36: PetscReal dp = 0.0,rnormprev = 0.0,rnormcurr = 0.0,tol = PETSC_SQRT_MACHINE_EPSILON,eps = PETSC_MACHINE_EPSILON,ds = 0.0,dz = 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;
37: Vec X,B,Z,P,W,Q,U,M,N,R,S;
38: Mat Amat,Pmat;
39: PetscBool diagonalscale;
42: PCGetDiagonalScale(ksp->pc,&diagonalscale);
43: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
45: X = ksp->vec_sol;
46: B = ksp->vec_rhs;
47: M = ksp->work[0];
48: Z = ksp->work[1];
49: P = ksp->work[2];
50: N = ksp->work[3];
51: W = ksp->work[4];
52: Q = ksp->work[5];
53: U = ksp->work[6];
54: R = ksp->work[7];
55: S = ksp->work[8];
57: PCGetOperators(ksp->pc,&Amat,&Pmat);
59: ksp->its = 0;
60: if (!ksp->guess_zero) {
61: KSP_MatMult(ksp,Amat,X,R); /* r <- b - Ax */
62: VecAYPX(R,-1.0,B);
63: } else {
64: VecCopy(B,R); /* r <- b (x is 0) */
65: }
67: KSP_PCApply(ksp,R,U); /* u <- Br */
69: switch (ksp->normtype) {
70: case KSP_NORM_PRECONDITIONED:
71: VecNormBegin(U,NORM_2,&dp); /* dp <- u'*u = e'*A'*B'*B*A'*e' */
72: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)U));
73: KSP_MatMult(ksp,Amat,U,W); /* w <- Au */
74: VecNormEnd(U,NORM_2,&dp);
75: break;
76: case KSP_NORM_UNPRECONDITIONED:
77: VecNormBegin(R,NORM_2,&dp); /* dp <- r'*r = e'*A'*A*e */
78: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));
79: KSP_MatMult(ksp,Amat,U,W); /* w <- Au */
80: VecNormEnd(R,NORM_2,&dp);
81: break;
82: case KSP_NORM_NATURAL:
83: VecDotBegin(R,U,&gamma); /* gamma <- u'*r */
84: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));
85: KSP_MatMult(ksp,Amat,U,W); /* w <- Au */
86: VecDotEnd(R,U,&gamma);
87: KSPCheckDot(ksp,gamma);
88: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- r'*u = r'*B*r = e'*A'*B*A*e */
89: break;
90: case KSP_NORM_NONE:
91: KSP_MatMult(ksp,Amat,U,W);
92: dp = 0.0;
93: break;
94: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
95: }
96: KSPLogResidualHistory(ksp,dp);
97: KSPMonitor(ksp,0,dp);
98: ksp->rnorm = dp;
99: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP); /* test for convergence */
100: if (ksp->reason) return(0);
102: i = 0;
103: replace = 0;
104: totreplaces = 0;
105: do {
106: rnormprev = dp;
108: if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
109: VecNormBegin(R,NORM_2,&dp);
110: } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
111: VecNormBegin(U,NORM_2,&dp);
112: }
113: if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) {
114: VecDotBegin(R,U,&gamma);
115: }
116: VecDotBegin(W,U,&delta);
117:
118: if (i > 0) {
119: VecNormBegin(S,NORM_2,&ds);
120: VecNormBegin(Z,NORM_2,&dz);
121: }
123: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));
124: KSP_PCApply(ksp,W,M); /* m <- Bw */
125: KSP_MatMult(ksp,Amat,M,N); /* n <- Am */
127: if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
128: VecNormEnd(R,NORM_2,&dp);
129: } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
130: VecNormEnd(U,NORM_2,&dp);
131: }
132: if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) {
133: VecDotEnd(R,U,&gamma);
134: }
135: VecDotEnd(W,U,&delta);
136:
137: if (i > 0) {
138: VecNormEnd(S,NORM_2,&ds);
139: VecNormEnd(Z,NORM_2,&dz);
140: }
142: if (i > 0) {
143: if (ksp->normtype == KSP_NORM_NATURAL) dp = PetscSqrtReal(PetscAbsScalar(gamma));
144: else if (ksp->normtype == KSP_NORM_NONE) dp = 0.0;
146: ksp->rnorm = dp;
147: KSPLogResidualHistory(ksp,dp);
148: KSPMonitor(ksp,i,dp);
149: (*ksp->converged)(ksp,i,dp,&ksp->reason,ksp->cnvP);
150: if (ksp->reason) break;
151: }
152: rnormcurr = dp;
154: if (i == 0) {
155: alpha = gamma / delta;
156: VecCopy(N,Z); /* z <- n */
157: VecCopy(M,Q); /* q <- m */
158: VecCopy(U,P); /* p <- u */
159: VecCopy(W,S); /* s <- w */
160: } else {
161: beta = gamma / gammaold;
162: alpha = gamma / (delta - beta / alpha * gamma);
163: VecAYPX(Z,beta,N); /* z <- n + beta * z */
164: VecAYPX(Q,beta,M); /* q <- m + beta * q */
165: VecAYPX(P,beta,U); /* p <- u + beta * p */
166: VecAYPX(S,beta,W); /* s <- w + beta * s */
167: }
168: VecAXPY(X, alpha,P); /* x <- x + alpha * p */
169: VecAXPY(U,-alpha,Q); /* u <- u - alpha * q */
170: VecAXPY(W,-alpha,Z); /* w <- w - alpha * z */
171: VecAXPY(R,-alpha,S); /* r <- r - alpha * s */
172: gammaold = gamma;
174: if (i > 0) {
175: errncr = 2.0 * PetscAbsScalar(alpha) * ds * eps;
176: errncs = 2.0 * PetscAbsScalar(beta) * ds * eps;
177: errncw = 2.0 * PetscAbsScalar(alpha) * dz * eps;
178: errncz = 2.0 * PetscAbsScalar(beta) * dz * eps;
179:
180: if (i == 1) {
181: errr = errncr;
182: errs = errncs;
183: errw = errncw;
184: errz = errncz;
185: } else if (replace == 1) {
186: errrprev = errr;
187: errr = errncr;
188: errs = errncs;
189: errw = errncw;
190: errz = errncz;
191: replace = 0;
192: } else {
193: errrprev = errr;
194: errr = errr + PetscAbsScalar(alpha) * errs + errncr;
195: errs = PetscAbsScalar(beta) * errs + errw + PetscAbsScalar(alpha) * errz + errncs + errncw;
196: errw = errw + PetscAbsScalar(alpha) * errz + errncw;
197: errz = PetscAbsScalar(beta) * errz + errncz;
198: }
200: if (i > 1 && errrprev <= (tol * rnormprev) && errr > (tol * rnormcurr)) {
201: KSP_MatMult(ksp,Amat,X,R); /* r <- Ax - b */
202: VecAYPX(R,-1.0,B);
203: KSP_PCApply(ksp,R,U); /* u <- Br */
204: KSP_MatMult(ksp,Amat,U,W); /* w <- Au */
205:
206: KSP_MatMult(ksp,Amat,P,S); /* s <- Ap */
207: KSP_PCApply(ksp,S,Q); /* q <- Bs */
208: KSP_MatMult(ksp,Amat,Q,Z); /* z <- Aq */
210: replace = 1;
211: totreplaces++;
212: }
213: }
215: i++;
216: ksp->its = i;
218: } while (i<ksp->max_it);
219: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
220: return(0);
221: }
224: /*MC
225: KSPPIPECGRR - Pipelined conjugate gradient method with automated residual replacements.
227: This method has only a single non-blocking reduction per iteration, compared to 2 blocking for standard CG. The
228: non-blocking reduction is overlapped by the matrix-vector product and preconditioner application.
230: KSPPIPECGRR improves the robustness of KSPPIPECG by adding an automated residual replacement strategy.
231: True residual and auxiliary variables are computed explicitly in a number of dynamically determined
232: iterations to reset accumulated rounding errors.
234: See also KSPPIPECG, which is identical to KSPPIPECGRR without residual replacements.
235: See also KSPPIPECR, where the reduction is only overlapped with the matrix-vector product.
237: Level: intermediate
239: Notes:
240: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for
241: performance of pipelined methods. See the FAQ on the PETSc website for details.
243: Contributed by:
244: Siegfried Cools, Universiteit Antwerpen,
245: EXA2CT European Project on EXascale Algorithms and Advanced Computational Techniques
247: Reference:
248: S. Cools, E.F. Yetkin, E. Agullo, L. Giraud, W. Vanroose, "Analysis of rounding error accumulation in the
249: conjugate gradients method to improve the maximal attainable accuracy of pipelined CG".
250: Submitted to SIAM J. Sci. Comp., February 2016.
252: .seealso: KSPCreate(), KSPSetType(), KSPPIPECR, KSPGROPPCG, KSPPIPECG, KSPPGMRES, KSPCG, KSPCGUseSingleReduction()
253: M*/
256: PETSC_EXTERN PetscErrorCode KSPCreate_PIPECGRR(KSP ksp)
257: {
261: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
262: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
263: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
265: ksp->ops->setup = KSPSetUp_PIPECGRR;
266: ksp->ops->solve = KSPSolve_PIPECGRR;
267: ksp->ops->destroy = KSPDestroyDefault;
268: ksp->ops->view = 0;
269: ksp->ops->setfromoptions = 0;
270: ksp->ops->buildsolution = KSPBuildSolutionDefault;
271: ksp->ops->buildresidual = KSPBuildResidualDefault;
272: return(0);
273: }