Actual source code: bjacobi.c


  2: /*
  3:    Defines a block Jacobi preconditioner.
  4: */

  6: #include <../src/ksp/pc/impls/bjacobi/bjacobi.h>

  8: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC,Mat,Mat);
  9: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC,Mat,Mat);
 10: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC);

 12: static PetscErrorCode PCSetUp_BJacobi(PC pc)
 13: {
 14:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
 15:   Mat            mat  = pc->mat,pmat = pc->pmat;
 16:   PetscBool      hasop;
 17:   PetscInt       N,M,start,i,sum,end;
 18:   PetscInt       bs,i_start=-1,i_end=-1;
 19:   PetscMPIInt    rank,size;

 21:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
 22:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
 23:   MatGetLocalSize(pc->pmat,&M,&N);
 24:   MatGetBlockSize(pc->pmat,&bs);

 26:   if (jac->n > 0 && jac->n < size) {
 27:     PCSetUp_BJacobi_Multiproc(pc);
 28:     return 0;
 29:   }

 31:   /* --------------------------------------------------------------------------
 32:       Determines the number of blocks assigned to each processor
 33:   -----------------------------------------------------------------------------*/

 35:   /*   local block count  given */
 36:   if (jac->n_local > 0 && jac->n < 0) {
 37:     MPIU_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
 38:     if (jac->l_lens) { /* check that user set these correctly */
 39:       sum = 0;
 40:       for (i=0; i<jac->n_local; i++) {
 42:         sum += jac->l_lens[i];
 43:       }
 45:     } else {
 46:       PetscMalloc1(jac->n_local,&jac->l_lens);
 47:       for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
 48:     }
 49:   } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
 50:     /* global blocks given: determine which ones are local */
 51:     if (jac->g_lens) {
 52:       /* check if the g_lens is has valid entries */
 53:       for (i=0; i<jac->n; i++) {
 56:       }
 57:       if (size == 1) {
 58:         jac->n_local = jac->n;
 59:         PetscMalloc1(jac->n_local,&jac->l_lens);
 60:         PetscArraycpy(jac->l_lens,jac->g_lens,jac->n_local);
 61:         /* check that user set these correctly */
 62:         sum = 0;
 63:         for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
 65:       } else {
 66:         MatGetOwnershipRange(pc->pmat,&start,&end);
 67:         /* loop over blocks determing first one owned by me */
 68:         sum = 0;
 69:         for (i=0; i<jac->n+1; i++) {
 70:           if (sum == start) { i_start = i; goto start_1;}
 71:           if (i < jac->n) sum += jac->g_lens[i];
 72:         }
 73:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
 74: start_1:
 75:         for (i=i_start; i<jac->n+1; i++) {
 76:           if (sum == end) { i_end = i; goto end_1; }
 77:           if (i < jac->n) sum += jac->g_lens[i];
 78:         }
 79:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
 80: end_1:
 81:         jac->n_local = i_end - i_start;
 82:         PetscMalloc1(jac->n_local,&jac->l_lens);
 83:         PetscArraycpy(jac->l_lens,jac->g_lens+i_start,jac->n_local);
 84:       }
 85:     } else { /* no global blocks given, determine then using default layout */
 86:       jac->n_local = jac->n/size + ((jac->n % size) > rank);
 87:       PetscMalloc1(jac->n_local,&jac->l_lens);
 88:       for (i=0; i<jac->n_local; i++) {
 89:         jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
 91:       }
 92:     }
 93:   } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
 94:     jac->n         = size;
 95:     jac->n_local   = 1;
 96:     PetscMalloc1(1,&jac->l_lens);
 97:     jac->l_lens[0] = M;
 98:   } else { /* jac->n > 0 && jac->n_local > 0 */
 99:     if (!jac->l_lens) {
100:       PetscMalloc1(jac->n_local,&jac->l_lens);
101:       for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
102:     }
103:   }

106:   /* -------------------------
107:       Determines mat and pmat
108:   ---------------------------*/
109:   MatHasOperation(pc->mat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
110:   if (!hasop && size == 1) {
111:     mat  = pc->mat;
112:     pmat = pc->pmat;
113:   } else {
114:     if (pc->useAmat) {
115:       /* use block from Amat matrix, not Pmat for local MatMult() */
116:       MatGetDiagonalBlock(pc->mat,&mat);
117:     }
118:     if (pc->pmat != pc->mat || !pc->useAmat) {
119:       MatGetDiagonalBlock(pc->pmat,&pmat);
120:     } else pmat = mat;
121:   }

123:   /* ------
124:      Setup code depends on the number of blocks
125:   */
126:   if (jac->n_local == 1) {
127:     PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
128:   } else {
129:     PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
130:   }
131:   return 0;
132: }

134: /* Default destroy, if it has never been setup */
135: static PetscErrorCode PCDestroy_BJacobi(PC pc)
136: {
137:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

139:   PetscFree(jac->g_lens);
140:   PetscFree(jac->l_lens);
141:   PetscFree(pc->data);
142:   return 0;
143: }

