Actual source code: bjacobi.c

petsc-main 2021-04-20
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);

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

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

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

 33:   /* --------------------------------------------------------------------------
 34:       Determines the number of blocks assigned to each processor
 35:   -----------------------------------------------------------------------------*/

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

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

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

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

143:   PetscFree(jac->g_lens);
144:   PetscFree(jac->l_lens);
145:   PetscFree(pc->data);
146:   return(0);
147: }


150: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
151: {
152:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
154:   PetscInt       blocks,i;
155:   PetscBool      flg;

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

174: #include <petscdraw.h>
175: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
176: {
177:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
178:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
179:   PetscErrorCode       ierr;
180:   PetscMPIInt          rank;
181:   PetscInt             i;
182:   PetscBool            iascii,isstring,isdraw;
183:   PetscViewer          sviewer;
184:   PetscViewerFormat    format;
185:   const char           *prefix;

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

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

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

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

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

280:   if (n_local) *n_local = jac->n_local;
281:   if (first_local) *first_local = jac->first_local;
282:   if (ksp) *ksp                 = jac->ksp;
283:   return(0);
284: }

286: static PetscErrorCode  PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
287: {
288:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

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

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

308:   *blocks = jac->n;
309:   if (lens) *lens = jac->g_lens;
310:   return(0);
311: }

313: static PetscErrorCode  PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
314: {
315:   PC_BJacobi     *jac;

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

321:   jac->n_local = blocks;
322:   if (!lens) jac->l_lens = NULL;
323:   else {
324:     PetscMalloc1(blocks,&jac->l_lens);
325:     PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
326:     PetscArraycpy(jac->l_lens,lens,blocks);
327:   }
328:   return(0);
329: }

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

336:   *blocks = jac->n_local;
337:   if (lens) *lens = jac->l_lens;
338:   return(0);
339: }

341: /* -------------------------------------------------------------------------------------*/

343: /*@C
344:    PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
345:    this processor.

347:    Not Collective

349:    Input Parameter:
350: .  pc - the preconditioner context

352:    Output Parameters:
353: +  n_local - the number of blocks on this processor, or NULL
354: .  first_local - the global number of the first block on this processor, or NULL
355: -  ksp - the array of KSP contexts

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

360:    Currently for some matrix implementations only 1 block per processor
361:    is supported.

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

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

369:    Level: advanced

371: .seealso: PCASMGetSubKSP()
372: @*/
373: PetscErrorCode  PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
374: {

379:   PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
380:   return(0);
381: }

383: /*@
384:    PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
385:    Jacobi preconditioner.

387:    Collective on PC

389:    Input Parameters:
390: +  pc - the preconditioner context
391: .  blocks - the number of blocks
392: -  lens - [optional] integer array containing the size of each block

394:    Options Database Key:
395: .  -pc_bjacobi_blocks <blocks> - Sets the number of global blocks

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

401:    Level: intermediate

403: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
404: @*/
405: PetscErrorCode  PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
406: {

411:   if (blocks <= 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
412:   PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
413:   return(0);
414: }

416: /*@C
417:    PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
418:    Jacobi preconditioner.

420:    Not Collective

422:    Input Parameter:
423: .  pc - the preconditioner context

425:    Output parameters:
426: +  blocks - the number of blocks
427: -  lens - integer array containing the size of each block

429:    Level: intermediate

431: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
432: @*/
433: PetscErrorCode  PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
434: {

440:   PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
441:   return(0);
442: }

444: /*@
445:    PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
446:    Jacobi preconditioner.

448:    Not Collective

450:    Input Parameters:
451: +  pc - the preconditioner context
452: .  blocks - the number of blocks
453: -  lens - [optional] integer array containing size of each block

455:    Options Database Key:
456: .  -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks

458:    Note:
459:    Currently only a limited number of blocking configurations are supported.

461:    Level: intermediate

463: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
464: @*/
465: PetscErrorCode  PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
466: {

471:   if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
472:   PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
473:   return(0);
474: }

476: /*@C
477:    PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
478:    Jacobi preconditioner.

480:    Not Collective

482:    Input Parameters:
483: +  pc - the preconditioner context
484: .  blocks - the number of blocks
485: -  lens - [optional] integer array containing size of each block

487:    Note:
488:    Currently only a limited number of blocking configurations are supported.

490:    Level: intermediate

492: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
493: @*/
494: PetscErrorCode  PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
495: {

501:   PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
502:   return(0);
503: }

505: /* -----------------------------------------------------------------------------------*/

507: /*MC
508:    PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
509:            its own KSP object.

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

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

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

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

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

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

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

533:    Level: beginner

535: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
536:            PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
537:            PCBJacobiSetLocalBlocks(), PCSetModifySubMatrices()
538: M*/

540: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
541: {
543:   PetscMPIInt    rank;
544:   PC_BJacobi     *jac;

547:   PetscNewLog(pc,&jac);
548:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);

550:   pc->ops->apply           = NULL;
551:   pc->ops->matapply        = NULL;
552:   pc->ops->applytranspose  = NULL;
553:   pc->ops->setup           = PCSetUp_BJacobi;
554:   pc->ops->destroy         = PCDestroy_BJacobi;
555:   pc->ops->setfromoptions  = PCSetFromOptions_BJacobi;
556:   pc->ops->view            = PCView_BJacobi;
557:   pc->ops->applyrichardson = NULL;

559:   pc->data               = (void*)jac;
560:   jac->n                 = -1;
561:   jac->n_local           = -1;
562:   jac->first_local       = rank;
563:   jac->ksp               = NULL;
564:   jac->g_lens            = NULL;
565:   jac->l_lens            = NULL;
566:   jac->psubcomm          = NULL;

568:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
569:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
570:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
571:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
572:   PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
573:   return(0);
574: }

