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