Actual source code: ibcgs.c


  2: #include <petsc/private/kspimpl.h>
  3: #include <petsc/private/vecimpl.h>

  5: static PetscErrorCode KSPSetUp_IBCGS(KSP ksp)
  6: {
  7:   PetscBool      diagonalscale;

  9:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
 11:   KSPSetWorkVecs(ksp,9);
 12:   return 0;
 13: }

 15: /*
 16:     The code below "cheats" from PETSc style
 17:        1) VecRestoreArray() is called immediately after VecGetArray() and the array values are still accessed; the reason for the immediate
 18:           restore is that Vec operations are done on some of the vectors during the solve and if we did not restore immediately it would
 19:           generate two VecGetArray() (the second one inside the Vec operation) calls without a restore between them.
 20:        2) The vector operations on done directly on the arrays instead of with VecXXXX() calls

 22:        For clarity in the code we name single VECTORS with two names, for example, Rn_1 and R, but they actually always
 23:      the exact same memory. We do this with macro defines so that compiler won't think they are
 24:      two different variables.

 26: */
 27: #define Xn_1 Xn
 28: #define xn_1 xn
 29: #define Rn_1 Rn
 30: #define rn_1 rn
 31: #define Un_1 Un
 32: #define un_1 un
 33: #define Vn_1 Vn
 34: #define vn_1 vn
 35: #define Qn_1 Qn
 36: #define qn_1 qn
 37: #define Zn_1 Zn
 38: #define zn_1 zn
 39: static PetscErrorCode  KSPSolve_IBCGS(KSP ksp)
 40: {
 41:   PetscInt       i,N;
 42:   PetscReal      rnorm = 0.0,rnormin = 0.0;
 43: #if defined(PETSC_HAVE_MPI_LONG_DOUBLE) && !defined(PETSC_USE_COMPLEX) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
 44:   /* Because of possible instabilities in the algorithm (as indicated by different residual histories for the same problem
 45:      on the same number of processes  with different runs) we support computing the inner products using Intel's 80 bit arithmetic
 46:      rather than just 64 bit. Thus we copy our double precision values into long doubles (hoping this keeps the 16 extra bits)
 47:      and tell MPI to do its ALlreduces with MPI_LONG_DOUBLE.

 49:      Note for developers that does not effect the code. Intel's long double is implemented by storing the 80 bits of extended double
 50:      precision into a 16 byte space (the rest of the space is ignored)  */
 51:   long double insums[7],outsums[7];
 52: #else
 53:   PetscScalar insums[7],outsums[7];
 54: #endif
 55:   PetscScalar                       sigman_2, sigman_1, sigman, pin_1, pin, phin_1, phin,tmp1,tmp2;
 56:   PetscScalar                       taun_1, taun, rhon, alphan_1, alphan, omegan_1, omegan;
 57:   const PetscScalar *PETSC_RESTRICT r0, *PETSC_RESTRICT f0, *PETSC_RESTRICT qn, *PETSC_RESTRICT b, *PETSC_RESTRICT un;
 58:   PetscScalar *PETSC_RESTRICT       rn, *PETSC_RESTRICT xn, *PETSC_RESTRICT vn, *PETSC_RESTRICT zn;
 59:   /* the rest do not have to keep n_1 values */
 60:   PetscScalar                       kappan, thetan, etan, gamman, betan, deltan;
 61:   const PetscScalar *PETSC_RESTRICT tn;
 62:   PetscScalar *PETSC_RESTRICT       sn;
 63:   Vec                               R0,Rn,Xn,F0,Vn,Zn,Qn,Tn,Sn,B,Un;
 64:   Mat                               A;


 68:  #if defined(PETSC_HAVE_MPI_LONG_DOUBLE) && !defined(PETSC_USE_COMPLEX) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
 69:   /* since 80 bit long doubls do not fill the upper bits, we fill them initially so that
 70:      valgrind won't detect MPI_Allreduce() with uninitialized data */
 71:   PetscMemzero(insums,sizeof(insums));
 72:   PetscMemzero(insums,sizeof(insums));
 73: #endif

 75:   PCGetOperators(ksp->pc,&A,NULL);
 76:   VecGetLocalSize(ksp->vec_sol,&N);
 77:   Xn   = ksp->vec_sol; VecGetArray(Xn_1,(PetscScalar**)&xn_1)); PetscCall(VecRestoreArray(Xn_1,NULL);
 78:   B    = ksp->vec_rhs; VecGetArrayRead(B,(const PetscScalar**)&b)); PetscCall(VecRestoreArrayRead(B,NULL);
 79:   R0   = ksp->work[0]; VecGetArrayRead(R0,(const PetscScalar**)&r0)); PetscCall(VecRestoreArrayRead(R0,NULL);
 80:   Rn   = ksp->work[1]; VecGetArray(Rn_1,(PetscScalar**)&rn_1)); PetscCall(VecRestoreArray(Rn_1,NULL);
 81:   Un   = ksp->work[2]; VecGetArrayRead(Un_1,(const PetscScalar**)&un_1)); PetscCall(VecRestoreArrayRead(Un_1,NULL);
 82:   F0   = ksp->work[3]; VecGetArrayRead(F0,(const PetscScalar**)&f0)); PetscCall(VecRestoreArrayRead(F0,NULL);
 83:   Vn   = ksp->work[4]; VecGetArray(Vn_1,(PetscScalar**)&vn_1)); PetscCall(VecRestoreArray(Vn_1,NULL);
 84:   Zn   = ksp->work[5]; VecGetArray(Zn_1,(PetscScalar**)&zn_1)); PetscCall(VecRestoreArray(Zn_1,NULL);
 85:   Qn   = ksp->work[6]; VecGetArrayRead(Qn_1,(const PetscScalar**)&qn_1)); PetscCall(VecRestoreArrayRead(Qn_1,NULL);
 86:   Tn   = ksp->work[7]; VecGetArrayRead(Tn,(const PetscScalar**)&tn)); PetscCall(VecRestoreArrayRead(Tn,NULL);
 87:   Sn   = ksp->work[8]; VecGetArrayRead(Sn,(const PetscScalar**)&sn)); PetscCall(VecRestoreArrayRead(Sn,NULL);

 89:   /* r0 = rn_1 = b - A*xn_1; */
 90:   /* KSP_PCApplyBAorAB(ksp,Xn_1,Rn_1,Tn);
 91:      VecAYPX(Rn_1,-1.0,B); */
 92:   KSPInitialResidual(ksp,Xn_1,Tn,Sn,Rn_1,B);
 93:   if (ksp->normtype != KSP_NORM_NONE) {
 94:     VecNorm(Rn_1,NORM_2,&rnorm);
 95:     KSPCheckNorm(ksp,rnorm);
 96:   }
 97:   KSPMonitor(ksp,0,rnorm);
 98:   (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
 99:   if (ksp->reason) return 0;

101:   VecCopy(Rn_1,R0);

103:   /* un_1 = A*rn_1; */
104:   KSP_PCApplyBAorAB(ksp,Rn_1,Un_1,Tn);

106:   /* f0   = A'*rn_1; */
107:   if (ksp->pc_side == PC_RIGHT) { /* B' A' */
108:     KSP_MatMultTranspose(ksp,A,R0,Tn);
109:     KSP_PCApplyTranspose(ksp,Tn,F0);
110:   } else if (ksp->pc_side == PC_LEFT) { /* A' B' */
111:     KSP_PCApplyTranspose(ksp,R0,Tn);
112:     KSP_MatMultTranspose(ksp,A,Tn,F0);
113:   }

115:   /*qn_1 = vn_1 = zn_1 = 0.0; */
116:   VecSet(Qn_1,0.0);
117:   VecSet(Vn_1,0.0);
118:   VecSet(Zn_1,0.0);

120:   sigman_2 = pin_1 = taun_1 = 0.0;

122:   /* the paper says phin_1 should be initialized to zero, it is actually R0'R0 */
123:   VecDot(R0,R0,&phin_1);
124:   KSPCheckDot(ksp,phin_1);

126:   /* sigman_1 = rn_1'un_1  */
127:   VecDot(R0,Un_1,&sigman_1);

129:   alphan_1 = omegan_1 = 1.0;

131:   for (ksp->its = 1; ksp->its<ksp->max_it+1; ksp->its++) {
132:     rhon = phin_1 - omegan_1*sigman_2 + omegan_1*alphan_1*pin_1;
133:     if (ksp->its == 1) deltan = rhon;
134:     else deltan = rhon/taun_1;
135:     betan = deltan/omegan_1;
136:     taun  = sigman_1 + betan*taun_1  - deltan*pin_1;
137:     if (taun == 0.0) {
139:       else {
140:         ksp->reason = KSP_DIVERGED_NANORINF;
141:         return 0;
142:       }
143:     }
144:     alphan = rhon/taun;
145:     PetscLogFlops(15.0);

147:     /*
148:         zn = alphan*rn_1 + (alphan/alphan_1)betan*zn_1 - alphan*deltan*vn_1
149:         vn = un_1 + betan*vn_1 - deltan*qn_1
150:         sn = rn_1 - alphan*vn

152:        The algorithm in the paper is missing the alphan/alphan_1 term in the zn update
153:     */
154:     PetscLogEventBegin(VEC_Ops,0,0,0,0);
155:     tmp1 = (alphan/alphan_1)*betan;
156:     tmp2 = alphan*deltan;
157:     for (i=0; i<N; i++) {
158:       zn[i] = alphan*rn_1[i] + tmp1*zn_1[i] - tmp2*vn_1[i];
159:       vn[i] = un_1[i] + betan*vn_1[i] - deltan*qn_1[i];
160:       sn[i] = rn_1[i] - alphan*vn[i];
161:     }
162:     PetscLogFlops(3.0+11.0*N);
163:     PetscLogEventEnd(VEC_Ops,0,0,0,0);

165:     /*
166:         qn = A*vn
167:     */
168:     KSP_PCApplyBAorAB(ksp,Vn,Qn,Tn);

170:     /*
171:         tn = un_1 - alphan*qn
172:     */
173:     VecWAXPY(Tn,-alphan,Qn,Un_1);

175:     /*
176:         phin = r0'sn
177:         pin  = r0'qn
178:         gamman = f0'sn
179:         etan   = f0'tn
180:         thetan = sn'tn
181:         kappan = tn'tn
182:     */
183:     PetscLogEventBegin(VEC_ReduceArithmetic,0,0,0,0);
184:     phin = pin = gamman = etan = thetan = kappan = 0.0;
185:     for (i=0; i<N; i++) {
186:       phin   += r0[i]*sn[i];
187:       pin    += r0[i]*qn[i];
188:       gamman += f0[i]*sn[i];
189:       etan   += f0[i]*tn[i];
190:       thetan += sn[i]*tn[i];
191:       kappan += tn[i]*tn[i];
192:     }
193:     PetscLogFlops(12.0*N);
194:     PetscLogEventEnd(VEC_ReduceArithmetic,0,0,0,0);

196:     insums[0] = phin;
197:     insums[1] = pin;
198:     insums[2] = gamman;
199:     insums[3] = etan;
200:     insums[4] = thetan;
201:     insums[5] = kappan;
202:     insums[6] = rnormin;

204:     PetscLogEventBegin(VEC_ReduceCommunication,0,0,0,0);
205: #if defined(PETSC_HAVE_MPI_LONG_DOUBLE) && !defined(PETSC_USE_COMPLEX) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
206:     if (ksp->lagnorm && ksp->its > 1) {
207:       MPIU_Allreduce(insums,outsums,7,MPI_LONG_DOUBLE,MPI_SUM,PetscObjectComm((PetscObject)ksp));
208:     } else {
209:       MPIU_Allreduce(insums,outsums,6,MPI_LONG_DOUBLE,MPI_SUM,PetscObjectComm((PetscObject)ksp));
210:     }
211: #else
212:     if (ksp->lagnorm && ksp->its > 1 && ksp->normtype != KSP_NORM_NONE) {
213:       MPIU_Allreduce(insums,outsums,7,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
214:     } else {
215:       MPIU_Allreduce(insums,outsums,6,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
216:     }
217: #endif
218:     PetscLogEventEnd(VEC_ReduceCommunication,0,0,0,0);
219:     phin   = outsums[0];
220:     pin    = outsums[1];
221:     gamman = outsums[2];
222:     etan   = outsums[3];
223:     thetan = outsums[4];
224:     kappan = outsums[5];
225:     if (ksp->lagnorm && ksp->its > 1 && ksp->normtype != KSP_NORM_NONE) rnorm = PetscSqrtReal(PetscRealPart(outsums[6]));

227:     if (kappan == 0.0) {
229:       else {
230:         ksp->reason = KSP_DIVERGED_NANORINF;
231:         return 0;
232:       }
233:     }
234:     if (thetan == 0.0) {
236:       else {
237:         ksp->reason = KSP_DIVERGED_NANORINF;
238:         return 0;
239:       }
240:     }
241:     omegan = thetan/kappan;
242:     sigman = gamman - omegan*etan;

244:     /*
245:         rn = sn - omegan*tn
246:         xn = xn_1 + zn + omegan*sn
247:     */
248:     PetscLogEventBegin(VEC_Ops,0,0,0,0);
249:     rnormin = 0.0;
250:     for (i=0; i<N; i++) {
251:       rn[i]    = sn[i] - omegan*tn[i];
252:       rnormin += PetscRealPart(PetscConj(rn[i])*rn[i]);
253:       xn[i]   += zn[i] + omegan*sn[i];
254:     }
255:     PetscObjectStateIncrease((PetscObject)Xn);
256:     PetscLogFlops(7.0*N);
257:     PetscLogEventEnd(VEC_Ops,0,0,0,0);

259:     if (!ksp->lagnorm && ksp->chknorm < ksp->its && ksp->normtype != KSP_NORM_NONE) {
260:       PetscLogEventBegin(VEC_ReduceCommunication,0,0,0,0);
261:       MPIU_Allreduce(&rnormin,&rnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
262:       PetscLogEventEnd(VEC_ReduceCommunication,0,0,0,0);
263:       rnorm = PetscSqrtReal(rnorm);
264:     }

266:     /* Test for convergence */
267:     KSPMonitor(ksp,ksp->its,rnorm);
268:     (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
269:     if (ksp->reason) {
270:       KSPUnwindPreconditioner(ksp,Xn,Tn);
271:       return 0;
272:     }

274:     /* un = A*rn */
275:     KSP_PCApplyBAorAB(ksp,Rn,Un,Tn);

277:     /* Update n-1 locations with n locations */
278:     sigman_2 = sigman_1;
279:     sigman_1 = sigman;
280:     pin_1    = pin;
281:     phin_1   = phin;
282:     alphan_1 = alphan;
283:     taun_1   = taun;
284:     omegan_1 = omegan;
285:   }
286:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
287:   KSPUnwindPreconditioner(ksp,Xn,Tn);
288:   return 0;
289: }

291: /*MC
292:      KSPIBCGS - Implements the IBiCGStab (Improved Stabilized version of BiConjugate Gradient) method
293:             in an alternative form to have only a single global reduction operation instead of the usual 3 (or 4)

295:    Options Database Keys:
296:     see KSPSolve()

298:    Level: beginner

300:    Notes:
301:     Supports left and right preconditioning

303:           See KSPBCGSL for additional stabilization

305:           Unlike the Bi-CG-stab algorithm, this requires one multiplication be the transpose of the operator
306:            before the iteration starts.

308:           The paper has two errors in the algorithm presented, they are fixed in the code in KSPSolve_IBCGS()

310:           For maximum reduction in the number of global reduction operations, this solver should be used with
311:           KSPSetLagNorm().

313:           This is not supported for complex numbers.

315:    Reference: The Improved BiCGStab Method for Large and Sparse Unsymmetric Linear Systems on Parallel Distributed Memory
316:                      Architectures. L. T. Yang and R. Brent, Proceedings of the Fifth International Conference on Algorithms and
317:                      Architectures for Parallel Processing, 2002, IEEE.

319: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPBCGSL, KSPIBCGS, KSPSetLagNorm()
320: M*/

322: PETSC_EXTERN PetscErrorCode KSPCreate_IBCGS(KSP ksp)
323: {

325:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
326:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
327:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

329:   ksp->ops->setup          = KSPSetUp_IBCGS;
330:   ksp->ops->solve          = KSPSolve_IBCGS;
331:   ksp->ops->destroy        = KSPDestroyDefault;
332:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
333:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
334:   ksp->ops->setfromoptions = NULL;
335:   ksp->ops->view           = NULL;
336: #if defined(PETSC_USE_COMPLEX)
337:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"This is not supported for complex numbers");
338: #else
339:   return 0;
340: #endif
341: }