145: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
146: {
147:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
148:   PetscInt       blocks,i;
149:   PetscBool      flg;

151:   PetscOptionsHead(PetscOptionsObject,"Block Jacobi options");
152:   PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
153:   if (flg) PCBJacobiSetTotalBlocks(pc,blocks,NULL);
154:   PetscOptionsInt("-pc_bjacobi_local_blocks","Local number of blocks","PCBJacobiSetLocalBlocks",jac->n_local,&blocks,&flg);
155:   if (flg) PCBJacobiSetLocalBlocks(pc,blocks,NULL);
156:   if (jac->ksp) {
157:     /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
158:      * unless we had already been called. */
159:     for (i=0; i<jac->n_local; i++) {
160:       KSPSetFromOptions(jac->ksp[i]);
161:     }
162:   }
163:   PetscOptionsTail();
164:   return 0;
165: }

167: #include <petscdraw.h>
168: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
169: {
170:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
171:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
172:   PetscErrorCode       ierr;
173:   PetscMPIInt          rank;
174:   PetscInt             i;
175:   PetscBool            iascii,isstring,isdraw;
176:   PetscViewer          sviewer;
177:   PetscViewerFormat    format;
178:   const char           *prefix;

180:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
181:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
182:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
183:   if (iascii) {
184:     if (pc->useAmat) {
185:       PetscViewerASCIIPrintf(viewer,"  using Amat local matrix, number of blocks = %D\n",jac->n);
186:     }
187:     PetscViewerASCIIPrintf(viewer,"  number of blocks = %D\n",jac->n);
188:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
189:     PetscViewerGetFormat(viewer,&format);
190:     if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
191:       PetscViewerASCIIPrintf(viewer,"  Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
192:       PCGetOptionsPrefix(pc,&prefix);
193:       PetscViewerASCIIPrintf(viewer,"  Use -%sksp_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
194:       if (jac->ksp && !jac->psubcomm) {
195:         PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
196:         if (rank == 0) {
197:           PetscViewerASCIIPushTab(viewer);
198:           KSPView(jac->ksp[0],sviewer);
199:           PetscViewerASCIIPopTab(viewer);
200:         }
201:         PetscViewerFlush(sviewer);
202:         PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
203:         PetscViewerFlush(viewer);
204:         /*  extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
205:         PetscViewerASCIIPopSynchronized(viewer);
206:       } else if (mpjac && jac->ksp && mpjac->psubcomm) {
207:         PetscViewerGetSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
208:         if (!mpjac->psubcomm->color) {
209:           PetscViewerASCIIPushTab(viewer);
210:           KSPView(*(jac->ksp),sviewer);
211:           PetscViewerASCIIPopTab(viewer);
212:         }
213:         PetscViewerFlush(sviewer);
214:         PetscViewerRestoreSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
215:         PetscViewerFlush(viewer);
216:         /*  extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
217:         PetscViewerASCIIPopSynchronized(viewer);
218:       } else {
219:         PetscViewerFlush(viewer);
220:       }
221:     } else {
222:       PetscInt n_global;
223:       MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
224:       PetscViewerASCIIPushSynchronized(viewer);
225:       PetscViewerASCIIPrintf(viewer,"  Local solver information for each block is in the following KSP and PC objects:\n");
226:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
227:                                                 rank,jac->n_local,jac->first_local);
228:       PetscViewerASCIIPushTab(viewer);
229:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
230:       for (i=0; i<jac->n_local; i++) {
231:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
232:         KSPView(jac->ksp[i],sviewer);
233:         PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
234:       }
235:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
236:       PetscViewerASCIIPopTab(viewer);
237:       PetscViewerFlush(viewer);
238:       PetscViewerASCIIPopSynchronized(viewer);
239:     }
240:   } else if (isstring) {
241:     PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
242:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
243:     if (jac->ksp) KSPView(jac->ksp[0],sviewer);
244:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
245:   } else if (isdraw) {
246:     PetscDraw draw;
247:     char      str[25];
248:     PetscReal x,y,bottom,h;

250:     PetscViewerDrawGetDraw(viewer,0,&draw);
251:     PetscDrawGetCurrentPoint(draw,&x,&y);
252:     PetscSNPrintf(str,25,"Number blocks %D",jac->n);
253:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
254:     bottom = y - h;
255:     PetscDrawPushCurrentPoint(draw,x,bottom);
256:     /* warning the communicator on viewer is different then on ksp in parallel */
257:     if (jac->ksp) KSPView(jac->ksp[0],viewer);
258:     PetscDrawPopCurrentPoint(draw);
259:   }
260:   return 0;
261: }

263: /* -------------------------------------------------------------------------------------*/

265: static PetscErrorCode  PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
266: {
267:   PC_BJacobi *jac = (PC_BJacobi*)pc->data;


271:   if (n_local) *n_local = jac->n_local;
272:   if (first_local) *first_local = jac->first_local;
273:   if (ksp) *ksp                 = jac->ksp;
274:   return 0;
275: }

277: static PetscErrorCode  PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
278: {
279:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

282:   jac->n = blocks;
283:   if (!lens) jac->g_lens = NULL;
284:   else {
285:     PetscMalloc1(blocks,&jac->g_lens);
286:     PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
287:     PetscArraycpy(jac->g_lens,lens,blocks);
288:   }
289:   return 0;
290: }

292: static PetscErrorCode  PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
293: {
294:   PC_BJacobi *jac = (PC_BJacobi*) pc->data;

296:   *blocks = jac->n;
297:   if (lens) *lens = jac->g_lens;
298:   return 0;
299: }

301: static PetscErrorCode  PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
302: {
303:   PC_BJacobi     *jac;

305:   jac = (PC_BJacobi*)pc->data;

307:   jac->n_local = blocks;
308:   if (!lens) jac->l_lens = NULL;
309:   else {
310:     PetscMalloc1(blocks,&jac->l_lens);
311:     PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
312:     PetscArraycpy(jac->l_lens,lens,blocks);
313:   }
314:   return 0;
315: }

317: static PetscErrorCode  PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
318: {
319:   PC_BJacobi *jac = (PC_BJacobi*) pc->data;

321:   *blocks = jac->n_local;
322:   if (lens) *lens = jac->l_lens;
323:   return 0;
324: }

326: /* -------------------------------------------------------------------------------------*/

328: /*@C
329:    PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
330:    this processor.

332:    Not Collective

334:    Input Parameter:
335: .  pc - the preconditioner context

337:    Output Parameters:
338: +  n_local - the number of blocks on this processor, or NULL
339: .  first_local - the global number of the first block on this processor, or NULL
340: -  ksp - the array of KSP contexts

342:    Notes:
343:    After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.

345:    Currently for some matrix implementations only 1 block per processor
346:    is supported.

348:    You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().

350:    Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
351:       You can call PCBJacobiGetSubKSP(pc,nlocal,firstlocal,PETSC_NULL_KSP,ierr) to determine how large the
352:       KSP array must be.

354:    Level: advanced

356: .seealso: PCASMGetSubKSP()
357: @*/
358: PetscErrorCode  PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
359: {
361:   PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
362:   return 0;
363: }

365: /*@
366:    PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
367:    Jacobi preconditioner.

369:    Collective on PC

371:    Input Parameters:
372: +  pc - the preconditioner context
373: .  blocks - the number of blocks
374: -  lens - [optional] integer array containing the size of each block

376:    Options Database Key:
377: .  -pc_bjacobi_blocks <blocks> - Sets the number of global blocks

379:    Notes:
380:    Currently only a limited number of blocking configurations are supported.
381:    All processors sharing the PC must call this routine with the same data.

383:    Level: intermediate

385: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
386: @*/
387: PetscErrorCode  PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
388: {
391:   PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
392:   return 0;
393: }

395: /*@C
396:    PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
397:    Jacobi preconditioner.

399:    Not Collective

401:    Input Parameter:
402: .  pc - the preconditioner context

404:    Output parameters:
405: +  blocks - the number of blocks
406: -  lens - integer array containing the size of each block

408:    Level: intermediate

410: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
411: @*/
412: PetscErrorCode  PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
413: {
416:   PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
417:   return 0;
418: }

420: /*@
421:    PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
422:    Jacobi preconditioner.

424:    Not Collective

426:    Input Parameters:
427: +  pc - the preconditioner context
428: .  blocks - the number of blocks
429: -  lens - [optional] integer array containing size of each block

431:    Options Database Key:
432: .  -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks

434:    Note:
435:    Currently only a limited number of blocking configurations are supported.

437:    Level: intermediate

439: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
440: @*/
441: PetscErrorCode  PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
442: {
445:   PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
446:   return 0;
447: }

449: /*@C
450:    PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
451:    Jacobi preconditioner.

453:    Not Collective

455:    Input Parameters:
456: +  pc - the preconditioner context
457: .  blocks - the number of blocks
458: -  lens - [optional] integer array containing size of each block

460:    Note:
461:    Currently only a limited number of blocking configurations are supported.

463:    Level: intermediate

465: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
466: @*/
467: PetscErrorCode  PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
468: {
471:   PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
472:   return 0;
473: }

475: /* -----------------------------------------------------------------------------------*/

477: /*MC
478:    PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
479:            its own KSP object.

481:    Options Database Keys:
482: +  -pc_use_amat - use Amat to apply block of operator in inner Krylov method
483: -  -pc_bjacobi_blocks <n> - use n total blocks

485:    Notes:
486:      Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.

488:      To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
489:         options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly

491:      To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
492:          and set the options directly on the resulting KSP object (you can access its PC
493:          KSPGetPC())

495:      For GPU-based vectors (CUDA, ViennaCL) it is recommended to use exactly one block per MPI process for best
496:          performance.  Different block partitioning may lead to additional data transfers
497:          between host and GPU that lead to degraded performance.

499:      The options prefix for each block is sub_, for example -sub_pc_type lu.

501:      When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.

503:      See PCJACOBI for point Jacobi preconditioning, PCVPBJACOBI for variable size point block Jacobi and PCPBJACOBI for large blocks

505:    Level: beginner

507: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
508:            PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
509:            PCBJacobiSetLocalBlocks(), PCSetModifySubMatrices(), PCJACOBI, PCVPBJACOBI, PCPBJACOBI
510: M*/

512: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
513: {
514:   PetscMPIInt    rank;
515:   PC_BJacobi     *jac;

517:   PetscNewLog(pc,&jac);
518:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);

520:   pc->ops->apply           = NULL;
521:   pc->ops->matapply        = NULL;
522:   pc->ops->applytranspose  = NULL;
523:   pc->ops->setup           = PCSetUp_BJacobi;
524:   pc->ops->destroy         = PCDestroy_BJacobi;
525:   pc->ops->setfromoptions  = PCSetFromOptions_BJacobi;
526:   pc->ops->view            = PCView_BJacobi;
527:   pc->ops->applyrichardson = NULL;

529:   pc->data               = (void*)jac;
530:   jac->n                 = -1;
531:   jac->n_local           = -1;
532:   jac->first_local       = rank;
533:   jac->ksp               = NULL;
534:   jac->g_lens            = NULL;
535:   jac->l_lens            = NULL;
536:   jac->psubcomm          = NULL;

538:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
539:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
540:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
541:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
542:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
543:   return 0;
544: }

546: /* --------------------------------------------------------------------------------------------*/
547: /*
548:         These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
549: */
550: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
551: {
552:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
553:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;

555:   KSPReset(jac->ksp[0]);
556:   VecDestroy(&bjac->x);
557:   VecDestroy(&bjac->y);
558:   return 0;
559: }

561: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
562: {
563:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
564:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;

566:   PCReset_BJacobi_Singleblock(pc);
567:   KSPDestroy(&jac->ksp[0]);
568:   PetscFree(jac->ksp);
569:   PetscFree(jac->l_lens);
570:   PetscFree(jac->g_lens);
571:   PetscFree(bjac);
572:   PetscFree(pc->data);
573:   return 0;
574: }

576: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
577: {
578:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
579:   KSP                subksp = jac->ksp[0];
580:   KSPConvergedReason reason;

582:   KSPSetUp(subksp);
583:   KSPGetConvergedReason(subksp,&reason);
584:   if (reason == KSP_DIVERGED_PC_FAILED) {
585:     pc->failedreason = PC_SUBPC_ERROR;
586:   }
587:   return 0;
588: }

590: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
591: {
592:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
593:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;

595:   VecGetLocalVectorRead(x, bjac->x);
596:   VecGetLocalVector(y, bjac->y);
597:   /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
598:      matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
599:      of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
600:   KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
601:   KSPSolve(jac->ksp[0],bjac->x,bjac->y);
602:   KSPCheckSolve(jac->ksp[0],pc,bjac->y);
603:   VecRestoreLocalVectorRead(x, bjac->x);
604:   VecRestoreLocalVector(y, bjac->y);
605:   return 0;
606: }

608: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc,Mat X,Mat Y)
609: {
610:   PC_BJacobi     *jac  = (PC_BJacobi*)pc->data;
611:   Mat            sX,sY;

613:   /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
614:      matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
615:      of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
616:   KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
617:   MatDenseGetLocalMatrix(X,&sX);
618:   MatDenseGetLocalMatrix(Y,&sY);
619:   KSPMatSolve(jac->ksp[0],sX,sY);
620:   return 0;
621: }

623: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
624: {
625:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
626:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
627:   PetscScalar            *y_array;
628:   const PetscScalar      *x_array;
629:   PC                     subpc;

631:   /*
632:       The VecPlaceArray() is to avoid having to copy the
633:     y vector into the bjac->x vector. The reason for
634:     the bjac->x vector is that we need a sequential vector
635:     for the sequential solve.
636:   */
637:   VecGetArrayRead(x,&x_array);
638:   VecGetArray(y,&y_array);
639:   VecPlaceArray(bjac->x,x_array);
640:   VecPlaceArray(bjac->y,y_array);
641:   /* apply the symmetric left portion of the inner PC operator */
642:   /* note this by-passes the inner KSP and its options completely */
643:   KSPGetPC(jac->ksp[0],&subpc);
644:   PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
645:   VecResetArray(bjac->x);
646:   VecResetArray(bjac->y);
647:   VecRestoreArrayRead(x,&x_array);
648:   VecRestoreArray(y,&y_array);
649:   return 0;
650: }

652: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
653: {
654:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
655:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
656:   PetscScalar            *y_array;
657:   const PetscScalar      *x_array;
658:   PC                     subpc;

660:   /*
661:       The VecPlaceArray() is to avoid having to copy the
662:     y vector into the bjac->x vector. The reason for
663:     the bjac->x vector is that we need a sequential vector
664:     for the sequential solve.
665:   */
666:   VecGetArrayRead(x,&x_array);
667:   VecGetArray(y,&y_array);
668:   VecPlaceArray(bjac->x,x_array);
669:   VecPlaceArray(bjac->y,y_array);

671:   /* apply the symmetric right portion of the inner PC operator */
672:   /* note this by-passes the inner KSP and its options completely */

674:   KSPGetPC(jac->ksp[0],&subpc);
675:   PCApplySymmetricRight(subpc,bjac->x,bjac->y);

677:   VecRestoreArrayRead(x,&x_array);
678:   VecRestoreArray(y,&y_array);
679:   return 0;
680: }

682: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
683: {
684:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
685:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
686:   PetscScalar            *y_array;
687:   const PetscScalar      *x_array;

689:   /*
690:       The VecPlaceArray() is to avoid having to copy the
691:     y vector into the bjac->x vector. The reason for
692:     the bjac->x vector is that we need a sequential vector
693:     for the sequential solve.
694:   */
695:   VecGetArrayRead(x,&x_array);
696:   VecGetArray(y,&y_array);
697:   VecPlaceArray(bjac->x,x_array);
698:   VecPlaceArray(bjac->y,y_array);
699:   KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
700:   KSPCheckSolve(jac->ksp[0],pc,bjac->y);
701:   VecResetArray(bjac->x);
702:   VecResetArray(bjac->y);
703:   VecRestoreArrayRead(x,&x_array);
704:   VecRestoreArray(y,&y_array);
705:   return 0;
706: }

708: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
709: {
710:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
711:   PetscInt               m;
712:   KSP                    ksp;
713:   PC_BJacobi_Singleblock *bjac;
714:   PetscBool              wasSetup = PETSC_TRUE;
715:   VecType                vectype;
716:   const char             *prefix;

718:   if (!pc->setupcalled) {
719:     if (!jac->ksp) {
720:       wasSetup = PETSC_FALSE;

722:       KSPCreate(PETSC_COMM_SELF,&ksp);
723:       KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
724:       PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
725:       PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
726:       KSPSetType(ksp,KSPPREONLY);
727:       PCGetOptionsPrefix(pc,&prefix);
728:       KSPSetOptionsPrefix(ksp,prefix);
729:       KSPAppendOptionsPrefix(ksp,"sub_");

731:       pc->ops->reset               = PCReset_BJacobi_Singleblock;
732:       pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
733:       pc->ops->apply               = PCApply_BJacobi_Singleblock;
734:       pc->ops->matapply            = PCMatApply_BJacobi_Singleblock;
735:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
736:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
737:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
738:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

740:       PetscMalloc1(1,&jac->ksp);
741:       jac->ksp[0] = ksp;

743:       PetscNewLog(pc,&bjac);
744:       jac->data = (void*)bjac;
745:     } else {
746:       ksp  = jac->ksp[0];
747:       bjac = (PC_BJacobi_Singleblock*)jac->data;
748:     }

750:     /*
751:       The reason we need to generate these vectors is to serve
752:       as the right-hand side and solution vector for the solve on the
753:       block. We do not need to allocate space for the vectors since
754:       that is provided via VecPlaceArray() just before the call to
755:       KSPSolve() on the block.
756:     */
757:     MatGetSize(pmat,&m,&m);
758:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
759:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
760:     MatGetVecType(pmat,&vectype);
761:     VecSetType(bjac->x,vectype);
762:     VecSetType(bjac->y,vectype);
763:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
764:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
765:   } else {
766:     ksp  = jac->ksp[0];
767:     bjac = (PC_BJacobi_Singleblock*)jac->data;
768:   }
769:   KSPGetOptionsPrefix(ksp,&prefix);
770:   if (pc->useAmat) {
771:     KSPSetOperators(ksp,mat,pmat);
772:     MatSetOptionsPrefix(mat,prefix);
773:   } else {
774:     KSPSetOperators(ksp,pmat,pmat);
775:   }
776:   MatSetOptionsPrefix(pmat,prefix);
777:   if (!wasSetup && pc->setfromoptionscalled) {
778:     /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
779:     KSPSetFromOptions(ksp);
780:   }
781:   return 0;
782: }

784: /* ---------------------------------------------------------------------------------------------*/
785: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
786: {
787:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
788:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
789:   PetscInt              i;

791:   if (bjac && bjac->pmat) {
792:     MatDestroyMatrices(jac->n_local,&bjac->pmat);
793:     if (pc->useAmat) {
794:       MatDestroyMatrices(jac->n_local,&bjac->mat);
795:     }
796:   }

798:   for (i=0; i<jac->n_local; i++) {
799:     KSPReset(jac->ksp[i]);
800:     if (bjac && bjac->x) {
801:       VecDestroy(&bjac->x[i]);
802:       VecDestroy(&bjac->y[i]);
803:       ISDestroy(&bjac->is[i]);
804:     }
805:   }
806:   PetscFree(jac->l_lens);
807:   PetscFree(jac->g_lens);
808:   return 0;
809: }

811: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
812: {
813:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
814:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
815:   PetscInt              i;

817:   PCReset_BJacobi_Multiblock(pc);
818:   if (bjac) {
819:     PetscFree2(bjac->x,bjac->y);
820:     PetscFree(bjac->starts);
821:     PetscFree(bjac->is);
822:   }
823:   PetscFree(jac->data);
824:   for (i=0; i<jac->n_local; i++) {
825:     KSPDestroy(&jac->ksp[i]);
826:   }
827:   PetscFree(jac->ksp);
828:   PetscFree(pc->data);
829:   return 0;
830: }

832: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
833: {
834:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
835:   PetscInt           i,n_local = jac->n_local;
836:   KSPConvergedReason reason;

838:   for (i=0; i<n_local; i++) {
839:     KSPSetUp(jac->ksp[i]);
840:     KSPGetConvergedReason(jac->ksp[i],&reason);
841:     if (reason == KSP_DIVERGED_PC_FAILED) {
842:       pc->failedreason = PC_SUBPC_ERROR;
843:     }
844:   }
845:   return 0;
846: }

848: /*
849:       Preconditioner for block Jacobi
850: */
851: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
852: {
853:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
854:   PetscInt              i,n_local = jac->n_local;
855:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
856:   PetscScalar           *yin;
857:   const PetscScalar     *xin;

859:   VecGetArrayRead(x,&xin);
860:   VecGetArray(y,&yin);
861:   for (i=0; i<n_local; i++) {
862:     /*
863:        To avoid copying the subvector from x into a workspace we instead
864:        make the workspace vector array point to the subpart of the array of
865:        the global vector.
866:     */
867:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
868:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

870:     PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
871:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
872:     KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
873:     PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

875:     VecResetArray(bjac->x[i]);
876:     VecResetArray(bjac->y[i]);
877:   }
878:   VecRestoreArrayRead(x,&xin);
879:   VecRestoreArray(y,&yin);
880:   return 0;
881: }

883: /*
884:       Preconditioner for block Jacobi
885: */
886: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
887: {
888:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
889:   PetscInt              i,n_local = jac->n_local;
890:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
891:   PetscScalar           *yin;
892:   const PetscScalar     *xin;

894:   VecGetArrayRead(x,&xin);
895:   VecGetArray(y,&yin);
896:   for (i=0; i<n_local; i++) {
897:     /*
898:        To avoid copying the subvector from x into a workspace we instead
899:        make the workspace vector array point to the subpart of the array of
900:        the global vector.
901:     */
902:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
903:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

905:     PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
906:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
907:     KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
908:     PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

910:     VecResetArray(bjac->x[i]);
911:     VecResetArray(bjac->y[i]);
912:   }
913:   VecRestoreArrayRead(x,&xin);
914:   VecRestoreArray(y,&yin);
915:   return 0;
916: }

918: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
919: {
920:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
921:   PetscInt              m,n_local,N,M,start,i;
922:   const char            *prefix;
923:   KSP                   ksp;
924:   Vec                   x,y;
925:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
926:   PC                    subpc;
927:   IS                    is;
928:   MatReuse              scall;
929:   VecType               vectype;

931:   MatGetLocalSize(pc->pmat,&M,&N);

933:   n_local = jac->n_local;

935:   if (pc->useAmat) {
936:     PetscBool same;
937:     PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
939:   }

941:   if (!pc->setupcalled) {
942:     scall = MAT_INITIAL_MATRIX;

944:     if (!jac->ksp) {
945:       pc->ops->reset         = PCReset_BJacobi_Multiblock;
946:       pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
947:       pc->ops->apply         = PCApply_BJacobi_Multiblock;
948:       pc->ops->matapply      = NULL;
949:       pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
950:       pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

952:       PetscNewLog(pc,&bjac);
953:       PetscMalloc1(n_local,&jac->ksp);
954:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
955:       PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
956:       PetscMalloc1(n_local,&bjac->starts);
957:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));

959:       jac->data = (void*)bjac;
960:       PetscMalloc1(n_local,&bjac->is);
961:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));