576: /* --------------------------------------------------------------------------------------------*/
577: /*
578:         These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
579: */
580: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
581: {
582:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
583:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
584:   PetscErrorCode         ierr;

587:   KSPReset(jac->ksp[0]);
588:   VecDestroy(&bjac->x);
589:   VecDestroy(&bjac->y);
590:   return(0);
591: }

593: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
594: {
595:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
596:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
597:   PetscErrorCode         ierr;

600:   PCReset_BJacobi_Singleblock(pc);
601:   KSPDestroy(&jac->ksp[0]);
602:   PetscFree(jac->ksp);
603:   PetscFree(jac->l_lens);
604:   PetscFree(jac->g_lens);
605:   PetscFree(bjac);
606:   PetscFree(pc->data);
607:   return(0);
608: }

610: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
611: {
612:   PetscErrorCode     ierr;
613:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
614:   KSP                subksp = jac->ksp[0];
615:   KSPConvergedReason reason;

618:   KSPSetUp(subksp);
619:   KSPGetConvergedReason(subksp,&reason);
620:   if (reason == KSP_DIVERGED_PC_FAILED) {
621:     pc->failedreason = PC_SUBPC_ERROR;
622:   }
623:   return(0);
624: }

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

633:   VecGetLocalVectorRead(x, bjac->x);
634:   VecGetLocalVector(y, bjac->y);
635:  /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
636:      matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
637:      of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
638:   KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
639:   KSPSolve(jac->ksp[0],bjac->x,bjac->y);
640:   KSPCheckSolve(jac->ksp[0],pc,bjac->y);
641:   VecRestoreLocalVectorRead(x, bjac->x);
642:   VecRestoreLocalVector(y, bjac->y);
643:   return(0);
644: }

646: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc,Mat X,Mat Y)
647: {
648:   PC_BJacobi     *jac  = (PC_BJacobi*)pc->data;
649:   Mat            sX,sY;

653:  /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
654:      matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
655:      of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
656:   KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
657:   MatDenseGetLocalMatrix(X,&sX);
658:   MatDenseGetLocalMatrix(Y,&sY);
659:   KSPMatSolve(jac->ksp[0],sX,sY);
660:   return(0);
661: }

663: static PetscErrorCode PCApplySymmetricLeft_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;
668:   PetscScalar            *y_array;
669:   const PetscScalar      *x_array;
670:   PC                     subpc;

673:   /*
674:       The VecPlaceArray() is to avoid having to copy the
675:     y vector into the bjac->x vector. The reason for
676:     the bjac->x vector is that we need a sequential vector
677:     for the sequential solve.
678:   */
679:   VecGetArrayRead(x,&x_array);
680:   VecGetArray(y,&y_array);
681:   VecPlaceArray(bjac->x,x_array);
682:   VecPlaceArray(bjac->y,y_array);
683:   /* apply the symmetric left portion of the inner PC operator */
684:   /* note this by-passes the inner KSP and its options completely */
685:   KSPGetPC(jac->ksp[0],&subpc);
686:   PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
687:   VecResetArray(bjac->x);
688:   VecResetArray(bjac->y);
689:   VecRestoreArrayRead(x,&x_array);
690:   VecRestoreArray(y,&y_array);
691:   return(0);
692: }

