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: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
77: }
78: KSPLogResidualHistory(ksp,dp);
79: KSPMonitor(ksp,0,dp);
80: ksp->rnorm = dp;
81: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP); /* test for convergence */
82: if (ksp->reason) return 0;
84: i = 0;
85: do {
86: ksp->its = i+1;
87: i++;
89: VecDotBegin(p,s,&t);
90: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));
92: KSP_PCApply(ksp,s,S); /* S <- Bs */
94: VecDotEnd(p,s,&t);
96: alpha = gamma / t;
97: VecAXPY(x, alpha,p); /* x <- x + alpha * p */
98: VecAXPY(r,-alpha,s); /* r <- r - alpha * s */
99: VecAXPY(z,-alpha,S); /* z <- z - alpha * S */
101: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
102: VecNormBegin(r,NORM_2,&dp);
103: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
104: VecNormBegin(z,NORM_2,&dp);
105: }
106: VecDotBegin(r,z,&gammaNew);
107: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
109: KSP_MatMult(ksp,Amat,z,Z); /* Z <- Az */
111: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
112: VecNormEnd(r,NORM_2,&dp);
113: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
114: VecNormEnd(z,NORM_2,&dp);
115: }
116: VecDotEnd(r,z,&gammaNew);
118: if (ksp->normtype == KSP_NORM_NATURAL) {
119: KSPCheckDot(ksp,gammaNew);
120: dp = PetscSqrtReal(PetscAbsScalar(gammaNew)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
121: } else if (ksp->normtype == KSP_NORM_NONE) {
122: dp = 0.0;
123: }
124: ksp->rnorm = dp;
125: KSPLogResidualHistory(ksp,dp);
126: KSPMonitor(ksp,i,dp);
127: (*ksp->converged)(ksp,i,dp,&ksp->reason,ksp->cnvP);
128: if (ksp->reason) return 0;
130: beta = gammaNew / gamma;
131: gamma = gammaNew;
132: VecAYPX(p,beta,z); /* p <- z + beta * p */
133: VecAYPX(s,beta,Z); /* s <- Z + beta * s */
135: } while (i<ksp->max_it);
137: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
138: return 0;
139: }
141: PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP,Vec,Vec,Vec*);
143: /*MC
144: KSPGROPPCG - A pipelined conjugate gradient method from Bill Gropp
146: This method has two reductions, one of which is overlapped with the matrix-vector product and one of which is
147: overlapped with the preconditioner.
149: See also KSPPIPECG, which has only a single reduction that overlaps both the matrix-vector product and the preconditioner.
151: Level: intermediate
153: Notes:
154: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
155: See the FAQ on the PETSc website for details.
157: Contributed by:
158: Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders
160: Reference:
161: http://www.cs.uiuc.edu/~wgropp/bib/talks/tdata/2012/icerm.pdf
163: .seealso: KSPCreate(), KSPSetType(), KSPPIPECG, KSPPIPECR, KSPPGMRES, KSPCG, KSPCGUseSingleReduction()
164: M*/
166: PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
167: {
168: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
169: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
170: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
171: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
173: ksp->ops->setup = KSPSetUp_GROPPCG;
174: ksp->ops->solve = KSPSolve_GROPPCG;
175: ksp->ops->destroy = KSPDestroyDefault;
176: ksp->ops->view = NULL;
177: ksp->ops->setfromoptions = NULL;
178: ksp->ops->buildsolution = KSPBuildSolutionDefault;
179: ksp->ops->buildresidual = KSPBuildResidual_CG;
180: return 0;
181: }