963:       for (i=0; i<n_local; i++) {
964:         KSPCreate(PETSC_COMM_SELF,&ksp);
965:         KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
966:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
967:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
968:         KSPSetType(ksp,KSPPREONLY);
969:         KSPGetPC(ksp,&subpc);
970:         PCGetOptionsPrefix(pc,&prefix);
971:         KSPSetOptionsPrefix(ksp,prefix);
972:         KSPAppendOptionsPrefix(ksp,"sub_");

974:         jac->ksp[i] = ksp;
975:       }
976:     } else {
977:       bjac = (PC_BJacobi_Multiblock*)jac->data;
978:     }

980:     start = 0;
981:     MatGetVecType(pmat,&vectype);
982:     for (i=0; i<n_local; i++) {
983:       m = jac->l_lens[i];
984:       /*
985:       The reason we need to generate these vectors is to serve
986:       as the right-hand side and solution vector for the solve on the
987:       block. We do not need to allocate space for the vectors since
988:       that is provided via VecPlaceArray() just before the call to
989:       KSPSolve() on the block.

991:       */
992:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
993:       VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
994:       VecSetType(x,vectype);
995:       VecSetType(y,vectype);
996:       PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
997:       PetscLogObjectParent((PetscObject)pc,(PetscObject)y);

999:       bjac->x[i]      = x;
1000:       bjac->y[i]      = y;
1001:       bjac->starts[i] = start;

1003:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1004:       bjac->is[i] = is;
1005:       PetscLogObjectParent((PetscObject)pc,(PetscObject)is);

1007:       start += m;
1008:     }
1009:   } else {
1010:     bjac = (PC_BJacobi_Multiblock*)jac->data;
1011:     /*
1012:        Destroy the blocks from the previous iteration
1013:     */
1014:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1015:       MatDestroyMatrices(n_local,&bjac->pmat);
1016:       if (pc->useAmat) {
1017:         MatDestroyMatrices(n_local,&bjac->mat);
1018:       }
1019:       scall = MAT_INITIAL_MATRIX;
1020:     } else scall = MAT_REUSE_MATRIX;
1021:   }

