Actual source code: bjacobi.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: /*
  3:    Defines a block Jacobi preconditioner.
  4: */
  5: #include <petsc/private/pcimpl.h>              /*I "petscpc.h" I*/
  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 (jac->psubcomm && !jac->psubcomm->color) {
220:         PetscViewerASCIIGetStdout(PetscSubcommChild(mpjac->psubcomm),&sviewer);
221:         PetscViewerASCIIPushTab(viewer);
222:         KSPView(*(jac->ksp),sviewer);
223:         PetscViewerASCIIPopTab(viewer);
224:       }
225:     } else {
226:       PetscInt n_global;
227:       MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
228:       PetscViewerASCIIPushSynchronized(viewer);
229:       PetscViewerASCIIPrintf(viewer,"  Local solve info for each block is in the following KSP and PC objects:\n");
230:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
231:                                                 rank,jac->n_local,jac->first_local);
232:       PetscViewerASCIIPushTab(viewer);
233:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
234:       for (i=0; i<jac->n_local; i++) {
235:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
236:         KSPView(jac->ksp[i],sviewer);
237:         PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
238:       }
239:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
240:       PetscViewerASCIIPopTab(viewer);
241:       PetscViewerFlush(viewer);
242:       PetscViewerASCIIPopSynchronized(viewer);
243:     }
244:   } else if (isstring) {
245:     PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
246:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
247:     if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
248:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
249:   } else if (isdraw) {
250:     PetscDraw draw;
251:     char      str[25];
252:     PetscReal x,y,bottom,h;

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

267: /* -------------------------------------------------------------------------------------*/

271: static PetscErrorCode  PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
272: {
273:   PC_BJacobi *jac = (PC_BJacobi*)pc->data;;

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

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

289: static PetscErrorCode  PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
290: {
291:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

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

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

313:   *blocks = jac->n;
314:   if (lens) *lens = jac->g_lens;
315:   return(0);
316: }

320: static PetscErrorCode  PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
321: {
322:   PC_BJacobi     *jac;

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

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

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

345:   *blocks = jac->n_local;
346:   if (lens) *lens = jac->l_lens;
347:   return(0);
348: }

350: /* -------------------------------------------------------------------------------------*/

354: /*@C
355:    PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
356:    this processor.

358:    Note Collective

360:    Input Parameter:
361: .  pc - the preconditioner context

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

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

371:    Currently for some matrix implementations only 1 block per processor
372:    is supported.

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

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

380:    Level: advanced

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

384: .seealso: PCBJacobiGetSubKSP()
385: @*/
386: PetscErrorCode  PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
387: {

392:   PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
393:   return(0);
394: }

398: /*@
399:    PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
400:    Jacobi preconditioner.

402:    Collective on PC

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

409:    Options Database Key:
410: .  -pc_bjacobi_blocks <blocks> - Sets the number of global blocks

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

416:    Level: intermediate

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

420: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
421: @*/
422: PetscErrorCode  PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
423: {

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

435: /*@C
436:    PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
437:    Jacobi preconditioner.

439:    Not Collective

441:    Input Parameter:
442: .  pc - the preconditioner context

444:    Output parameters:
445: +  blocks - the number of blocks
446: -  lens - integer array containing the size of each block

448:    Level: intermediate

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

452: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
453: @*/
454: PetscErrorCode  PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
455: {

461:   PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
462:   return(0);
463: }

467: /*@
468:    PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
469:    Jacobi preconditioner.

471:    Not Collective

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

478:    Note:
479:    Currently only a limited number of blocking configurations are supported.

481:    Level: intermediate

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

485: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
486: @*/
487: PetscErrorCode  PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
488: {

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

500: /*@C
501:    PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
502:    Jacobi preconditioner.

504:    Not Collective

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

511:    Note:
512:    Currently only a limited number of blocking configurations are supported.

514:    Level: intermediate

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

518: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
519: @*/
520: PetscErrorCode  PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
521: {

527:   PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
528:   return(0);
529: }

531: /* -----------------------------------------------------------------------------------*/

533: /*MC
534:    PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
535:            its own KSP object.

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

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

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

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

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

554:    Level: beginner

556:    Concepts: block Jacobi


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

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

573:   PetscNewLog(pc,&jac);
574:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);

576:   pc->ops->apply           = 0;
577:   pc->ops->applytranspose  = 0;
578:   pc->ops->setup           = PCSetUp_BJacobi;
579:   pc->ops->destroy         = PCDestroy_BJacobi;
580:   pc->ops->setfromoptions  = PCSetFromOptions_BJacobi;
581:   pc->ops->view            = PCView_BJacobi;
582:   pc->ops->applyrichardson = 0;

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

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

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

615:   KSPReset(jac->ksp[0]);
616:   VecDestroy(&bjac->x);
617:   VecDestroy(&bjac->y);
618:   return(0);
619: }

623: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
624: {
625:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
626:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
627:   PetscErrorCode         ierr;

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

640: #include <petsc/private/kspimpl.h>  
643: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
644: {
646:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
647:   KSP            subksp = jac->ksp[0];

650:   KSPSetUp(subksp);
651:   if (subksp->reason == KSP_DIVERGED_PCSETUP_FAILED) {
652:     pc->failedreason = PC_SUBPC_ERROR;
653:   }
654:   return(0);
655: }

659: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
660: {
661:   PetscErrorCode         ierr;
662:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
663:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;

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

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

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

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

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

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

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

740:   VecRestoreArrayRead(x,&x_array);
741:   VecRestoreArray(y,&y_array);
742:   return(0);
743: }

747: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
748: {
749:   PetscErrorCode         ierr;
750:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
751:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
752:   PetscScalar            *y_array;
753:   const PetscScalar      *x_array;

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

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

786:   if (!pc->setupcalled) {
787:     const char *prefix;

789:     if (!jac->ksp) {
790:       wasSetup = PETSC_FALSE;

792:       KSPCreate(PETSC_COMM_SELF,&ksp);
793:       KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
794:       PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
795:       PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
796:       KSPSetType(ksp,KSPPREONLY);
797:       PCGetOptionsPrefix(pc,&prefix);
798:       KSPSetOptionsPrefix(ksp,prefix);
799:       KSPAppendOptionsPrefix(ksp,"sub_");

801:       pc->ops->reset               = PCReset_BJacobi_Singleblock;
802:       pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
803:       pc->ops->apply               = PCApply_BJacobi_Singleblock;
804:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
805:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
806:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
807:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

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

812:       PetscNewLog(pc,&bjac);
813:       jac->data = (void*)bjac;
814:     } else {
815:       ksp  = jac->ksp[0];
816:       bjac = (PC_BJacobi_Singleblock*)jac->data;
817:     }

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

854: /* ---------------------------------------------------------------------------------------------*/
857: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
858: {
859:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
860:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
861:   PetscErrorCode        ierr;
862:   PetscInt              i;

865:   if (bjac && bjac->pmat) {
866:     MatDestroyMatrices(jac->n_local,&bjac->pmat);
867:     if (pc->useAmat) {
868:       MatDestroyMatrices(jac->n_local,&bjac->mat);
869:     }
870:   }

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

887: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
888: {
889:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
890:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
891:   PetscErrorCode        ierr;
892:   PetscInt              i;

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

912: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
913: {
914:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
916:   PetscInt       i,n_local = jac->n_local;

919:   for (i=0; i<n_local; i++) {
920:     KSPSetUp(jac->ksp[i]);
921:     if (jac->ksp[i]->reason == KSP_DIVERGED_PCSETUP_FAILED) {
922:       pc->failedreason = PC_SUBPC_ERROR;
923:     }
924:   }
925:   return(0);
926: }

928: /*
929:       Preconditioner for block Jacobi
930: */
933: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
934: {
935:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
936:   PetscErrorCode        ierr;
937:   PetscInt              i,n_local = jac->n_local;
938:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
939:   PetscScalar           *yin;
940:   const PetscScalar     *xin;

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

954:     PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
955:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
956:     PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

958:     VecResetArray(bjac->x[i]);
959:     VecResetArray(bjac->y[i]);
960:   }
961:   VecRestoreArrayRead(x,&xin);
962:   VecRestoreArray(y,&yin);
963:   return(0);
964: }

966: /*
967:       Preconditioner for block Jacobi
968: */
971: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
972: {
973:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
974:   PetscErrorCode        ierr;
975:   PetscInt              i,n_local = jac->n_local;
976:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
977:   PetscScalar           *yin;
978:   const PetscScalar     *xin;

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

992:     PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
993:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
994:     PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

996:     VecResetArray(bjac->x[i]);
997:     VecResetArray(bjac->y[i]);
998:   }
999:   VecRestoreArrayRead(x,&xin);
1000:   VecRestoreArray(y,&yin);
1001:   return(0);
1002: }

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

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

1022:   n_local = jac->n_local;

1024:   if (pc->useAmat) {
1025:     PetscBool same;
1026:     PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1027:     if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1028:   }

1030:   if (!pc->setupcalled) {
1031:     scall = MAT_INITIAL_MATRIX;

1033:     if (!jac->ksp) {
1034:       pc->ops->reset         = PCReset_BJacobi_Multiblock;
1035:       pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
1036:       pc->ops->apply         = PCApply_BJacobi_Multiblock;
1037:       pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1038:       pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

1040:       PetscNewLog(pc,&bjac);
1041:       PetscMalloc1(n_local,&jac->ksp);
1042:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1043:       PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1044:       PetscMalloc1(n_local,&bjac->starts);
1045:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));

1047:       jac->data = (void*)bjac;
1048:       PetscMalloc1(n_local,&bjac->is);
1049:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));

1051:       for (i=0; i<n_local; i++) {
1052:         KSPCreate(PETSC_COMM_SELF,&ksp);
1053:         KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
1054:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1055:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
1056:         KSPSetType(ksp,KSPPREONLY);
1057:         KSPGetPC(ksp,&subpc);
1058:         PCGetOptionsPrefix(pc,&prefix);
1059:         KSPSetOptionsPrefix(ksp,prefix);
1060:         KSPAppendOptionsPrefix(ksp,"sub_");

1062:         jac->ksp[i] = ksp;
1063:       }
1064:     } else {
1065:       bjac = (PC_BJacobi_Multiblock*)jac->data;
1066:     }

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

1078:       */
1079:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
1080:       VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1081: #if defined(PETSC_HAVE_CUSP)
1082:       VecSetType(x,VECCUSP);
1083:       VecSetType(y,VECCUSP);
1084: #elif defined(PETSC_HAVE_VECCUDA)
1085:       VecSetType(x,VECCUDA);
1086:       VecSetType(y,VECCUDA);
1087: #endif
1088:       PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1089:       PetscLogObjectParent((PetscObject)pc,(PetscObject)y);

1091:       bjac->x[i]      = x;
1092:       bjac->y[i]      = y;
1093:       bjac->starts[i] = start;

1095:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1096:       bjac->is[i] = is;
1097:       PetscLogObjectParent((PetscObject)pc,(PetscObject)is);

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

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

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

1142: /* ---------------------------------------------------------------------------------------------*/
1143: /*
1144:       These are for a single block with multiple processes;
1145: */
1148: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1149: {
1150:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1151:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1152:   PetscErrorCode       ierr;

1155:   VecDestroy(&mpjac->ysub);
1156:   VecDestroy(&mpjac->xsub);
1157:   MatDestroy(&mpjac->submats);
1158:   if (jac->ksp) {KSPReset(jac->ksp[0]);}
1159:   return(0);
1160: }

1164: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1165: {
1166:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1167:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1168:   PetscErrorCode       ierr;

1171:   PCReset_BJacobi_Multiproc(pc);
1172:   KSPDestroy(&jac->ksp[0]);
1173:   PetscFree(jac->ksp);
1174:   PetscSubcommDestroy(&mpjac->psubcomm);

1176:   PetscFree(mpjac);
1177:   PetscFree(pc->data);
1178:   return(0);
1179: }

1183: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1184: {
1185:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1186:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1187:   PetscErrorCode       ierr;
1188:   PetscScalar          *yarray;
1189:   const PetscScalar    *xarray;

1192:   /* place x's and y's local arrays into xsub and ysub */
1193:   VecGetArrayRead(x,&xarray);
1194:   VecGetArray(y,&yarray);
1195:   VecPlaceArray(mpjac->xsub,xarray);
1196:   VecPlaceArray(mpjac->ysub,yarray);

1198:   /* apply preconditioner on each matrix block */
1199:   PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1200:   KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1201:   PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);

1203:   if (jac->ksp[0]->reason == KSP_DIVERGED_PCSETUP_FAILED) {
1204:     pc->failedreason = PC_SUBPC_ERROR;
1205:   }

1207:   VecResetArray(mpjac->xsub);
1208:   VecResetArray(mpjac->ysub);
1209:   VecRestoreArrayRead(x,&xarray);
1210:   VecRestoreArray(y,&yarray);
1211:   return(0);
1212: }

