Actual source code: pbjacobi.c

petsc-3.4.5 2014-06-29
  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/matimpl.h>
  8: #include <petsc-private/pcimpl.h>   /*I "petscpc.h" I*/

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


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

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

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

 51:   VecGetArray(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:   VecRestoreArray(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,*xx,*yy;

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

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

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

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

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

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

156:   VecGetArray(x,&xx);
157:   VecGetArray(y,&yy);
158:   for (i=0; i<m; i++) {
159:     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];

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

185:   VecGetArray(x,&xx);
186:   VecGetArray(y,&yy);
187:   for (i=0; i<m; i++) {
188:     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];

190:     yy[7*i]   = diag[0]*x0 + diag[7]*x1  + diag[14]*x2  + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
191:     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;
192:     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;
193:     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;
194:     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;
195:     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;
196:     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;
197:     diag     += 49;
198:   }
199:   VecRestoreArray(x,&xx);
200:   VecRestoreArray(y,&yy);
201:   PetscLogFlops(80.0*m);
202:   return(0);
203: }
204: /* -------------------------------------------------------------------------- */
207: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
208: {
209:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
211:   Mat            A = pc->pmat;

214:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage");

216:   MatInvertBlockDiagonal(A,&jac->diag);
217:   jac->bs  = A->rmap->bs;
218:   jac->mbs = A->rmap->n/A->rmap->bs;
219:   switch (jac->bs) {
220:   case 1:
221:     pc->ops->apply = PCApply_PBJacobi_1;
222:     break;
223:   case 2:
224:     pc->ops->apply = PCApply_PBJacobi_2;
225:     break;
226:   case 3:
227:     pc->ops->apply = PCApply_PBJacobi_3;
228:     break;
229:   case 4:
230:     pc->ops->apply = PCApply_PBJacobi_4;
231:     break;
232:   case 5:
233:     pc->ops->apply = PCApply_PBJacobi_5;
234:     break;
235:   case 6:
236:     pc->ops->apply = PCApply_PBJacobi_6;
237:     break;
238:   case 7:
239:     pc->ops->apply = PCApply_PBJacobi_7;
240:     break;
241:   default:
242:     SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
243:   }
244:   return(0);
245: }
246: /* -------------------------------------------------------------------------- */
249: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
250: {

254:   /*
255:       Free the private data structure that was hanging off the PC
256:   */
257:   PetscFree(pc->data);
258:   return(0);
259: }

263: static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer)
264: {
266:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
267:   PetscBool      iascii;

270:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
271:   if (iascii) {
272:     PetscViewerASCIIPrintf(viewer,"  point-block Jacobi: block size %D\n",jac->bs);
273:   }
274:   return(0);
275: }

277: /* -------------------------------------------------------------------------- */
278: /*MC
279:      PCPBJACOBI - Point block Jacobi

281:    Level: beginner

283:   Concepts: point block Jacobi


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

288: M*/

292: PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc)
293: {
294:   PC_PBJacobi    *jac;

298:   /*
299:      Creates the private data structure for this preconditioner and
300:      attach it to the PC object.
301:   */
302:   PetscNewLog(pc,PC_PBJacobi,&jac);
303:   pc->data = (void*)jac;

305:   /*
306:      Initialize the pointers to vectors to ZERO; these will be used to store
307:      diagonal entries of the matrix for fast preconditioner application.
308:   */
309:   jac->diag = 0;

311:   /*
312:       Set the pointers for the functions that are provided above.
313:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
314:       are called, they will automatically call these functions.  Note we
315:       choose not to provide a couple of these functions since they are
316:       not needed.
317:   */
318:   pc->ops->apply               = 0; /*set depending on the block size */
319:   pc->ops->applytranspose      = 0;
320:   pc->ops->setup               = PCSetUp_PBJacobi;
321:   pc->ops->destroy             = PCDestroy_PBJacobi;
322:   pc->ops->setfromoptions      = 0;
323:   pc->ops->view                = PCView_PBJacobi;
324:   pc->ops->applyrichardson     = 0;
325:   pc->ops->applysymmetricleft  = 0;
326:   pc->ops->applysymmetricright = 0;
327:   return(0);
328: }