Actual source code: bjacobi.c

petsc-master 2016-12-08
Report Typos and Errors

  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);

 14: static PetscErrorCode PCSetUp_BJacobi(PC pc)
 15: {
 16:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
 17:   Mat            mat  = pc->mat,pmat = pc->pmat;
 19:   void           (*f)(void);
 20:   PetscInt       N,M,start,i,sum,end;
 21:   PetscInt       bs,i_start=-1,i_end=-1;
 22:   PetscMPIInt    rank,size;
 23:   const char     *pprefix,*mprefix;

 26:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
 27:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
 28:   MatGetLocalSize(pc->pmat,&M,&N);
 29:   MatGetBlockSize(pc->pmat,&bs);

 31:   if (jac->n > 0 && jac->n < size) {
 32:     PCSetUp_BJacobi_Multiproc(pc);
 33:     return(0);
 34:   }

 36:   /* --------------------------------------------------------------------------
 37:       Determines the number of blocks assigned to each processor
 38:   -----------------------------------------------------------------------------*/

 40:   /*   local block count  given */
 41:   if (jac->n_local > 0 && jac->n < 0) {
 42:     MPIU_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
 43:     if (jac->l_lens) { /* check that user set these correctly */
 44:       sum = 0;
 45:       for (i=0; i<jac->n_local; i++) {
 46:         if (jac->l_lens[i]/bs*bs !=jac->l_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat blocksize doesn't match block Jacobi layout");
 47:         sum += jac->l_lens[i];
 48:       }
 49:       if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local lens set incorrectly");
 50:     } else {
 51:       PetscMalloc1(jac->n_local,&jac->l_lens);
 52:       for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
 53:     }
 54:   } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
 55:     /* global blocks given: determine which ones are local */
 56:     if (jac->g_lens) {
 57:       /* check if the g_lens is has valid entries */
 58:       for (i=0; i<jac->n; i++) {
 59:         if (!jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Zero block not allowed");
 60:         if (jac->g_lens[i]/bs*bs != jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat blocksize doesn't match block Jacobi layout");
 61:       }
 62:       if (size == 1) {
 63:         jac->n_local = jac->n;
 64:         PetscMalloc1(jac->n_local,&jac->l_lens);
 65:         PetscMemcpy(jac->l_lens,jac->g_lens,jac->n_local*sizeof(PetscInt));
 66:         /* check that user set these correctly */
 67:         sum = 0;
 68:         for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
 69:         if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Global lens set incorrectly");
 70:       } else {
 71:         MatGetOwnershipRange(pc->pmat,&start,&end);
 72:         /* loop over blocks determing first one owned by me */
 73:         sum = 0;
 74:         for (i=0; i<jac->n+1; i++) {
 75:           if (sum == start) { i_start = i; goto start_1;}
 76:           if (i < jac->n) sum += jac->g_lens[i];
 77:         }
 78:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
 79: start_1:
 80:         for (i=i_start; i<jac->n+1; i++) {
 81:           if (sum == end) { i_end = i; goto end_1; }
 82:           if (i < jac->n) sum += jac->g_lens[i];
 83:         }
 84:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
 85: end_1:
 86:         jac->n_local = i_end - i_start;
 87:         PetscMalloc1(jac->n_local,&jac->l_lens);
 88:         PetscMemcpy(jac->l_lens,jac->g_lens+i_start,jac->n_local*sizeof(PetscInt));
 89:       }
 90:     } else { /* no global blocks given, determine then using default layout */
 91:       jac->n_local = jac->n/size + ((jac->n % size) > rank);
 92:       PetscMalloc1(jac->n_local,&jac->l_lens);
 93:       for (i=0; i<jac->n_local; i++) {
 94:         jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
 95:         if (!jac->l_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Too many blocks given");
 96:       }
 97:     }
 98:   } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
 99:     jac->n         = size;
100:     jac->n_local   = 1;
101:     PetscMalloc1(1,&jac->l_lens);
102:     jac->l_lens[0] = M;
103:   } else { /* jac->n > 0 && jac->n_local > 0 */
104:     if (!jac->l_lens) {
105:       PetscMalloc1(jac->n_local,&jac->l_lens);
106:       for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
107:     }
108:   }
109:   if (jac->n_local < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of blocks is less than number of processors");

111:   /* -------------------------
112:       Determines mat and pmat
113:   ---------------------------*/
114:   MatShellGetOperation(pc->mat,MATOP_GET_DIAGONAL_BLOCK,&f);
115:   if (!f && size == 1) {
116:     mat  = pc->mat;
117:     pmat = pc->pmat;
118:   } else {
119:     if (pc->useAmat) {
120:       /* use block from Amat matrix, not Pmat for local MatMult() */
121:       MatGetDiagonalBlock(pc->mat,&mat);
122:       /* make submatrix have same prefix as entire matrix */
123:       PetscObjectGetOptionsPrefix((PetscObject)pc->mat,&mprefix);
124:       PetscObjectSetOptionsPrefix((PetscObject)mat,mprefix);
125:     }
126:     if (pc->pmat != pc->mat || !pc->useAmat) {
127:       MatGetDiagonalBlock(pc->pmat,&pmat);
128:       /* make submatrix have same prefix as entire matrix */
129:       PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
130:       PetscObjectSetOptionsPrefix((PetscObject)pmat,pprefix);
131:     } else pmat = mat;
132:   }

134:   /* ------
135:      Setup code depends on the number of blocks
136:   */
137:   if (jac->n_local == 1) {
138:     PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
139:   } else {
140:     PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
141:   }
142:   return(0);
143: }

145: /* Default destroy, if it has never been setup */
148: static PetscErrorCode PCDestroy_BJacobi(PC pc)
149: {
150:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

154:   PetscFree(jac->g_lens);
155:   PetscFree(jac->l_lens);
156:   PetscFree(pc->data);
157:   return(0);
158: }


163: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
164: {
165:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
167:   PetscInt       blocks,i;
168:   PetscBool      flg;

171:   PetscOptionsHead(PetscOptionsObject,"Block Jacobi options");
172:   PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
173:   if (flg) {
174:     PCBJacobiSetTotalBlocks(pc,blocks,NULL);
175:   }
176:   if (jac->ksp) {
177:     /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
178:      * unless we had already been called. */
179:     for (i=0; i<jac->n_local; i++) {
180:       KSPSetFromOptions(jac->ksp[i]);
181:     }
182:   }
183:   PetscOptionsTail();
184:   return(0);
185: }

187:  #include <petscdraw.h>
190: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
191: {
192:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
193:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
194:   PetscErrorCode       ierr;
195:   PetscMPIInt          rank;
196:   PetscInt             i;
197:   PetscBool            iascii,isstring,isdraw;
198:   PetscViewer          sviewer;

201:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
202:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
203:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
204:   if (iascii) {
205:     if (pc->useAmat) {
206:       PetscViewerASCIIPrintf(viewer,"  block Jacobi: using Amat local matrix, number of blocks = %D\n",jac->n);
207:     }
208:     PetscViewerASCIIPrintf(viewer,"  block Jacobi: number of blocks = %D\n",jac->n);
209:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
210:     if (jac->same_local_solves) {
211:       PetscViewerASCIIPrintf(viewer,"  Local solve is same for all blocks, in the following KSP and PC objects:\n");
212:       if (jac->ksp && !jac->psubcomm) {
213:         PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
214:         if (!rank) {
215:           PetscViewerASCIIPushTab(viewer);
216:           KSPView(jac->ksp[0],sviewer);
217:           PetscViewerASCIIPopTab(viewer);
218:         }
219:         PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
220:       } else if (mpjac && jac->ksp && mpjac->psubcomm) {
221:         PetscViewerGetSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
222:         if (!mpjac->psubcomm->color) {
223:           PetscViewerASCIIPushTab(viewer);
224:           KSPView(*(jac->ksp),sviewer);
225:           PetscViewerASCIIPopTab(viewer);
226:         }
227:         PetscViewerRestoreSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
228:       }
229:     } else {
230:       PetscInt n_global;
231:       MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
232:       PetscViewerASCIIPushSynchronized(viewer);
233:       PetscViewerASCIIPrintf(viewer,"  Local solve info for each block is in the following KSP and PC objects:\n");
234:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
235:                                                 rank,jac->n_local,jac->first_local);
236:       PetscViewerASCIIPushTab(viewer);
237:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
238:       for (i=0; i<jac->n_local; i++) {
239:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
240:         KSPView(jac->ksp[i],sviewer);
241:         PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
242:       }
243:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
244:       PetscViewerASCIIPopTab(viewer);
245:       PetscViewerFlush(viewer);
246:       PetscViewerASCIIPopSynchronized(viewer);
247:     }
248:   } else if (isstring) {
249:     PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
250:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
251:     if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
252:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
253:   } else if (isdraw) {
254:     PetscDraw draw;
255:     char      str[25];
256:     PetscReal x,y,bottom,h;

258:     PetscViewerDrawGetDraw(viewer,0,&draw);
259:     PetscDrawGetCurrentPoint(draw,&x,&y);
260:     PetscSNPrintf(str,25,"Number blocks %D",jac->n);
261:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
262:     bottom = y - h;
263:     PetscDrawPushCurrentPoint(draw,x,bottom);
264:     /* warning the communicator on viewer is different then on ksp in parallel */
265:     if (jac->ksp) {KSPView(jac->ksp[0],viewer);}
266:     PetscDrawPopCurrentPoint(draw);
267:   }
268:   return(0);
269: }

271: /* -------------------------------------------------------------------------------------*/

275: static PetscErrorCode  PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
276: {
277:   PC_BJacobi *jac = (PC_BJacobi*)pc->data;;

280:   if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() first");

282:   if (n_local) *n_local = jac->n_local;
283:   if (first_local) *first_local = jac->first_local;
284:   *ksp                   = jac->ksp;
285:   jac->same_local_solves = PETSC_FALSE;        /* Assume that local solves are now different;
286:                                                   not necessarily true though!  This flag is
287:                                                   used only for PCView_BJacobi() */
288:   return(0);
289: }

293: static PetscErrorCode  PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
294: {
295:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

299:   if (pc->setupcalled > 0 && jac->n!=blocks) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
300:   jac->n = blocks;
301:   if (!lens) jac->g_lens = 0;
302:   else {
303:     PetscMalloc1(blocks,&jac->g_lens);
304:     PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
305:     PetscMemcpy(jac->g_lens,lens,blocks*sizeof(PetscInt));
306:   }
307:   return(0);
308: }

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

317:   *blocks = jac->n;
318:   if (lens) *lens = jac->g_lens;
319:   return(0);
320: }

324: static PetscErrorCode  PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
325: {
326:   PC_BJacobi     *jac;

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

332:   jac->n_local = blocks;
333:   if (!lens) jac->l_lens = 0;
334:   else {
335:     PetscMalloc1(blocks,&jac->l_lens);
336:     PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
337:     PetscMemcpy(jac->l_lens,lens,blocks*sizeof(PetscInt));
338:   }
339:   return(0);
340: }

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

349:   *blocks = jac->n_local;
350:   if (lens) *lens = jac->l_lens;
351:   return(0);
352: }

