Actual source code: bjacobi.c

petsc-master 2016-08-24
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;
 18:   PetscErrorCode ierr,(*f)(Mat,Mat*);
 19:   PetscInt       N,M,start,i,sum,end;
 20:   PetscInt       bs,i_start=-1,i_end=-1;
 21:   PetscMPIInt    rank,size;
 22:   const char     *pprefix,*mprefix;

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

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

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

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

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

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

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

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


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

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

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

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

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

270: /* -------------------------------------------------------------------------------------*/

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

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

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

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

298:   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");
299:   jac->n = blocks;
300:   if (!lens) jac->g_lens = 0;
301:   else {
302:     PetscMalloc1(blocks,&jac->g_lens);
303:     PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
304:     PetscMemcpy(jac->g_lens,lens,blocks*sizeof(PetscInt));
305:   }
306:   return(0);
307: }

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

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

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

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

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

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

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

353: /* -------------------------------------------------------------------------------------*/

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

361:    Note Collective

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

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

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

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

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

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

383:    Level: advanced

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

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

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

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

405:    Collective on PC

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

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

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

419:    Level: intermediate

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

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

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

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

442:    Not Collective

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

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

451:    Level: intermediate

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

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

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

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

474:    Not Collective

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

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

484:    Level: intermediate

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

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

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

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

507:    Not Collective

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

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

517:    Level: intermediate

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

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

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

534: /* -----------------------------------------------------------------------------------*/

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

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

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

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

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

553:      For CUSP vectors it is recommended to use exactly one block per MPI process for best
554:          performance.  Different block partitioning may lead to additional data transfers
555:          between host and GPU that lead to degraded performance.

557:    Level: beginner

559:    Concepts: block Jacobi


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

858: /* ---------------------------------------------------------------------------------------------*/
861: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
862: {
863:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
864:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
865:   PetscErrorCode        ierr;
866:   PetscInt              i;

869:   if (bjac && bjac->pmat) {
870:     MatDestroyMatrices(jac->n_local,&bjac->pmat);
871:     if (pc->useAmat) {
872:       MatDestroyMatrices(jac->n_local,&bjac->mat);
873:     }
874:   }

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

891: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
892: {
893:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
894:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
895:   PetscErrorCode        ierr;
896:   PetscInt              i;

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

916: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
917: {
918:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
919:   PetscErrorCode     ierr;
920:   PetscInt           i,n_local = jac->n_local;
921:   KSPConvergedReason reason;

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

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

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

960:     PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
961:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
962:     PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

964:     VecResetArray(bjac->x[i]);
965:     VecResetArray(bjac->y[i]);
966:   }
967:   VecRestoreArrayRead(x,&xin);
968:   VecRestoreArray(y,&yin);
969:   return(0);
970: }

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

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

998:     PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
999:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
1000:     PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

1002:     VecResetArray(bjac->x[i]);
1003:     VecResetArray(bjac->y[i]);
1004:   }
1005:   VecRestoreArrayRead(x,&xin);
1006:   VecRestoreArray(y,&yin);
1007:   return(0);
1008: }

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

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

1028:   n_local = jac->n_local;

1030:   if (pc->useAmat) {
1031:     PetscBool same;
1032:     PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1033:     if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1034:   }

1036:   if (!pc->setupcalled) {
1037:     scall = MAT_INITIAL_MATRIX;

1039:     if (!jac->ksp) {
1040:       pc->ops->reset         = PCReset_BJacobi_Multiblock;
1041:       pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
1042:       pc->ops->apply         = PCApply_BJacobi_Multiblock;
1043:       pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1044:       pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

1046:       PetscNewLog(pc,&bjac);
1047:       PetscMalloc1(n_local,&jac->ksp);
1048:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1049:       PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1050:       PetscMalloc1(n_local,&bjac->starts);
1051:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));

1053:       jac->data = (void*)bjac;
1054:       PetscMalloc1(n_local,&bjac->is);
1055:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));

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

1068:         jac->ksp[i] = ksp;
1069:       }
1070:     } else {
1071:       bjac = (PC_BJacobi_Multiblock*)jac->data;
1072:     }

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

1084:       */
1085:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
1086:       VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1087: #if defined(PETSC_HAVE_CUSP)
1088:       VecSetType(x,VECCUSP);
1089:       VecSetType(y,VECCUSP);
1090: #elif defined(PETSC_HAVE_VECCUDA)
1091:       VecSetType(x,VECCUDA);
1092:       VecSetType(y,VECCUDA);
1093: #endif
1094:       PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1095:       PetscLogObjectParent((PetscObject)pc,(PetscObject)y);

1097:       bjac->x[i]      = x;
1098:       bjac->y[i]      = y;
1099:       bjac->starts[i] = start;

1101:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1102:       bjac->is[i] = is;
1103:       PetscLogObjectParent((PetscObject)pc,(PetscObject)is);

1105:       start += m;
1106:     }
1107:   } else {
1108:     bjac = (PC_BJacobi_Multiblock*)jac->data;
1109:     /*
1110:        Destroy the blocks from the previous iteration
1111:     */
1112:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1113:       MatDestroyMatrices(n_local,&bjac->pmat);
1114:       if (pc->useAmat) {
1115:         MatDestroyMatrices(n_local,&bjac->mat);
1116:       }
1117:       scall = MAT_INITIAL_MATRIX;
1118:     } else scall = MAT_REUSE_MATRIX;
1119:   }

1121:   MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1122:   if (pc->useAmat) {
1123:     PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1124:     MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1125:   }
1126:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
1127:      different boundary conditions for the submatrices than for the global problem) */
1128:   PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);