1023:   MatCreateSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1024:   if (pc->useAmat) {
1025:     MatCreateSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1026:   }
1027:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
1028:      different boundary conditions for the submatrices than for the global problem) */
1029:   PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);

1031:   for (i=0; i<n_local; i++) {
1032:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1033:     KSPGetOptionsPrefix(jac->ksp[i],&prefix);
1034:     if (pc->useAmat) {
1035:       PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1036:       KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1037:       MatSetOptionsPrefix(bjac->mat[i],prefix);
1038:     } else {
1039:       KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1040:     }
1041:     MatSetOptionsPrefix(bjac->pmat[i],prefix);
1042:     if (pc->setfromoptionscalled) {
1043:       KSPSetFromOptions(jac->ksp[i]);
1044:     }
1045:   }
1046:   return 0;
1047: }

1049: /* ---------------------------------------------------------------------------------------------*/
1050: /*
1051:       These are for a single block with multiple processes
1052: */
1053: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiproc(PC pc)
1054: {
1055:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
1056:   KSP                subksp = jac->ksp[0];
1057:   KSPConvergedReason reason;

1059:   KSPSetUp(subksp);
1060:   KSPGetConvergedReason(subksp,&reason);
1061:   if (reason == KSP_DIVERGED_PC_FAILED) {
1062:     pc->failedreason = PC_SUBPC_ERROR;
1063:   }
1064:   return 0;
1065: }

