Actual source code: groppcg.c
petsc-3.9.4 2018-09-11
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: {
14: KSPSetWorkVecs(ksp,6);
15: return(0);
16: }
18: /*
19: KSPSolve_GROPPCG
21: Input Parameter:
22: . ksp - the Krylov space object that was set to use conjugate gradient, by, for
23: example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
24: */
25: static PetscErrorCode KSPSolve_GROPPCG(KSP ksp)
26: {
28: PetscInt i;
29: PetscScalar alpha,beta = 0.0,gamma,gammaNew,t;
30: PetscReal dp = 0.0;
31: Vec x,b,r,p,s,S,z,Z;
32: Mat Amat,Pmat;
33: PetscBool diagonalscale;
36: PCGetDiagonalScale(ksp->pc,&diagonalscale);
37: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
39: x = ksp->vec_sol;
40: b = ksp->vec_rhs;
41: r = ksp->work[0];
42: p = ksp->work[1];
43: s = ksp->work[2];
44: S = ksp->work[3];
45: z = ksp->work[4];
46: Z = ksp->work[5];
48: PCGetOperators(ksp->pc,&Amat,&Pmat);
50: ksp->its = 0;
51: if (!ksp->guess_zero) {
52: KSP_MatMult(ksp,Amat,x,r); /* r <- b - Ax */
53: VecAYPX(r,-1.0,b);
54: } else {
55: VecCopy(b,r); /* r <- b (x is 0) */
56: }
58: KSP_PCApply(ksp,r,z); /* z <- Br */
59: VecCopy(z,p); /* p <- z */
60: VecDotBegin(r,z,&gamma); /* gamma <- z'*r */
61: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
62: KSP_MatMult(ksp,Amat,p,s); /* s <- Ap */
63: VecDotEnd(r,z,&gamma); /* gamma <- z'*r */
65: switch (ksp->normtype) {
66: case KSP_NORM_PRECONDITIONED:
67: /* This could be merged with the computation of gamma above */
68: VecNorm(z,NORM_2,&dp); /* dp <- z'*z = e'*A'*B'*B*A'*e' */
69: break;
70: case KSP_NORM_UNPRECONDITIONED:
71: /* This could be merged with the computation of gamma above */
72: VecNorm(r,NORM_2,&dp); /* dp <- r'*r = e'*A'*A*e */
73: break;
74: case KSP_NORM_NATURAL:
75: KSPCheckDot(ksp,gamma);
76: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
77: break;
78: case KSP_NORM_NONE:
79: dp = 0.0;
80: break;
81: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
82: }
83: KSPLogResidualHistory(ksp,dp);
84: KSPMonitor(ksp,0,dp);
85: ksp->rnorm = dp;
86: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP); /* test for convergence */
87: if (ksp->reason) return(0);
89: i = 0;
90: do {
91: ksp->its = i+1;
92: i++;
94: VecDotBegin(p,s,&t);
95: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));
97: KSP_PCApply(ksp,s,S); /* S <- Bs */
99: VecDotEnd(p,s,&t);
101: alpha = gamma / t;
102: VecAXPY(x, alpha,p); /* x <- x + alpha * p */
103: VecAXPY(r,-alpha,s); /* r <- r - alpha * s */
104: VecAXPY(z,-alpha,S); /* z <- z - alpha * S */
106: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
107: VecNormBegin(r,NORM_2,&dp);
108: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
109: VecNormBegin(z,NORM_2,&dp);
110: }
111: VecDotBegin(r,z,&gammaNew);
112: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
114: KSP_MatMult(ksp,Amat,z,Z); /* Z <- Az */
116: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
117: VecNormEnd(r,NORM_2,&dp);
118: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
119: VecNormEnd(z,NORM_2,&dp);
120: }
121: VecDotEnd(r,z,&gammaNew);
123: if (ksp->normtype == KSP_NORM_NATURAL) {
124: KSPCheckDot(ksp,gammaNew);
125: dp = PetscSqrtReal(PetscAbsScalar(gammaNew)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
126: } else if (ksp->normtype == KSP_NORM_NONE) {
127: dp = 0.0;
128: }
129: ksp->rnorm = dp;
130: KSPLogResidualHistory(ksp,dp);
131: KSPMonitor(ksp,i,dp);
132: (*ksp->converged)(ksp,i,dp,&ksp->reason,ksp->cnvP);
133: if (ksp->reason) break;
135: beta = gammaNew / gamma;
136: gamma = gammaNew;
137: VecAYPX(p,beta,z); /* p <- z + beta * p */
138: VecAYPX(s,beta,Z); /* s <- Z + beta * s */
140: } while (i<ksp->max_it);
142: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
143: return(0);
144: }
146: /*MC
147: KSPGROPPCG - A pipelined conjugate gradient method from Bill Gropp
149: This method has two reductions, one of which is overlapped with the matrix-vector product and one of which is
150: overlapped with the preconditioner.
152: See also KSPPIPECG, which has only a single reduction that overlaps both the matrix-vector product and the preconditioner.
154: Level: intermediate
156: Notes:
157: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
158: See the FAQ on the PETSc website for details.
160: Contributed by:
161: Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders
163: Reference:
164: http://www.cs.uiuc.edu/~wgropp/bib/talks/tdata/2012/icerm.pdf
166: .seealso: KSPCreate(), KSPSetType(), KSPPIPECG, KSPPIPECR, KSPPGMRES, KSPCG, KSPCGUseSingleReduction()
167: M*/
169: PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
170: {
174: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
175: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
176: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
177: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
179: ksp->ops->setup = KSPSetUp_GROPPCG;
180: ksp->ops->solve = KSPSolve_GROPPCG;
181: ksp->ops->destroy = KSPDestroyDefault;
182: ksp->ops->view = 0;
183: ksp->ops->setfromoptions = 0;
184: ksp->ops->buildsolution = KSPBuildSolutionDefault;
185: ksp->ops->buildresidual = KSPBuildResidualDefault;
186: return(0);
187: }