1 #include <petsc/private/kspimpl.h>
2
3 /*
4 KSPSetUp_GROPPCG - Sets up the workspace needed by the GROPPCG method.
5
6 This is called once, usually automatically by KSPSolve() or KSPSetUp()
7 but can be called directly by KSPSetUp()
8 */
KSPSetUp_GROPPCG(KSP ksp)9 static PetscErrorCode KSPSetUp_GROPPCG(KSP ksp)
10 {
11 PetscErrorCode ierr;
12
13 PetscFunctionBegin;
14 ierr = KSPSetWorkVecs(ksp,6);CHKERRQ(ierr);
15 PetscFunctionReturn(0);
16 }
17
18 /*
19 KSPSolve_GROPPCG
20
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 */
KSPSolve_GROPPCG(KSP ksp)25 static PetscErrorCode KSPSolve_GROPPCG(KSP ksp)
26 {
27 PetscErrorCode ierr;
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;
34
35 PetscFunctionBegin;
36 ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
37 if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
38
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];
47
48 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
49
50 ksp->its = 0;
51 if (!ksp->guess_zero) {
52 ierr = KSP_MatMult(ksp,Amat,x,r);CHKERRQ(ierr); /* r <- b - Ax */
53 ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
54 } else {
55 ierr = VecCopy(b,r);CHKERRQ(ierr); /* r <- b (x is 0) */
56 }
57
58 ierr = KSP_PCApply(ksp,r,z);CHKERRQ(ierr); /* z <- Br */
59 ierr = VecCopy(z,p);CHKERRQ(ierr); /* p <- z */
60 ierr = VecDotBegin(r,z,&gamma);CHKERRQ(ierr); /* gamma <- z'*r */
61 ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));CHKERRQ(ierr);
62 ierr = KSP_MatMult(ksp,Amat,p,s);CHKERRQ(ierr); /* s <- Ap */
63 ierr = VecDotEnd(r,z,&gamma);CHKERRQ(ierr); /* gamma <- z'*r */
64
65 switch (ksp->normtype) {
66 case KSP_NORM_PRECONDITIONED:
67 /* This could be merged with the computation of gamma above */
68 ierr = VecNorm(z,NORM_2,&dp);CHKERRQ(ierr); /* 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 ierr = VecNorm(r,NORM_2,&dp);CHKERRQ(ierr); /* 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 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
84 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
85 ksp->rnorm = dp;
86 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */
87 if (ksp->reason) PetscFunctionReturn(0);
88
89 i = 0;
90 do {
91 ksp->its = i+1;
92 i++;
93
94 ierr = VecDotBegin(p,s,&t);CHKERRQ(ierr);
95 ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));CHKERRQ(ierr);
96
97 ierr = KSP_PCApply(ksp,s,S);CHKERRQ(ierr); /* S <- Bs */
98
99 ierr = VecDotEnd(p,s,&t);CHKERRQ(ierr);
100
101 alpha = gamma / t;
102 ierr = VecAXPY(x, alpha,p);CHKERRQ(ierr); /* x <- x + alpha * p */
103 ierr = VecAXPY(r,-alpha,s);CHKERRQ(ierr); /* r <- r - alpha * s */
104 ierr = VecAXPY(z,-alpha,S);CHKERRQ(ierr); /* z <- z - alpha * S */
105
106 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
107 ierr = VecNormBegin(r,NORM_2,&dp);CHKERRQ(ierr);
108 } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
109 ierr = VecNormBegin(z,NORM_2,&dp);CHKERRQ(ierr);
110 }
111 ierr = VecDotBegin(r,z,&gammaNew);CHKERRQ(ierr);
112 ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));CHKERRQ(ierr);
113
114 ierr = KSP_MatMult(ksp,Amat,z,Z);CHKERRQ(ierr); /* Z <- Az */
115
116 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
117 ierr = VecNormEnd(r,NORM_2,&dp);CHKERRQ(ierr);
118 } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
119 ierr = VecNormEnd(z,NORM_2,&dp);CHKERRQ(ierr);
120 }
121 ierr = VecDotEnd(r,z,&gammaNew);CHKERRQ(ierr);
122
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 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
131 ierr = KSPMonitor(ksp,i,dp);CHKERRQ(ierr);
132 ierr = (*ksp->converged)(ksp,i,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
133 if (ksp->reason) PetscFunctionReturn(0);
134
135 beta = gammaNew / gamma;
136 gamma = gammaNew;
137 ierr = VecAYPX(p,beta,z);CHKERRQ(ierr); /* p <- z + beta * p */
138 ierr = VecAYPX(s,beta,Z);CHKERRQ(ierr); /* s <- Z + beta * s */
139
140 } while (i<ksp->max_it);
141
142 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
143 PetscFunctionReturn(0);
144 }
145
146 PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP,Vec,Vec,Vec*);
147
148 /*MC
149 KSPGROPPCG - A pipelined conjugate gradient method from Bill Gropp
150
151 This method has two reductions, one of which is overlapped with the matrix-vector product and one of which is
152 overlapped with the preconditioner.
153
154 See also KSPPIPECG, which has only a single reduction that overlaps both the matrix-vector product and the preconditioner.
155
156 Level: intermediate
157
158 Notes:
159 MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
160 See the FAQ on the PETSc website for details.
161
162 Contributed by:
163 Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders
164
165 Reference:
166 http://www.cs.uiuc.edu/~wgropp/bib/talks/tdata/2012/icerm.pdf
167
168 .seealso: KSPCreate(), KSPSetType(), KSPPIPECG, KSPPIPECR, KSPPGMRES, KSPCG, KSPCGUseSingleReduction()
169 M*/
170
KSPCreate_GROPPCG(KSP ksp)171 PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
172 {
173 PetscErrorCode ierr;
174
175 PetscFunctionBegin;
176 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
177 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
178 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr);
179 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
180
181 ksp->ops->setup = KSPSetUp_GROPPCG;
182 ksp->ops->solve = KSPSolve_GROPPCG;
183 ksp->ops->destroy = KSPDestroyDefault;
184 ksp->ops->view = NULL;
185 ksp->ops->setfromoptions = NULL;
186 ksp->ops->buildsolution = KSPBuildSolutionDefault;
187 ksp->ops->buildresidual = KSPBuildResidual_CG;
188 PetscFunctionReturn(0);
189 }
190