Actual source code: pbjacobi.c

petsc-master 2016-08-23
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>   /*I "petscpc.h" I*/

  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;


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

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

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

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

 76:   VecGetArrayRead(x,&xx);
 77:   VecGetArray(y,&yy);
 78:   for (i=0; i<m; i++) {
 79:     x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2];

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

103:   VecGetArrayRead(x,&xx);
104:   VecGetArray(y,&yy);
105:   for (i=0; i<m; i++) {
106:     x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3];

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

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

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

160:   VecGetArrayRead(x,&xx);
161:   VecGetArray(y,&yy);
162:   for (i=0; i<m; i++) {
163:     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];

165:     yy[6*i]   = diag[0]*x0 + diag[6]*x1  + diag[12]*x2  + diag[18]*x3 + diag[24]*x4 + diag[30]*x5;
166:     yy[6*i+1] = diag[1]*x0 + diag[7]*x1  + diag[13]*x2  + diag[19]*x3 + diag[25]*x4 + diag[31]*x5;
167:     yy[6*i+2] = diag[2]*x0 + diag[8]*x1  + diag[14]*x2  + diag[20]*x3 + diag[26]*x4 + diag[32]*x5;
168:     yy[6*i+3] = diag[3]*x0 + diag[9]*x1  + diag[15]*x2  + diag[21]*x3 + diag[27]*x4 + diag[33]*x5;
169:     yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2  + diag[22]*x3 + diag[28]*x4 + diag[34]*x5;
170:     yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2  + diag[23]*x3 + diag[29]*x4 + diag[35]*x5;
171:     diag     += 36;
172:   }
173:   VecRestoreArrayRead(x,&xx);
174:   VecRestoreArray(y,&yy);
175:   PetscLogFlops(66.0*m);
176:   return(0);
177: }
180: static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y)
181: {
182:   PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
183:   PetscErrorCode  ierr;
184:   PetscInt        i,m = jac->mbs;
185:   const MatScalar *diag = jac->diag;
186:   PetscScalar     x0,x1,x2,x3,x4,x5,x6,*yy;
187:   const PetscScalar *xx;

190:   VecGetArrayRead(x,&xx);
191:   VecGetArray(y,&yy);
192:   for (i=0; i<m; i++) {
193:     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];

195:     yy[7*i]   = diag[0]*x0 + diag[7]*x1  + diag[14]*x2  + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
196:     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;
197:     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;
198:     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;
199:     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;
200:     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;
201:     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;
202:     diag     += 49;
203:   }
204:   VecRestoreArrayRead(x,&xx);
205:   VecRestoreArray(y,&yy);
206:   PetscLogFlops(80.0*m);
207:   return(0);
208: }
209: /* -------------------------------------------------------------------------- */
212: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
213: {
214:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
216:   Mat            A = pc->pmat;
217:   MatFactorError err;
218:   PetscInt       nlocal;
219: 
221:   MatInvertBlockDiagonal(A,&jac->diag);
222:   MatFactorGetError(A,&err);
223:   if (err) pc->failedreason = (PCFailedReason)err;
224: 
225:   MatGetBlockSize(A,&jac->bs);
226:   MatGetLocalSize(A,&nlocal,NULL);
227:   jac->mbs = nlocal/jac->bs;
228:   switch (jac->bs) {
229:   case 1:
230:     pc->ops->apply = PCApply_PBJacobi_1;
231:     break;
232:   case 2:
233:     pc->ops->apply = PCApply_PBJacobi_2;
234:     break;
235:   case 3:
236:     pc->ops->apply = PCApply_PBJacobi_3;
237:     break;
238:   case 4:
239:     pc->ops->apply = PCApply_PBJacobi_4;
240:     break;
241:   case 5:
242:     pc->ops->apply = PCApply_PBJacobi_5;
243:     break;
244:   case 6:
245:     pc->ops->apply = PCApply_PBJacobi_6;
246:     break;
247:   case 7:
248:     pc->ops->apply = PCApply_PBJacobi_7;
249:     break;
250:   default:
251:     SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
252:   }
253:   return(0);
254: }
255: /* -------------------------------------------------------------------------- */
258: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
259: {

263:   /*
264:       Free the private data structure that was hanging off the PC
265:   */
266:   PetscFree(pc->data);
267:   return(0);
268: }

272: static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer)
273: {
275:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
276:   PetscBool      iascii;

279:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
280:   if (iascii) {
281:     PetscViewerASCIIPrintf(viewer,"  point-block Jacobi: block size %D\n",jac->bs);
282:   }
283:   return(0);
284: }

286: /* -------------------------------------------------------------------------- */
287: /*MC
288:      PCPBJACOBI - Point block Jacobi preconditioner


291:    Notes: See PCJACOBI for point Jacobi preconditioning

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

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

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

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

306:    Level: beginner

308:   Concepts: point block Jacobi


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

313: M*/

317: PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc)
318: {
319:   PC_PBJacobi    *jac;

323:   /*
324:      Creates the private data structure for this preconditioner and
325:      attach it to the PC object.
326:   */
327:   PetscNewLog(pc,&jac);
328:   pc->data = (void*)jac;

330:   /*
331:      Initialize the pointers to vectors to ZERO; these will be used to store
332:      diagonal entries of the matrix for fast preconditioner application.
333:   */
334:   jac->diag = 0;

336:   /*
337:       Set the pointers for the functions that are provided above.
338:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
339:       are called, they will automatically call these functions.  Note we
340:       choose not to provide a couple of these functions since they are
341:       not needed.
342:   */
343:   pc->ops->apply               = 0; /*set depending on the block size */
344:   pc->ops->applytranspose      = 0;
345:   pc->ops->setup               = PCSetUp_PBJacobi;
346:   pc->ops->destroy             = PCDestroy_PBJacobi;
347:   pc->ops->setfromoptions      = 0;
348:   pc->ops->view                = PCView_PBJacobi;
349:   pc->ops->applyrichardson     = 0;
350:   pc->ops->applysymmetricleft  = 0;
351:   pc->ops->applysymmetricright = 0;
352:   return(0);
353: }