1067: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1068: {
1069:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1070:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;

1072:   VecDestroy(&mpjac->ysub);
1073:   VecDestroy(&mpjac->xsub);
1074:   MatDestroy(&mpjac->submats);
1075:   if (jac->ksp) KSPReset(jac->ksp[0]);
1076:   return 0;
1077: }

1079: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1080: {
1081:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1082:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;

1084:   PCReset_BJacobi_Multiproc(pc);
1085:   KSPDestroy(&jac->ksp[0]);
1086:   PetscFree(jac->ksp);
1087:   PetscSubcommDestroy(&mpjac->psubcomm);

1089:   PetscFree(mpjac);
1090:   PetscFree(pc->data);
1091:   return 0;
1092: }

1094: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1095: {
1096:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1097:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1098:   PetscScalar          *yarray;
1099:   const PetscScalar    *xarray;
1100:   KSPConvergedReason   reason;

1102:   /* place x's and y's local arrays into xsub and ysub */
1103:   VecGetArrayRead(x,&xarray);
1104:   VecGetArray(y,&yarray);
1105:   VecPlaceArray(mpjac->xsub,xarray);
1106:   VecPlaceArray(mpjac->ysub,yarray);

1108:   /* apply preconditioner on each matrix block */
1109:   PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1110:   KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1111:   KSPCheckSolve(jac->ksp[0],pc,mpjac->ysub);
1112:   PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1113:   KSPGetConvergedReason(jac->ksp[0],&reason);
1114:   if (reason == KSP_DIVERGED_PC_FAILED) {
1115:     pc->failedreason = PC_SUBPC_ERROR;
1116:   }

1118:   VecResetArray(mpjac->xsub);
1119:   VecResetArray(mpjac->ysub);
1120:   VecRestoreArrayRead(x,&xarray);
1121:   VecRestoreArray(y,&yarray);
1122:   return 0;
1123: }

