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