354: /* -------------------------------------------------------------------------------------*/

358: /*@C
359:    PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
360:    this processor.

362:    Note Collective

364:    Input Parameter:
365: .  pc - the preconditioner context

367:    Output Parameters:
368: +  n_local - the number of blocks on this processor, or NULL
369: .  first_local - the global number of the first block on this processor, or NULL
370: -  ksp - the array of KSP contexts

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

375:    Currently for some matrix implementations only 1 block per processor
376:    is supported.

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

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

384:    Level: advanced

386: .keywords:  block, Jacobi, get, sub, KSP, context

388: .seealso: PCBJacobiGetSubKSP()
389: @*/
390: PetscErrorCode  PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
391: {

396:   PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
397:   return(0);
398: }

402: /*@
403:    PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
404:    Jacobi preconditioner.

406:    Collective on PC

408:    Input Parameters:
409: +  pc - the preconditioner context
410: .  blocks - the number of blocks
411: -  lens - [optional] integer array containing the size of each block

413:    Options Database Key:
414: .  -pc_bjacobi_blocks <blocks> - Sets the number of global blocks

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

420:    Level: intermediate

422: .keywords:  set, number, Jacobi, global, total, blocks

424: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
425: @*/
426: PetscErrorCode  PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
427: {

432:   if (blocks <= 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
433:   PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
434:   return(0);
435: }

439: /*@C
440:    PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
441:    Jacobi preconditioner.

443:    Not Collective

445:    Input Parameter:
446: .  pc - the preconditioner context

448:    Output parameters:
449: +  blocks - the number of blocks
450: -  lens - integer array containing the size of each block

452:    Level: intermediate

454: .keywords:  get, number, Jacobi, global, total, blocks

456: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
457: @*/
458: PetscErrorCode  PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
459: {

465:   PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
466:   return(0);
467: }

471: /*@
472:    PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
473:    Jacobi preconditioner.

475:    Not Collective

477:    Input Parameters:
478: +  pc - the preconditioner context
479: .  blocks - the number of blocks
480: -  lens - [optional] integer array containing size of each block

482:    Note:
483:    Currently only a limited number of blocking configurations are supported.

485:    Level: intermediate

487: .keywords: PC, set, number, Jacobi, local, blocks

489: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
490: @*/
491: PetscErrorCode  PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
492: {

497:   if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
498:   PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
499:   return(0);
500: }

504: /*@C
505:    PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
506:    Jacobi preconditioner.

508:    Not Collective

510:    Input Parameters:
511: +  pc - the preconditioner context
512: .  blocks - the number of blocks
513: -  lens - [optional] integer array containing size of each block

515:    Note:
516:    Currently only a limited number of blocking configurations are supported.

518:    Level: intermediate

520: .keywords: PC, get, number, Jacobi, local, blocks

522: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
523: @*/
524: PetscErrorCode  PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
525: {

531:   PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
532:   return(0);
533: }

535: /* -----------------------------------------------------------------------------------*/

537: /*MC
538:    PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
539:            its own KSP object.

541:    Options Database Keys:
542: .  -pc_use_amat - use Amat to apply block of operator in inner Krylov method

544:    Notes: Each processor can have one or more blocks, but a block cannot be shared by more
545:      than one processor. Defaults to one block per processor.

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

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

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

558:    Level: beginner

560:    Concepts: block Jacobi


563: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
564:            PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
565:            PCBJacobiSetLocalBlocks(), PCSetModifySubmatrices()
566: M*/

570: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
571: {
573:   PetscMPIInt    rank;
574:   PC_BJacobi     *jac;

577:   PetscNewLog(pc,&jac);
578:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);

580:   pc->ops->apply           = 0;
581:   pc->ops->applytranspose  = 0;
582:   pc->ops->setup           = PCSetUp_BJacobi;
583:   pc->ops->destroy         = PCDestroy_BJacobi;
584:   pc->ops->setfromoptions  = PCSetFromOptions_BJacobi;
585:   pc->ops->view            = PCView_BJacobi;
586:   pc->ops->applyrichardson = 0;

588:   pc->data               = (void*)jac;
589:   jac->n                 = -1;
590:   jac->n_local           = -1;
591:   jac->first_local       = rank;
592:   jac->ksp               = 0;
593:   jac->same_local_solves = PETSC_TRUE;
594:   jac->g_lens            = 0;
595:   jac->l_lens            = 0;
596:   jac->psubcomm          = 0;

598:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
599:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
600:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
601:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
602:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
603:   return(0);
604: }

