Actual source code: pbjacobi.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  2: /*
  3:    Include files needed for the PBJacobi preconditioner:
  4:      pcimpl.h - private include file intended for use by all preconditioners
  5: */

  7:  #include <petsc/private/pcimpl.h>

  9: /*
 10:    Private context (data structure) for the PBJacobi preconditioner.
 11: */
 12: typedef struct {
 13:   const MatScalar *diag;
 14:   PetscInt        bs,mbs;
 15: } PC_PBJacobi;


 18: static PetscErrorCode PCApply_PBJacobi_1(PC pc,Vec x,Vec y)
 19: {
 20:   PC_PBJacobi       *jac = (PC_PBJacobi*)pc->data;
 21:   PetscErrorCode    ierr;
 22:   PetscInt          i,m = jac->mbs;
 23:   const MatScalar   *diag = jac->diag;
 24:   const PetscScalar *xx;
 25:   PetscScalar       *yy;

 28:   VecGetArrayRead(x,&xx);
 29:   VecGetArray(y,&yy);
 30:   for (i=0; i<m; i++) yy[i] = diag[i]*xx[i];
 31:   VecRestoreArrayRead(x,&xx);
 32:   VecRestoreArray(y,&yy);
 33:   PetscLogFlops(m);
 34:   return(0);
 35: }

 37: static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y)
 38: {
 39:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
 40:   PetscErrorCode  ierr;
 41:   PetscInt        i,m = jac->mbs;
 42:   const MatScalar *diag = jac->diag;
 43:   PetscScalar     x0,x1,*yy;
 44:   const PetscScalar *xx;

 47:   VecGetArrayRead(x,&xx);
 48:   VecGetArray(y,&yy);
 49:   for (i=0; i<m; i++) {
 50:     x0        = xx[2*i]; x1 = xx[2*i+1];
 51:     yy[2*i]   = diag[0]*x0 + diag[2]*x1;
 52:     yy[2*i+1] = diag[1]*x0 + diag[3]*x1;
 53:     diag     += 4;
 54:   }
 55:   VecRestoreArrayRead(x,&xx);
 56:   VecRestoreArray(y,&yy);
 57:   PetscLogFlops(6.0*m);
 58:   return(0);
 59: }
 60: static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y)
 61: {
 62:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
 63:   PetscErrorCode  ierr;
 64:   PetscInt        i,m = jac->mbs;
 65:   const MatScalar *diag = jac->diag;
 66:   PetscScalar     x0,x1,x2,*yy;
 67:   const PetscScalar *xx;

 70:   VecGetArrayRead(x,&xx);
 71:   VecGetArray(y,&yy);
 72:   for (i=0; i<m; i++) {
 73:     x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2];

 75:     yy[3*i]   = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
 76:     yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
 77:     yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
 78:     diag     += 9;
 79:   }
 80:   VecRestoreArrayRead(x,&xx);
 81:   VecRestoreArray(y,&yy);
 82:   PetscLogFlops(15.0*m);
 83:   return(0);
 84: }
 85: static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y)
 86: {
 87:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
 88:   PetscErrorCode  ierr;
 89:   PetscInt        i,m = jac->mbs;
 90:   const MatScalar *diag = jac->diag;
 91:   PetscScalar     x0,x1,x2,x3,*yy;
 92:   const PetscScalar *xx;

 95:   VecGetArrayRead(x,&xx);
 96:   VecGetArray(y,&yy);
 97:   for (i=0; i<m; i++) {
 98:     x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3];

100:     yy[4*i]   = diag[0]*x0 + diag[4]*x1 + diag[8]*x2  + diag[12]*x3;
101:     yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2  + diag[13]*x3;
102:     yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
103:     yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
104:     diag     += 16;
105:   }
106:   VecRestoreArrayRead(x,&xx);
107:   VecRestoreArray(y,&yy);
108:   PetscLogFlops(28.0*m);
109:   return(0);
110: }
111: static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y)
112: {
113:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
114:   PetscErrorCode  ierr;
115:   PetscInt        i,m = jac->mbs;
116:   const MatScalar *diag = jac->diag;
117:   PetscScalar     x0,x1,x2,x3,x4,*yy;
118:   const PetscScalar *xx;

121:   VecGetArrayRead(x,&xx);
122:   VecGetArray(y,&yy);
123:   for (i=0; i<m; i++) {
124:     x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4];

126:     yy[5*i]   = diag[0]*x0 + diag[5]*x1 + diag[10]*x2  + diag[15]*x3 + diag[20]*x4;
127:     yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2  + diag[16]*x3 + diag[21]*x4;
128:     yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
129:     yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
130:     yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
131:     diag     += 25;
132:   }
133:   VecRestoreArrayRead(x,&xx);
134:   VecRestoreArray(y,&yy);
135:   PetscLogFlops(45.0*m);
136:   return(0);
137: }
138: static PetscErrorCode PCApply_PBJacobi_6(PC pc,Vec x,Vec y)
139: {
140:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
141:   PetscErrorCode  ierr;
142:   PetscInt        i,m = jac->mbs;
143:   const MatScalar *diag = jac->diag;
144:   PetscScalar     x0,x1,x2,x3,x4,x5,*yy;
145:   const PetscScalar *xx;

148:   VecGetArrayRead(x,&xx);
149:   VecGetArray(y,&yy);
150:   for (i=0; i<m; i++) {
151:     x0 = xx[6*i]; x1 = xx[6*i+1]; x2 = xx[6*i+2]; x3 = xx[6*i+3]; x4 = xx[6*i+4]; x5 = xx[6*i+5];

153:     yy[6*i]   = diag[0]*x0 + diag[6]*x1  + diag[12]*x2  + diag[18]*x3 + diag[24]*x4 + diag[30]*x5;
154:     yy[6*i+1] = diag[1]*x0 + diag[7]*x1  + diag[13]*x2  + diag[19]*x3 + diag[25]*x4 + diag[31]*x5;
155:     yy[6*i+2] = diag[2]*x0 + diag[8]*x1  + diag[14]*x2  + diag[20]*x3 + diag[26]*x4 + diag[32]*x5;
156:     yy[6*i+3] = diag[3]*x0 + diag[9]*x1  + diag[15]*x2  + diag[21]*x3 + diag[27]*x4 + diag[33]*x5;
157:     yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2  + diag[22]*x3 + diag[28]*x4 + diag[34]*x5;
158:     yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2  + diag[23]*x3 + diag[29]*x4 + diag[35]*x5;
159:     diag     += 36;
160:   }
161:   VecRestoreArrayRead(x,&xx);
162:   VecRestoreArray(y,&yy);
163:   PetscLogFlops(66.0*m);
164:   return(0);
165: }
166: static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y)
167: {
168:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
169:   PetscErrorCode  ierr;
170:   PetscInt        i,m = jac->mbs;
171:   const MatScalar *diag = jac->diag;
172:   PetscScalar     x0,x1,x2,x3,x4,x5,x6,*yy;
173:   const PetscScalar *xx;

176:   VecGetArrayRead(x,&xx);
177:   VecGetArray(y,&yy);
178:   for (i=0; i<m; i++) {
179:     x0 = xx[7*i]; x1 = xx[7*i+1]; x2 = xx[7*i+2]; x3 = xx[7*i+3]; x4 = xx[7*i+4]; x5 = xx[7*i+5]; x6 = xx[7*i+6];

181:     yy[7*i]   = diag[0]*x0 + diag[7]*x1  + diag[14]*x2  + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
182:     yy[7*i+1] = diag[1]*x0 + diag[8]*x1  + diag[15]*x2  + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6;
183:     yy[7*i+2] = diag[2]*x0 + diag[9]*x1  + diag[16]*x2  + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6;
184:     yy[7*i+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2  + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6;
185:     yy[7*i+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2  + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6;
186:     yy[7*i+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2  + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6;
187:     yy[7*i+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2  + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6;
188:     diag     += 49;
189:   }
190:   VecRestoreArrayRead(x,&xx);
191:   VecRestoreArray(y,&yy);
192:   PetscLogFlops(91.0*m); /* 2*bs2 - bs */
193:   return(0);
194: }
195: static PetscErrorCode PCApply_PBJacobi_N(PC pc,Vec x,Vec y)
196: {
197:   PC_PBJacobi       *jac = (PC_PBJacobi*)pc->data;
198:   PetscErrorCode    ierr;
199:   PetscInt          i,ib,jb;
200:   const PetscInt    m = jac->mbs;
201:   const PetscInt    bs = jac->bs;
202:   const MatScalar   *diag = jac->diag;
203:   PetscScalar       *yy;
204:   const PetscScalar *xx;

207:   VecGetArrayRead(x,&xx);
208:   VecGetArray(y,&yy);
209:   for (i=0; i<m; i++) {
210:     for (ib=0; ib<bs; ib++){
211:       PetscScalar rowsum = 0;
212:       for (jb=0; jb<bs; jb++){
213:         rowsum += diag[ib+jb*bs] * xx[bs*i+jb];
214:       }
215:       yy[bs*i+ib] = rowsum;
216:     }
217:     diag += bs*bs;
218:   }
219:   VecRestoreArrayRead(x,&xx);
220:   VecRestoreArray(y,&yy);
221:   PetscLogFlops((2.0*bs*bs-bs)*m); /* 2*bs2 - bs */
222:   return(0);
223: }
224: /* -------------------------------------------------------------------------- */
225: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
226: {
227:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
229:   Mat            A = pc->pmat;
230:   MatFactorError err;
231:   PetscInt       nlocal;
232: 
234:   MatInvertBlockDiagonal(A,&jac->diag);
235:   MatFactorGetError(A,&err);
236:   if (err) pc->failedreason = (PCFailedReason)err;
237: 
238:   MatGetBlockSize(A,&jac->bs);
239:   MatGetLocalSize(A,&nlocal,NULL);
240:   jac->mbs = nlocal/jac->bs;
241:   switch (jac->bs) {
242:   case 1:
243:     pc->ops->apply = PCApply_PBJacobi_1;
244:     break;
245:   case 2:
246:     pc->ops->apply = PCApply_PBJacobi_2;
247:     break;
248:   case 3:
249:     pc->ops->apply = PCApply_PBJacobi_3;
250:     break;
251:   case 4:
252:     pc->ops->apply = PCApply_PBJacobi_4;
253:     break;
254:   case 5:
255:     pc->ops->apply = PCApply_PBJacobi_5;
256:     break;
257:   case 6:
258:     pc->ops->apply = PCApply_PBJacobi_6;
259:     break;
260:   case 7:
261:     pc->ops->apply = PCApply_PBJacobi_7;
262:     break;
263:   default:
264:     pc->ops->apply = PCApply_PBJacobi_N;
265:     break;
266:   }
267:   return(0);
268: }
269: /* -------------------------------------------------------------------------- */
270: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
271: {

275:   /*
276:       Free the private data structure that was hanging off the PC
277:   */
278:   PetscFree(pc->data);
279:   return(0);
280: }

282: static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer)
283: {
285:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
286:   PetscBool      iascii;

289:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
290:   if (iascii) {
291:     PetscViewerASCIIPrintf(viewer,"  point-block size %D\n",jac->bs);
292:   }
293:   return(0);
294: }

296: /* -------------------------------------------------------------------------- */
297: /*MC
298:      PCPBJACOBI - Point block Jacobi preconditioner


301:    Notes:
302:     See PCJACOBI for point Jacobi preconditioning

304:    This works for AIJ and BAIJ matrices and uses the blocksize provided to the matrix

306:    Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot
307:    is detected a PETSc error is generated.

309:    Developer Notes:
310:     This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow
311:    the factorization to continue even after a zero pivot is found resulting in a Nan and hence
312:    terminating KSP with a KSP_DIVERGED_NANORIF allowing
313:    a nonlinear solver/ODE integrator to recover without stopping the program as currently happens.

315:    Developer Note: Perhaps should provide an option that allows generation of a valid preconditioner
316:    even if a block is singular as the PCJACOBI does.

318:    Level: beginner

320:   Concepts: point block Jacobi


323: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI

325: M*/

327: PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc)
328: {
329:   PC_PBJacobi    *jac;

333:   /*
334:      Creates the private data structure for this preconditioner and
335:      attach it to the PC object.
336:   */
337:   PetscNewLog(pc,&jac);
338:   pc->data = (void*)jac;

340:   /*
341:      Initialize the pointers to vectors to ZERO; these will be used to store
342:      diagonal entries of the matrix for fast preconditioner application.
343:   */
344:   jac->diag = 0;

346:   /*
347:       Set the pointers for the functions that are provided above.
348:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
349:       are called, they will automatically call these functions.  Note we
350:       choose not to provide a couple of these functions since they are
351:       not needed.
352:   */
353:   pc->ops->apply               = 0; /*set depending on the block size */
354:   pc->ops->applytranspose      = 0;
355:   pc->ops->setup               = PCSetUp_PBJacobi;
356:   pc->ops->destroy             = PCDestroy_PBJacobi;
357:   pc->ops->setfromoptions      = 0;
358:   pc->ops->view                = PCView_PBJacobi;
359:   pc->ops->applyrichardson     = 0;
360:   pc->ops->applysymmetricleft  = 0;
361:   pc->ops->applysymmetricright = 0;
362:   return(0);
363: }