1125: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc,Mat X,Mat Y)
1126: {
1127:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1128:   KSPConvergedReason   reason;
1129:   Mat                  sX,sY;
1130:   const PetscScalar    *x;
1131:   PetscScalar          *y;
1132:   PetscInt             m,N,lda,ldb;

1134:   /* apply preconditioner on each matrix block */
1135:   MatGetLocalSize(X,&m,NULL);
1136:   MatGetSize(X,NULL,&N);
1137:   MatDenseGetLDA(X,&lda);
1138:   MatDenseGetLDA(Y,&ldb);
1139:   MatDenseGetArrayRead(X,&x);
1140:   MatDenseGetArrayWrite(Y,&y);
1141:   MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,(PetscScalar*)x,&sX);
1142:   MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,y,&sY);
1143:   MatDenseSetLDA(sX,lda);
1144:   MatDenseSetLDA(sY,ldb);
1145:   PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1146:   KSPMatSolve(jac->ksp[0],sX,sY);
1147:   KSPCheckSolve(jac->ksp[0],pc,NULL);
1148:   PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1149:   MatDestroy(&sY);
1150:   MatDestroy(&sX);
1151:   MatDenseRestoreArrayWrite(Y,&y);
1152:   MatDenseRestoreArrayRead(X,&x);
1153:   KSPGetConvergedReason(jac->ksp[0],&reason);
1154:   if (reason == KSP_DIVERGED_PC_FAILED) {
1155:     pc->failedreason = PC_SUBPC_ERROR;
1156:   }
1157:   return 0;
1158: }