606: /* --------------------------------------------------------------------------------------------*/
607: /*
608:         These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
609: */
612: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
613: {
614:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
615:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
616:   PetscErrorCode         ierr;

619:   KSPReset(jac->ksp[0]);
620:   VecDestroy(&bjac->x);
621:   VecDestroy(&bjac->y);
622:   return(0);
623: }

627: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
628: {
629:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
630:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
631:   PetscErrorCode         ierr;

634:   PCReset_BJacobi_Singleblock(pc);
635:   KSPDestroy(&jac->ksp[0]);
636:   PetscFree(jac->ksp);
637:   PetscFree(jac->l_lens);
638:   PetscFree(jac->g_lens);
639:   PetscFree(bjac);
640:   PetscFree(pc->data);
641:   return(0);
642: }

646: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
647: {
648:   PetscErrorCode     ierr;
649:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
650:   KSP                subksp = jac->ksp[0];
651:   KSPConvergedReason reason;

654:   KSPSetUp(subksp);
655:   KSPGetConvergedReason(subksp,&reason);
656:   if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
657:     pc->failedreason = PC_SUBPC_ERROR;
658:   }
659:   return(0);
660: }

664: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
665: {
666:   PetscErrorCode         ierr;
667:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
668:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;

671:   VecGetLocalVectorRead(x, bjac->x);
672:   VecGetLocalVector(y, bjac->y);
673:  /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
674:      matrix may change even if the outter KSP/PC has not updated the preconditioner, this will trigger a rebuild
675:      of the inner preconditioner automatically unless we pass down the outter preconditioners reuse flag.*/
676:   KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
677:   KSPSolve(jac->ksp[0],bjac->x,bjac->y);
678:   VecRestoreLocalVectorRead(x, bjac->x);
679:   VecRestoreLocalVector(y, bjac->y);
680:   return(0);
681: }