694: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
695: {
696:   PetscErrorCode         ierr;
697:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
698:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
699:   PetscScalar            *y_array;
700:   const PetscScalar      *x_array;
701:   PC                     subpc;

704:   /*
705:       The VecPlaceArray() is to avoid having to copy the
706:     y vector into the bjac->x vector. The reason for
707:     the bjac->x vector is that we need a sequential vector
708:     for the sequential solve.
709:   */
710:   VecGetArrayRead(x,&x_array);
711:   VecGetArray(y,&y_array);
712:   VecPlaceArray(bjac->x,x_array);
713:   VecPlaceArray(bjac->y,y_array);

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

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

721:   VecRestoreArrayRead(x,&x_array);
722:   VecRestoreArray(y,&y_array);
723:   return(0);
724: }

726: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
727: {
728:   PetscErrorCode         ierr;
729:   PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
730:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
731:   PetscScalar            *y_array;
732:   const PetscScalar      *x_array;

735:   /*
736:       The VecPlaceArray() is to avoid having to copy the
737:     y vector into the bjac->x vector. The reason for
738:     the bjac->x vector is that we need a sequential vector
739:     for the sequential solve.
740:   */
741:   VecGetArrayRead(x,&x_array);
742:   VecGetArray(y,&y_array);
743:   VecPlaceArray(bjac->x,x_array);
744:   VecPlaceArray(bjac->y,y_array);
745:   KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
746:   KSPCheckSolve(jac->ksp[0],pc,bjac->y);
747:   VecResetArray(bjac->x);
748:   VecResetArray(bjac->y);
749:   VecRestoreArrayRead(x,&x_array);
750:   VecRestoreArray(y,&y_array);
751:   return(0);
752: }

754: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
755: {
756:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
757:   PetscErrorCode         ierr;
758:   PetscInt               m;
759:   KSP                    ksp;
760:   PC_BJacobi_Singleblock *bjac;
761:   PetscBool              wasSetup = PETSC_TRUE;
762:   VecType                vectype;
763:   const char             *prefix;

766:   if (!pc->setupcalled) {

768:     if (!jac->ksp) {
769:       wasSetup = PETSC_FALSE;

771:       KSPCreate(PETSC_COMM_SELF,&ksp);
772:       KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
773:       PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
774:       PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
775:       KSPSetType(ksp,KSPPREONLY);
776:       PCGetOptionsPrefix(pc,&prefix);
777:       KSPSetOptionsPrefix(ksp,prefix);
778:       KSPAppendOptionsPrefix(ksp,"sub_");

780:       pc->ops->reset               = PCReset_BJacobi_Singleblock;
781:       pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
782:       pc->ops->apply               = PCApply_BJacobi_Singleblock;
783:       pc->ops->matapply            = PCMatApply_BJacobi_Singleblock;
784:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
785:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
786:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
787:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

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

792:       PetscNewLog(pc,&bjac);
793:       jac->data = (void*)bjac;
794:     } else {
795:       ksp  = jac->ksp[0];
796:       bjac = (PC_BJacobi_Singleblock*)jac->data;
797:     }

799:     /*
800:       The reason we need to generate these vectors is to serve
801:       as the right-hand side and solution vector for the solve on the
802:       block. We do not need to allocate space for the vectors since
803:       that is provided via VecPlaceArray() just before the call to
804:       KSPSolve() on the block.
805:     */
806:     MatGetSize(pmat,&m,&m);
807:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
808:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
809:     MatGetVecType(pmat,&vectype);
810:     VecSetType(bjac->x,vectype);
811:     VecSetType(bjac->y,vectype);
812:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
813:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
814:   } else {
815:     ksp  = jac->ksp[0];
816:     bjac = (PC_BJacobi_Singleblock*)jac->data;
817:   }
818:   KSPGetOptionsPrefix(ksp,&prefix);
819:   if (pc->useAmat) {
820:     KSPSetOperators(ksp,mat,pmat);
821:     MatSetOptionsPrefix(mat,prefix);
822:   } else {
823:     KSPSetOperators(ksp,pmat,pmat);
824:   }
825:   MatSetOptionsPrefix(pmat,prefix);
826:   if (!wasSetup && pc->setfromoptionscalled) {
827:     /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
828:     KSPSetFromOptions(ksp);
829:   }
830:   return(0);
831: }