1130:   PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1131:   for (i=0; i<n_local; i++) {
1132:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1133:     PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1134:     if (pc->useAmat) {
1135:       PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1136:       PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1137:       KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1138:     } else {
1139:       KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1140:     }
1141:     if (pc->setfromoptionscalled) {
1142:       KSPSetFromOptions(jac->ksp[i]);
1143:     }
1144:   }
1145:   return(0);
1146: }

1148: /* ---------------------------------------------------------------------------------------------*/
1149: /*
1150:       These are for a single block with multiple processes;
1151: */
1154: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1155: {
1156:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1157:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1158:   PetscErrorCode       ierr;

1161:   VecDestroy(&mpjac->ysub);
1162:   VecDestroy(&mpjac->xsub);
1163:   MatDestroy(&mpjac->submats);
1164:   if (jac->ksp) {KSPReset(jac->ksp[0]);}
1165:   return(0);
1166: }

1170: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1171: {
1172:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1173:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1174:   PetscErrorCode       ierr;

1177:   PCReset_BJacobi_Multiproc(pc);
1178:   KSPDestroy(&jac->ksp[0]);
1179:   PetscFree(jac->ksp);
1180:   PetscSubcommDestroy(&mpjac->psubcomm);

1182:   PetscFree(mpjac);
1183:   PetscFree(pc->data);
1184:   return(0);
1185: }

1189: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1190: {
1191:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1192:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1193:   PetscErrorCode       ierr;
1194:   PetscScalar          *yarray;
1195:   const PetscScalar    *xarray;
1196:   KSPConvergedReason   reason;

1199:   /* place x's and y's local arrays into xsub and ysub */
1200:   VecGetArrayRead(x,&xarray);
1201:   VecGetArray(y,&yarray);
1202:   VecPlaceArray(mpjac->xsub,xarray);
1203:   VecPlaceArray(mpjac->ysub,yarray);

1205:   /* apply preconditioner on each matrix block */
1206:   PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1207:   KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1208:   PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1209:   KSPGetConvergedReason(jac->ksp[0],&reason);
1210:   if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1211:     pc->failedreason = PC_SUBPC_ERROR;
1212:   }

1214:   VecResetArray(mpjac->xsub);
1215:   VecResetArray(mpjac->ysub);
1216:   VecRestoreArrayRead(x,&xarray);
1217:   VecRestoreArray(y,&yarray);
1218:   return(0);
1219: }

1223: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1224: {
1225:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1226:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1227:   PetscErrorCode       ierr;
1228:   PetscInt             m,n;
1229:   MPI_Comm             comm,subcomm=0;
1230:   const char           *prefix;
1231:   PetscBool            wasSetup = PETSC_TRUE;

1234:   PetscObjectGetComm((PetscObject)pc,&comm);
1235:   if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1236:   jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1237:   if (!pc->setupcalled) {
1238:     wasSetup  = PETSC_FALSE;
1239:     PetscNewLog(pc,&mpjac);
1240:     jac->data = (void*)mpjac;

1242:     /* initialize datastructure mpjac */
1243:     if (!jac->psubcomm) {
1244:       /* Create default contiguous subcommunicatiors if user does not provide them */
1245:       PetscSubcommCreate(comm,&jac->psubcomm);
1246:       PetscSubcommSetNumber(jac->psubcomm,jac->n);
1247:       PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1248:       PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1249:     }
1250:     mpjac->psubcomm = jac->psubcomm;
1251:     subcomm         = PetscSubcommChild(mpjac->psubcomm);

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

1256:     /* create a new PC that processors in each subcomm have copy of */
1257:     PetscMalloc1(1,&jac->ksp);
1258:     KSPCreate(subcomm,&jac->ksp[0]);
1259:     KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1260:     PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1261:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1262:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1263:     KSPGetPC(jac->ksp[0],&mpjac->pc);

1265:     PCGetOptionsPrefix(pc,&prefix);
1266:     KSPSetOptionsPrefix(jac->ksp[0],prefix);
1267:     KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1268:     /*
1269:       PetscMPIInt rank,subsize,subrank;
1270:       MPI_Comm_rank(comm,&rank);
1271:       MPI_Comm_size(subcomm,&subsize);
1272:       MPI_Comm_rank(subcomm,&subrank);

1274:       MatGetLocalSize(mpjac->submats,&m,NULL);
1275:       MatGetSize(mpjac->submats,&n,NULL);
1276:       PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1277:       PetscSynchronizedFlush(comm,PETSC_STDOUT);
1278:     */

1280:     /* create dummy vectors xsub and ysub */
1281:     MatGetLocalSize(mpjac->submats,&m,&n);
1282:     VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1283:     VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1284: #if defined(PETSC_HAVE_CUSP)
1285:     VecSetType(mpjac->xsub,VECMPICUSP);
1286:     VecSetType(mpjac->ysub,VECMPICUSP);
1287: #elif defined(PETSC_HAVE_VECCUDA)
1288:     VecSetType(mpjac->xsub,VECMPICUDA);
1289:     VecSetType(mpjac->ysub,VECMPICUDA);
1290: #endif
1291:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1292:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);

1294:     pc->ops->reset   = PCReset_BJacobi_Multiproc;
1295:     pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1296:     pc->ops->apply   = PCApply_BJacobi_Multiproc;
1297:   } else { /* pc->setupcalled */
1298:     subcomm = PetscSubcommChild(mpjac->psubcomm);
1299:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1300:       /* destroy old matrix blocks, then get new matrix blocks */
1301:       if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1302:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1303:     } else {
1304:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1305:     }
1306:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1307:   }

1309:   if (!wasSetup && pc->setfromoptionscalled) {
1310:     KSPSetFromOptions(jac->ksp[0]);
1311:   }
1312:   return(0);
1313: }