685: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
686: {
687:   PetscErrorCode         ierr;
688:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
689:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
690:   PetscScalar            *y_array;
691:   const PetscScalar      *x_array;
692:   PC                     subpc;

695:   /*
696:       The VecPlaceArray() is to avoid having to copy the
697:     y vector into the bjac->x vector. The reason for
698:     the bjac->x vector is that we need a sequential vector
699:     for the sequential solve.
700:   */
701:   VecGetArrayRead(x,&x_array);
702:   VecGetArray(y,&y_array);
703:   VecPlaceArray(bjac->x,x_array);
704:   VecPlaceArray(bjac->y,y_array);
705:   /* apply the symmetric left portion of the inner PC operator */
706:   /* note this by-passes the inner KSP and its options completely */
707:   KSPGetPC(jac->ksp[0],&subpc);
708:   PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
709:   VecResetArray(bjac->x);
710:   VecResetArray(bjac->y);
711:   VecRestoreArrayRead(x,&x_array);
712:   VecRestoreArray(y,&y_array);
713:   return(0);
714: }

718: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
719: {
720:   PetscErrorCode         ierr;
721:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
722:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
723:   PetscScalar            *y_array;
724:   const PetscScalar      *x_array;
725:   PC                     subpc;

728:   /*
729:       The VecPlaceArray() is to avoid having to copy the
730:     y vector into the bjac->x vector. The reason for
731:     the bjac->x vector is that we need a sequential vector
732:     for the sequential solve.
733:   */
734:   VecGetArrayRead(x,&x_array);
735:   VecGetArray(y,&y_array);
736:   VecPlaceArray(bjac->x,x_array);
737:   VecPlaceArray(bjac->y,y_array);

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

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

745:   VecRestoreArrayRead(x,&x_array);
746:   VecRestoreArray(y,&y_array);
747:   return(0);
748: }