833: /* ---------------------------------------------------------------------------------------------*/
834: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
835: {
836:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
837:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
838:   PetscErrorCode        ierr;
839:   PetscInt              i;

842:   if (bjac && bjac->pmat) {
843:     MatDestroyMatrices(jac->n_local,&bjac->pmat);
844:     if (pc->useAmat) {
845:       MatDestroyMatrices(jac->n_local,&bjac->mat);
846:     }
847:   }

849:   for (i=0; i<jac->n_local; i++) {
850:     KSPReset(jac->ksp[i]);
851:     if (bjac && bjac->x) {
852:       VecDestroy(&bjac->x[i]);
853:       VecDestroy(&bjac->y[i]);
854:       ISDestroy(&bjac->is[i]);
855:     }
856:   }
857:   PetscFree(jac->l_lens);
858:   PetscFree(jac->g_lens);
859:   return(0);
860: }

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

870:   PCReset_BJacobi_Multiblock(pc);
871:   if (bjac) {
872:     PetscFree2(bjac->x,bjac->y);
873:     PetscFree(bjac->starts);
874:     PetscFree(bjac->is);
875:   }
876:   PetscFree(jac->data);
877:   for (i=0; i<jac->n_local; i++) {
878:     KSPDestroy(&jac->ksp[i]);
879:   }
880:   PetscFree(jac->ksp);
881:   PetscFree(pc->data);
882:   return(0);
883: }

885: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
886: {
887:   PC_BJacobi         *jac = (PC_BJacobi*)pc->data;
888:   PetscErrorCode     ierr;
889:   PetscInt           i,n_local = jac->n_local;
890:   KSPConvergedReason reason;

893:   for (i=0; i<n_local; i++) {
894:     KSPSetUp(jac->ksp[i]);
895:     KSPGetConvergedReason(jac->ksp[i],&reason);
896:     if (reason == KSP_DIVERGED_PC_FAILED) {
897:       pc->failedreason = PC_SUBPC_ERROR;
898:     }
899:   }
900:   return(0);
901: }

903: /*
904:       Preconditioner for block Jacobi
905: */
906: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
907: {
908:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
909:   PetscErrorCode        ierr;
910:   PetscInt              i,n_local = jac->n_local;
911:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
912:   PetscScalar           *yin;
913:   const PetscScalar     *xin;

916:   VecGetArrayRead(x,&xin);
917:   VecGetArray(y,&yin);
918:   for (i=0; i<n_local; i++) {
919:     /*
920:        To avoid copying the subvector from x into a workspace we instead
921:        make the workspace vector array point to the subpart of the array of
922:        the global vector.
923:     */
924:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
925:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

927:     PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
928:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
929:     KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
930:     PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

932:     VecResetArray(bjac->x[i]);
933:     VecResetArray(bjac->y[i]);
934:   }
935:   VecRestoreArrayRead(x,&xin);
936:   VecRestoreArray(y,&yin);
937:   return(0);
938: }

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

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

964:     PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
965:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
966:     KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
967:     PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

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

977: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
978: {
979:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
980:   PetscErrorCode        ierr;
981:   PetscInt              m,n_local,N,M,start,i;
982:   const char            *prefix;
983:   KSP                   ksp;
984:   Vec                   x,y;
985:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
986:   PC                    subpc;
987:   IS                    is;
988:   MatReuse              scall;
989:   VecType               vectype;

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

994:   n_local = jac->n_local;

996:   if (pc->useAmat) {
997:     PetscBool same;
998:     PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
999:     if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1000:   }

1002:   if (!pc->setupcalled) {
1003:     scall = MAT_INITIAL_MATRIX;

1005:     if (!jac->ksp) {
1006:       pc->ops->reset         = PCReset_BJacobi_Multiblock;
1007:       pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
1008:       pc->ops->apply         = PCApply_BJacobi_Multiblock;
1009:       pc->ops->matapply      = NULL;
1010:       pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1011:       pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

1013:       PetscNewLog(pc,&bjac);
1014:       PetscMalloc1(n_local,&jac->ksp);
1015:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1016:       PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1017:       PetscMalloc1(n_local,&bjac->starts);
1018:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));

1020:       jac->data = (void*)bjac;
1021:       PetscMalloc1(n_local,&bjac->is);
1022:       PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));