1160: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1161: {
1162:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1163:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1164:   PetscInt             m,n;
1165:   MPI_Comm             comm,subcomm=0;
1166:   const char           *prefix;
1167:   PetscBool            wasSetup = PETSC_TRUE;
1168:   VecType              vectype;

1170:   PetscObjectGetComm((PetscObject)pc,&comm);
1172:   jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1173:   if (!pc->setupcalled) {
1174:     wasSetup  = PETSC_FALSE;
1175:     PetscNewLog(pc,&mpjac);
1176:     jac->data = (void*)mpjac;

1178:     /* initialize datastructure mpjac */
1179:     if (!jac->psubcomm) {
1180:       /* Create default contiguous subcommunicatiors if user does not provide them */
1181:       PetscSubcommCreate(comm,&jac->psubcomm);
1182:       PetscSubcommSetNumber(jac->psubcomm,jac->n);
1183:       PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1184:       PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1185:     }
1186:     mpjac->psubcomm = jac->psubcomm;
1187:     subcomm         = PetscSubcommChild(mpjac->psubcomm);

1189:     /* Get matrix blocks of pmat */
1190:     MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);

1192:     /* create a new PC that processors in each subcomm have copy of */
1193:     PetscMalloc1(1,&jac->ksp);
1194:     KSPCreate(subcomm,&jac->ksp[0]);
1195:     KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1196:     PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1197:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1198:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1199:     KSPGetPC(jac->ksp[0],&mpjac->pc);

1201:     PCGetOptionsPrefix(pc,&prefix);
1202:     KSPSetOptionsPrefix(jac->ksp[0],prefix);
1203:     KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1204:     KSPGetOptionsPrefix(jac->ksp[0],&prefix);
1205:     MatSetOptionsPrefix(mpjac->submats,prefix);

1207:     /* create dummy vectors xsub and ysub */
1208:     MatGetLocalSize(mpjac->submats,&m,&n);
1209:     VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1210:     VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1211:     MatGetVecType(mpjac->submats,&vectype);
1212:     VecSetType(mpjac->xsub,vectype);
1213:     VecSetType(mpjac->ysub,vectype);
1214:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1215:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);

1217:     pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1218:     pc->ops->reset         = PCReset_BJacobi_Multiproc;
1219:     pc->ops->destroy       = PCDestroy_BJacobi_Multiproc;
1220:     pc->ops->apply         = PCApply_BJacobi_Multiproc;
1221:     pc->ops->matapply      = PCMatApply_BJacobi_Multiproc;
1222:   } else { /* pc->setupcalled */
1223:     subcomm = PetscSubcommChild(mpjac->psubcomm);
1224:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1225:       /* destroy old matrix blocks, then get new matrix blocks */
1226:       if (mpjac->submats) MatDestroy(&mpjac->submats);
1227:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1228:     } else {
1229:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1230:     }
1231:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1232:   }

1234:   if (!wasSetup && pc->setfromoptionscalled) {
1235:     KSPSetFromOptions(jac->ksp[0]);
1236:   }
1237:   return 0;
1238: }