752: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
753: {
754:   PetscErrorCode         ierr;
755:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
756:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
757:   PetscScalar            *y_array;
758:   const PetscScalar      *x_array;

761:   /*
762:       The VecPlaceArray() is to avoid having to copy the
763:     y vector into the bjac->x vector. The reason for
764:     the bjac->x vector is that we need a sequential vector
765:     for the sequential solve.
766:   */
767:   VecGetArrayRead(x,&x_array);
768:   VecGetArray(y,&y_array);
769:   VecPlaceArray(bjac->x,x_array);
770:   VecPlaceArray(bjac->y,y_array);
771:   KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
772:   VecResetArray(bjac->x);
773:   VecResetArray(bjac->y);
774:   VecRestoreArrayRead(x,&x_array);
775:   VecRestoreArray(y,&y_array);
776:   return(0);
777: }

781: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
782: {
783:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
784:   PetscErrorCode         ierr;
785:   PetscInt               m;
786:   KSP                    ksp;
787:   PC_BJacobi_Singleblock *bjac;
788:   PetscBool              wasSetup = PETSC_TRUE;

791:   if (!pc->setupcalled) {
792:     const char *prefix;

794:     if (!jac->ksp) {
795:       wasSetup = PETSC_FALSE;

797:       KSPCreate(PETSC_COMM_SELF,&ksp);
798:       KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
799:       PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
800:       PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
801:       KSPSetType(ksp,KSPPREONLY);
802:       PCGetOptionsPrefix(pc,&prefix);
803:       KSPSetOptionsPrefix(ksp,prefix);
804:       KSPAppendOptionsPrefix(ksp,"sub_");

806:       pc->ops->reset               = PCReset_BJacobi_Singleblock;
807:       pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
808:       pc->ops->apply               = PCApply_BJacobi_Singleblock;
809:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
810:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
811:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
812:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

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

817:       PetscNewLog(pc,&bjac);
818:       jac->data = (void*)bjac;
819:     } else {
820:       ksp  = jac->ksp[0];
821:       bjac = (PC_BJacobi_Singleblock*)jac->data;
822:     }

824:     /*
825:       The reason we need to generate these vectors is to serve
826:       as the right-hand side and solution vector for the solve on the
827:       block. We do not need to allocate space for the vectors since
828:       that is provided via VecPlaceArray() just before the call to
829:       KSPSolve() on the block.
830:     */
831:     MatGetSize(pmat,&m,&m);
832:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
833:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
834: #if defined(PETSC_HAVE_CUSP)
835:     VecSetType(bjac->x,VECCUSP);
836:     VecSetType(bjac->y,VECCUSP);
837: #elif defined(PETSC_HAVE_VECCUDA)
838:     VecSetType(bjac->x,VECCUDA);
839:     VecSetType(bjac->y,VECCUDA);
840: #elif defined(PETSC_HAVE_VIENNACL)
841:     VecSetType(bjac->x,VECVIENNACL);
842:     VecSetType(bjac->y,VECVIENNACL);
843: #endif
844:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
845:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
846:   } else {
847:     ksp  = jac->ksp[0];
848:     bjac = (PC_BJacobi_Singleblock*)jac->data;
849:   }
850:   if (pc->useAmat) {
851:     KSPSetOperators(ksp,mat,pmat);
852:   } else {
853:     KSPSetOperators(ksp,pmat,pmat);
854:   }
855:   if (!wasSetup && pc->setfromoptionscalled) {
856:     /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
857:     KSPSetFromOptions(ksp);
858:   }
859:   return(0);
860: }