1024:       for (i=0; i<n_local; i++) {
1025:         KSPCreate(PETSC_COMM_SELF,&ksp);
1026:         KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
1027:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1028:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
1029:         KSPSetType(ksp,KSPPREONLY);
1030:         KSPGetPC(ksp,&subpc);
1031:         PCGetOptionsPrefix(pc,&prefix);
1032:         KSPSetOptionsPrefix(ksp,prefix);
1033:         KSPAppendOptionsPrefix(ksp,"sub_");

1035:         jac->ksp[i] = ksp;
1036:       }
1037:     } else {
1038:       bjac = (PC_BJacobi_Multiblock*)jac->data;
1039:     }

1041:     start = 0;
1042:     MatGetVecType(pmat,&vectype);
1043:     for (i=0; i<n_local; i++) {
1044:       m = jac->l_lens[i];
1045:       /*
1046:       The reason we need to generate these vectors is to serve
1047:       as the right-hand side and solution vector for the solve on the
1048:       block. We do not need to allocate space for the vectors since
1049:       that is provided via VecPlaceArray() just before the call to
1050:       KSPSolve() on the block.

1052:       */
1053:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
1054:       VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1055:       VecSetType(x,vectype);
1056:       VecSetType(y,vectype);
1057:       PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1058:       PetscLogObjectParent((PetscObject)pc,(PetscObject)y);

1060:       bjac->x[i]      = x;
1061:       bjac->y[i]      = y;
1062:       bjac->starts[i] = start;

1064:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1065:       bjac->is[i] = is;
1066:       PetscLogObjectParent((PetscObject)pc,(PetscObject)is);

1068:       start += m;
1069:     }
1070:   } else {
1071:     bjac = (PC_BJacobi_Multiblock*)jac->data;
1072:     /*
1073:        Destroy the blocks from the previous iteration
1074:     */
1075:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1076:       MatDestroyMatrices(n_local,&bjac->pmat);
1077:       if (pc->useAmat) {
1078:         MatDestroyMatrices(n_local,&bjac->mat);
1079:       }
1080:       scall = MAT_INITIAL_MATRIX;
1081:     } else scall = MAT_REUSE_MATRIX;
1082:   }

1084:   MatCreateSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1085:   if (pc->useAmat) {
1086:     MatCreateSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1087:   }
1088:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
1089:      different boundary conditions for the submatrices than for the global problem) */
1090:   PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);

1092:   for (i=0; i<n_local; i++) {
1093:     PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1094:     KSPGetOptionsPrefix(jac->ksp[i],&prefix);
1095:     if (pc->useAmat) {
1096:       PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1097:       KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1098:       MatSetOptionsPrefix(bjac->mat[i],prefix);
1099:     } else {
1100:       KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1101:     }
1102:     MatSetOptionsPrefix(bjac->pmat[i],prefix);
1103:     if (pc->setfromoptionscalled) {
1104:       KSPSetFromOptions(jac->ksp[i]);
1105:     }
1106:   }
1107:   return(0);
1108: }

1110: /* ---------------------------------------------------------------------------------------------*/
1111: /*
1112:       These are for a single block with multiple processes;
1113: */
1114: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1115: {
1116:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1117:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1118:   PetscErrorCode       ierr;

1121:   VecDestroy(&mpjac->ysub);
1122:   VecDestroy(&mpjac->xsub);
1123:   MatDestroy(&mpjac->submats);
1124:   if (jac->ksp) {KSPReset(jac->ksp[0]);}
1125:   return(0);
1126: }

1128: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1129: {
1130:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1131:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1132:   PetscErrorCode       ierr;

1135:   PCReset_BJacobi_Multiproc(pc);
1136:   KSPDestroy(&jac->ksp[0]);
1137:   PetscFree(jac->ksp);
1138:   PetscSubcommDestroy(&mpjac->psubcomm);

1140:   PetscFree(mpjac);
1141:   PetscFree(pc->data);
1142:   return(0);
1143: }