1214: #include <petsc/private/matimpl.h>
1217: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1218: {
1219:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1220:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1221:   PetscErrorCode       ierr;
1222:   PetscInt             m,n;
1223:   MPI_Comm             comm,subcomm=0;
1224:   const char           *prefix;
1225:   PetscBool            wasSetup = PETSC_TRUE;

1228:   PetscObjectGetComm((PetscObject)pc,&comm);
1229:   if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1230:   jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1231:   if (!pc->setupcalled) {
1232:     wasSetup  = PETSC_FALSE;
1233:     PetscNewLog(pc,&mpjac);
1234:     jac->data = (void*)mpjac;

1236:     /* initialize datastructure mpjac */
1237:     if (!jac->psubcomm) {
1238:       /* Create default contiguous subcommunicatiors if user does not provide them */
1239:       PetscSubcommCreate(comm,&jac->psubcomm);
1240:       PetscSubcommSetNumber(jac->psubcomm,jac->n);
1241:       PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1242:       PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1243:     }
1244:     mpjac->psubcomm = jac->psubcomm;
1245:     subcomm         = PetscSubcommChild(mpjac->psubcomm);

1247:     /* Get matrix blocks of pmat */
1248:     if (!pc->pmat->ops->getmultiprocblock) SETERRQ(PetscObjectComm((PetscObject)pc->pmat),PETSC_ERR_SUP,"No support for the requested operation");
1249:     (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);

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

1260:     PCGetOptionsPrefix(pc,&prefix);
1261:     KSPSetOptionsPrefix(jac->ksp[0],prefix);
1262:     KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1263:     /*
1264:       PetscMPIInt rank,subsize,subrank;
1265:       MPI_Comm_rank(comm,&rank);
1266:       MPI_Comm_size(subcomm,&subsize);
1267:       MPI_Comm_rank(subcomm,&subrank);

1269:       MatGetLocalSize(mpjac->submats,&m,NULL);
1270:       MatGetSize(mpjac->submats,&n,NULL);
1271:       PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1272:       PetscSynchronizedFlush(comm,PETSC_STDOUT);
1273:     */

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

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

1304:   if (!wasSetup && pc->setfromoptionscalled) {
1305:     KSPSetFromOptions(jac->ksp[0]);
1306:   }
1307:   return(0);
1308: }