862: /* ---------------------------------------------------------------------------------------------*/
865: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
866: {
867:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
868:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
869:   PetscErrorCode        ierr;
870:   PetscInt              i;

873:   if (bjac && bjac->pmat) {
874:     MatDestroyMatrices(jac->n_local,&bjac->pmat);
875:     if (pc->useAmat) {
876:       MatDestroyMatrices(jac->n_local,&bjac->mat);
877:     }
878:   }

880:   for (i=0; i<jac->n_local; i++) {
881:     KSPReset(jac->ksp[i]);
882:     if (bjac && bjac->x) {
883:       VecDestroy(&bjac->x[i]);
884:       VecDestroy(&bjac->y[i]);
885:       ISDestroy(&bjac->is[i]);
886:     }
887:   }
888:   PetscFree(jac->l_lens);
889:   PetscFree(jac->g_lens);
890:   return(0);
891: }

895: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
896: {
897:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
898:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
899:   PetscErrorCode        ierr;
900:   PetscInt              i;

903:   PCReset_BJacobi_Multiblock(pc);
904:   if (bjac) {
905:     PetscFree2(bjac->x,bjac->y);
906:     PetscFree(bjac->starts);
907:     PetscFree(bjac->is);
908:   }
909:   PetscFree(jac->data);
910:   for (i=0; i<jac->n_local; i++) {
911:     KSPDestroy(&jac->ksp[i]);
912:   }
913:   PetscFree(jac->ksp);
914:   PetscFree(pc->data);
915:   return(0);
916: }

920: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
921: {
922:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
923:   PetscErrorCode     ierr;
924:   PetscInt           i,n_local = jac->n_local;
925:   KSPConvergedReason reason;

928:   for (i=0; i<n_local; i++) {
929:     KSPSetUp(jac->ksp[i]);
930:     KSPGetConvergedReason(jac->ksp[i],&reason);
931:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
932:       pc->failedreason = PC_SUBPC_ERROR;
933:     }
934:   }
935:   return(0);
936: }

938: /*
939:       Preconditioner for block Jacobi
940: */
943: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
944: {
945:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
946:   PetscErrorCode        ierr;
947:   PetscInt              i,n_local = jac->n_local;
948:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
949:   PetscScalar           *yin;
950:   const PetscScalar     *xin;

953:   VecGetArrayRead(x,&xin);
954:   VecGetArray(y,&yin);
955:   for (i=0; i<n_local; i++) {
956:     /*
957:        To avoid copying the subvector from x into a workspace we instead
958:        make the workspace vector array point to the subpart of the array of
959:        the global vector.
960:     */
961:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
962:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

964:     PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
965:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
966:     PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

968:     VecResetArray(bjac->x[i]);
969:     VecResetArray(bjac->y[i]);
970:   }
971:   VecRestoreArrayRead(x,&xin);
972:   VecRestoreArray(y,&yin);
973:   return(0);
974: }

976: /*
977:       Preconditioner for block Jacobi
978: */
981: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
982: {
983:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
984:   PetscErrorCode        ierr;
985:   PetscInt              i,n_local = jac->n_local;
986:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
987:   PetscScalar           *yin;
988:   const PetscScalar     *xin;

991:   VecGetArrayRead(x,&xin);
992:   VecGetArray(y,&yin);
993:   for (i=0; i<n_local; i++) {
994:     /*
995:        To avoid copying the subvector from x into a workspace we instead
996:        make the workspace vector array point to the subpart of the array of
997:        the global vector.
998:     */
999:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
1000:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

1002:     PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
1003:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
1004:     PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

1006:     VecResetArray(bjac->x[i]);
1007:     VecResetArray(bjac->y[i]);
1008:   }
1009:   VecRestoreArrayRead(x,&xin);
1010:   VecRestoreArray(y,&yin);
1011:   return(0);
1012: }

1016: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
1017: {
1018:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
1019:   PetscErrorCode        ierr;
1020:   PetscInt              m,n_local,N,M,start,i;
1021:   const char            *prefix,*pprefix,*mprefix;
1022:   KSP                   ksp;
1023:   Vec                   x,y;
1024:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
1025:   PC                    subpc;
1026:   IS                    is;
1027:   MatReuse              scall;

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

1032:   n_local = jac->n_local;

1034:   if (pc->useAmat) {
1035:     PetscBool same;
1036:     PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1037:     if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1038:   }

1040:   if (!pc->setupcalled) {
1041:     scall = MAT_INITIAL_MATRIX;

1043:     if (!jac->ksp) {
1044:       pc->ops->reset         = PCReset_BJacobi_Multiblock;
1045:       pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
1046:       pc->ops->apply         = PCApply_BJacobi_Multiblock;
1047:       pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1048:       pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

1050:       PetscNewLog(pc,&bjac);
1051:       PetscMalloc1(n_local,&jac->ksp);
1052:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1053:       PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1054:       PetscMalloc1(n_local,&bjac->starts);
1055:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));

1057:       jac->data = (void*)bjac;
1058:       PetscMalloc1(n_local,&bjac->is);
1059:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));