1145: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1146: {
1147:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1148:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1149:   PetscErrorCode       ierr;
1150:   PetscScalar          *yarray;
1151:   const PetscScalar    *xarray;
1152:   KSPConvergedReason   reason;

1155:   /* place x's and y's local arrays into xsub and ysub */
1156:   VecGetArrayRead(x,&xarray);
1157:   VecGetArray(y,&yarray);
1158:   VecPlaceArray(mpjac->xsub,xarray);
1159:   VecPlaceArray(mpjac->ysub,yarray);

1161:   /* apply preconditioner on each matrix block */
1162:   PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1163:   KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1164:   KSPCheckSolve(jac->ksp[0],pc,mpjac->ysub);
1165:   PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1166:   KSPGetConvergedReason(jac->ksp[0],&reason);
1167:   if (reason == KSP_DIVERGED_PC_FAILED) {
1168:     pc->failedreason = PC_SUBPC_ERROR;
1169:   }

1171:   VecResetArray(mpjac->xsub);
1172:   VecResetArray(mpjac->ysub);
1173:   VecRestoreArrayRead(x,&xarray);
1174:   VecRestoreArray(y,&yarray);
1175:   return(0);
1176: }

1178: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc,Mat X,Mat Y)
1179: {
1180:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1181:   KSPConvergedReason   reason;
1182:   Mat                  sX,sY;
1183:   const PetscScalar    *x;
1184:   PetscScalar          *y;
1185:   PetscInt             m,N,lda,ldb;
1186:   PetscErrorCode       ierr;

1189:   /* apply preconditioner on each matrix block */
1190:   MatGetLocalSize(X,&m,NULL);
1191:   MatGetSize(X,NULL,&N);
1192:   MatDenseGetLDA(X,&lda);
1193:   MatDenseGetLDA(Y,&ldb);
1194:   MatDenseGetArrayRead(X,&x);
1195:   MatDenseGetArrayWrite(Y,&y);
1196:   MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,(PetscScalar*)x,&sX);
1197:   MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,y,&sY);
1198:   MatDenseSetLDA(sX,lda);
1199:   MatDenseSetLDA(sY,ldb);
1200:   PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1201:   KSPMatSolve(jac->ksp[0],sX,sY);
1202:   KSPCheckSolve(jac->ksp[0],pc,NULL);
1203:   PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1204:   MatDestroy(&sY);
1205:   MatDestroy(&sX);
1206:   MatDenseRestoreArrayWrite(Y,&y);
1207:   MatDenseRestoreArrayRead(X,&x);
1208:   KSPGetConvergedReason(jac->ksp[0],&reason);
1209:   if (reason == KSP_DIVERGED_PC_FAILED) {
1210:     pc->failedreason = PC_SUBPC_ERROR;
1211:   }
1212:   return(0);
1213: }

1215: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1216: {
1217:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1218:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1219:   PetscErrorCode       ierr;
1220:   PetscInt             m,n;
1221:   MPI_Comm             comm,subcomm=0;
1222:   const char           *prefix;
1223:   PetscBool            wasSetup = PETSC_TRUE;
1224:   VecType              vectype;


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:     MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);

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

1259:     PCGetOptionsPrefix(pc,&prefix);
1260:     KSPSetOptionsPrefix(jac->ksp[0],prefix);
1261:     KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1262:     KSPGetOptionsPrefix(jac->ksp[0],&prefix);
1263:     MatSetOptionsPrefix(mpjac->submats,prefix);

1265:     /* create dummy vectors xsub and ysub */
1266:     MatGetLocalSize(mpjac->submats,&m,&n);
1267:     VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1268:     VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1269:     MatGetVecType(mpjac->submats,&vectype);
1270:     VecSetType(mpjac->xsub,vectype);
1271:     VecSetType(mpjac->ysub,vectype);
1272:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1273:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);

1275:     pc->ops->reset   = PCReset_BJacobi_Multiproc;
1276:     pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1277:     pc->ops->apply   = PCApply_BJacobi_Multiproc;
1278:     pc->ops->matapply= PCMatApply_BJacobi_Multiproc;
1279:   } else { /* pc->setupcalled */
1280:     subcomm = PetscSubcommChild(mpjac->psubcomm);
1281:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1282:       /* destroy old matrix blocks, then get new matrix blocks */
1283:       if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1284:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1285:     } else {
1286:       MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1287:     }
1288:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1289:   }

1291:   if (!wasSetup && pc->setfromoptionscalled) {
1292:     KSPSetFromOptions(jac->ksp[0]);
1293:   }
1294:   return(0);
1295: }