1061:       for (i=0; i<n_local; i++) {
1062:         KSPCreate(PETSC_COMM_SELF,&ksp);
1063:         KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
1064:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1065:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
1066:         KSPSetType(ksp,KSPPREONLY);
1067:         KSPGetPC(ksp,&subpc);
1068:         PCGetOptionsPrefix(pc,&prefix);
1069:         KSPSetOptionsPrefix(ksp,prefix);
1070:         KSPAppendOptionsPrefix(ksp,"sub_");

1072:         jac->ksp[i] = ksp;
1073:       }
1074:     } else {
1075:       bjac = (PC_BJacobi_Multiblock*)jac->data;
1076:     }

1078:     start = 0;
1079:     for (i=0; i<n_local; i++) {
1080:       m = jac->l_lens[i];
1081:       /*
1082:       The reason we need to generate these vectors is to serve
1083:       as the right-hand side and solution vector for the solve on the
1084:       block. We do not need to allocate space for the vectors since
1085:       that is provided via VecPlaceArray() just before the call to
1086:       KSPSolve() on the block.

1088:       */
1089:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
1090:       VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1091: #if defined(PETSC_HAVE_CUSP)
1092:       VecSetType(x,VECCUSP);
1093:       VecSetType(y,VECCUSP);
1094: #elif defined(PETSC_HAVE_VECCUDA)
1095:       VecSetType(x,VECCUDA);
1096:       VecSetType(y,VECCUDA);
1097: #elif defined(PETSC_HAVE_VIENNACL)
1098:       VecSetType(x,VECVIENNACL);
1099:       VecSetType(y,VECVIENNACL);
1100: #endif
1101:       PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1102:       PetscLogObjectParent((PetscObject)pc,(PetscObject)y);

1104:       bjac->x[i]      = x;
1105:       bjac->y[i]      = y;
1106:       bjac->starts[i] = start;

1108:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1109:       bjac->is[i] = is;
1110:       PetscLogObjectParent((PetscObject)pc,(PetscObject)is);

1112:       start += m;
1113:     }
1114:   } else {
1115:     bjac = (PC_BJacobi_Multiblock*)jac->data;
1116:     /*
1117:        Destroy the blocks from the previous iteration
1118:     */
1119:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1120:       MatDestroyMatrices(n_local,&bjac->pmat);
1121:       if (pc->useAmat) {
1122:         MatDestroyMatrices(n_local,&bjac->mat);
1123:       }
1124:       scall = MAT_INITIAL_MATRIX;
1125:     } else scall = MAT_REUSE_MATRIX;
1126:   }

1128:   MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1129:   if (pc->useAmat) {
1130:     PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1131:     MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1132:   }
1133:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
1134:      different boundary conditions for the submatrices than for the global problem) */
1135:   PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);

1137:   PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1138:   for (i=0; i<n_local; i++) {
1139:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1140:     PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1141:     if (pc->useAmat) {
1142:       PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1143:       PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1144:       KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1145:     } else {
1146:       KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1147:     }
1148:     if (pc->setfromoptionscalled) {
1149:       KSPSetFromOptions(jac->ksp[i]);
1150:     }
1151:   }
1152:   return(0);
1153: }

1155: /* ---------------------------------------------------------------------------------------------*/
1156: /*
1157:       These are for a single block with multiple processes;
1158: */
1161: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1162: {
1163:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1164:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1165:   PetscErrorCode       ierr;

1168:   VecDestroy(&mpjac->ysub);
1169:   VecDestroy(&mpjac->xsub);
1170:   MatDestroy(&mpjac->submats);
1171:   if (jac->ksp) {KSPReset(jac->ksp[0]);}
1172:   return(0);
1173: }

1177: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1178: {
1179:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1180:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1181:   PetscErrorCode       ierr;

1184:   PCReset_BJacobi_Multiproc(pc);
1185:   KSPDestroy(&jac->ksp[0]);
1186:   PetscFree(jac->ksp);
1187:   PetscSubcommDestroy(&mpjac->psubcomm);

1189:   PetscFree(mpjac);
1190:   PetscFree(pc->data);
1191:   return(0);
1192: }

1196: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1197: {
1198:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1199:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1200:   PetscErrorCode       ierr;
1201:   PetscScalar          *yarray;
1202:   const PetscScalar    *xarray;
1203:   KSPConvergedReason   reason;

1206:   /* place x's and y's local arrays into xsub and ysub */
1207:   VecGetArrayRead(x,&xarray);
1208:   VecGetArray(y,&yarray);
1209:   VecPlaceArray(mpjac->xsub,xarray);
1210:   VecPlaceArray(mpjac->ysub,yarray);

1212:   /* apply preconditioner on each matrix block */
1213:   PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1214:   KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1215:   PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1216:   KSPGetConvergedReason(jac->ksp[0],&reason);
1217:   if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1218:     pc->failedreason = PC_SUBPC_ERROR;
1219:   }

1221:   VecResetArray(mpjac->xsub);
1222:   VecResetArray(mpjac->ysub);
1223:   VecRestoreArrayRead(x,&xarray);
1224:   VecRestoreArray(y,&yarray);
1225:   return(0);
1226: }

1230: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1231: {
1232:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1233:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1234:   PetscErrorCode       ierr;
1235:   PetscInt             m,n;
1236:   MPI_Comm             comm,subcomm=0;
1237:   const char           *prefix;
1238:   PetscBool            wasSetup = PETSC_TRUE;

1241:   PetscObjectGetComm((PetscObject)pc,&comm);
1242:   if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1243:   jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1244:   if (!pc->setupcalled) {
1245:     wasSetup  = PETSC_FALSE;
1246:     PetscNewLog(pc,&mpjac);
1247:     jac->data = (void*)mpjac;

1249:     /* initialize datastructure mpjac */
1250:     if (!jac->psubcomm) {
1251:       /* Create default contiguous subcommunicatiors if user does not provide them */
1252:       PetscSubcommCreate(comm,&jac->psubcomm);
1253:       PetscSubcommSetNumber(jac->psubcomm,jac->n);
1254:       PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1255:       PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1256:     }
1257:     mpjac->psubcomm = jac->psubcomm;
1258:     subcomm         = PetscSubcommChild(mpjac->psubcomm);

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

1263:     /* create a new PC that processors in each subcomm have copy of */
1264:     PetscMalloc1(1,&jac->ksp);
1265:     KSPCreate(subcomm,&jac->ksp[0]);
1266:     KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1267:     PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1268:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1269:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1270:     KSPGetPC(jac->ksp[0],&mpjac->pc);

1272:     PCGetOptionsPrefix(pc,&prefix);
1273:     KSPSetOptionsPrefix(jac->ksp[0],prefix);
1274:     KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1275:     /*
1276:       PetscMPIInt rank,subsize,subrank;
1277:       MPI_Comm_rank(comm,&rank);
1278:       MPI_Comm_size(subcomm,&subsize);
1279:       MPI_Comm_rank(subcomm,&subrank);

1281:       MatGetLocalSize(mpjac->submats,&m,NULL);
1282:       MatGetSize(mpjac->submats,&n,NULL);
1283:       PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1284:       PetscSynchronizedFlush(comm,PETSC_STDOUT);
1285:     */

1287:     /* create dummy vectors xsub and ysub */
1288:     MatGetLocalSize(mpjac->submats,&m,&n);
1289:     VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1290:     VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1291: #if defined(PETSC_HAVE_CUSP)
1292:     VecSetType(mpjac->xsub,VECMPICUSP);
1293:     VecSetType(mpjac->ysub,VECMPICUSP);
1294: #elif defined(PETSC_HAVE_VECCUDA)
1295:     VecSetType(mpjac->xsub,VECMPICUDA);
1296:     VecSetType(mpjac->ysub,VECMPICUDA);
1297: #elif defined(PETSC_HAVE_VECVIENNACL)
1298:     VecSetType(mpjac->xsub,VECMPIVIENNACL);
1299:     VecSetType(mpjac->ysub,VECMPIVIENNACL);
1300: #endif
1301:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1302:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);

1304:     pc->ops->reset   = PCReset_BJacobi_Multiproc;
1305:     pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1306:     pc->ops->apply   = PCApply_BJacobi_Multiproc;
1307:   } else { /* pc->setupcalled */
1308:     subcomm = PetscSubcommChild(mpjac->psubcomm);
1309:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1310:       /* destroy old matrix blocks, then get new matrix blocks */
1311:       if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1312:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1313:     } else {
1314:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1315:     }
1316:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1317:   }

1319:   if (!wasSetup && pc->setfromoptionscalled) {
1320:     KSPSetFromOptions(jac->ksp[0]);
1321:   }
1322:   return(0);
1323: }