Actual source code: bjacobi.c

petsc-3.4.4 2014-03-13
  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:     MPI_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:       PetscMalloc(jac->n_local*sizeof(PetscInt),&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:         PetscMalloc(jac->n_local*sizeof(PetscInt),&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:         PetscMalloc(jac->n_local*sizeof(PetscInt),&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:       PetscMalloc(jac->n_local*sizeof(PetscInt),&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:     PetscMalloc(sizeof(PetscInt),&jac->l_lens);
101:     jac->l_lens[0] = M;
102:   } else { /* jac->n > 0 && jac->n_local > 0 */
103:     if (!jac->l_lens) {
104:       PetscMalloc(jac->n_local*sizeof(PetscInt),&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(PC pc)
163: {
164:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
166:   PetscInt       blocks,i;
167:   PetscBool      flg;

170:   PetscOptionsHead("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:         PetscViewerGetSingleton(viewer,&sviewer);
213:         if (!rank) {
214:           PetscViewerASCIIPushTab(viewer);
215:           KSPView(jac->ksp[0],sviewer);
216:           PetscViewerASCIIPopTab(viewer);
217:         }
218:         PetscViewerRestoreSingleton(viewer,&sviewer);
219:       } else if (jac->psubcomm && !jac->psubcomm->color) {
220:         PetscViewerASCIIGetStdout(mpjac->psubcomm->comm,&sviewer);
221:         PetscViewerASCIIPushTab(viewer);
222:         KSPView(*(jac->ksp),sviewer);
223:         PetscViewerASCIIPopTab(viewer);
224:       }
225:     } else {
226:       PetscInt n_global;
227:       MPI_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
228:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
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:       for (i=0; i<n_global; i++) {
234:         PetscViewerGetSingleton(viewer,&sviewer);
235:         if (i < jac->n_local) {
236:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
237:           KSPView(jac->ksp[i],sviewer);
238:           PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
239:         }
240:         PetscViewerRestoreSingleton(viewer,&sviewer);
241:       }
242:       PetscViewerASCIIPopTab(viewer);
243:       PetscViewerFlush(viewer);
244:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
245:     }
246:   } else if (isstring) {
247:     PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
248:     PetscViewerGetSingleton(viewer,&sviewer);
249:     if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
250:     PetscViewerRestoreSingleton(viewer,&sviewer);
251:   } else if (isdraw) {
252:     PetscDraw draw;
253:     char      str[25];
254:     PetscReal x,y,bottom,h;

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

269: /* -------------------------------------------------------------------------------------*/

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:   *ksp                   = jac->ksp;
283:   jac->same_local_solves = PETSC_FALSE;        /* Assume that local solves are now different;
284:                                                   not necessarily true though!  This flag is
285:                                                   used only for PCView_BJacobi() */
286:   return(0);
287: }

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

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

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

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

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

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

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

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

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

352: /* -------------------------------------------------------------------------------------*/

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

360:    Note Collective

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

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

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

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

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

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

382:    Level: advanced

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

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

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

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

404:    Collective on PC

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

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

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

418:    Level: intermediate

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

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

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

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

441:    Not Collective

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

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

450:    Level: intermediate

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

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

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

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

473:    Not Collective

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

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

483:    Level: intermediate

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

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

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

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

506:    Not Collective

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

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

516:    Level: intermediate

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

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

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

533: /* -----------------------------------------------------------------------------------*/

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

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

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

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

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

552:    Level: beginner

554:    Concepts: block Jacobi

556:    Developer Notes: This preconditioner does not currently work with CUDA/CUSP for a couple of reasons.
557:        (1) It creates seq vectors as work vectors that should be cusp
558:        (2) The use of VecPlaceArray() is not handled properly by CUSP (that is it will not know where
559:            the ownership of the vector is so may use wrong values) even if it did know the ownership
560:            it may induce extra copy ups and downs. Satish suggests a VecTransplantArray() to handle two
561:            vectors sharing the same pointer and handling the CUSP side as well instead of VecGetArray()/VecPlaceArray().


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

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

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

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

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

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

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

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

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

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

647: PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
648: {
650:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

653:   KSPSetUp(jac->ksp[0]);
654:   return(0);
655: }

659: 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;
664:   PetscScalar            *x_array,*y_array;

667:   /*
668:       The VecPlaceArray() is to avoid having to copy the
669:     y vector into the bjac->x vector. The reason for
670:     the bjac->x vector is that we need a sequential vector
671:     for the sequential solve.
672:   */
673:   VecGetArray(x,&x_array);
674:   VecGetArray(y,&y_array);
675:   VecPlaceArray(bjac->x,x_array);
676:   VecPlaceArray(bjac->y,y_array);
677:   KSPSolve(jac->ksp[0],bjac->x,bjac->y);
678:   VecResetArray(bjac->x);
679:   VecResetArray(bjac->y);
680:   VecRestoreArray(x,&x_array);
681:   VecRestoreArray(y,&y_array);
682:   return(0);
683: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

812:       PetscMalloc(sizeof(KSP),&jac->ksp);
813:       jac->ksp[0] = ksp;

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

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

850: /* ---------------------------------------------------------------------------------------------*/
853: PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
854: {
855:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
856:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
857:   PetscErrorCode        ierr;
858:   PetscInt              i;

861:   if (bjac && bjac->pmat) {
862:     MatDestroyMatrices(jac->n_local,&bjac->pmat);
863:     if (pc->useAmat) {
864:       MatDestroyMatrices(jac->n_local,&bjac->mat);
865:     }
866:   }

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

883: PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
884: {
885:   PC_BJacobi            *jac  = (PC_BJacobi*)pc->data;
886:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
887:   PetscErrorCode        ierr;
888:   PetscInt              i;

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

908: PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
909: {
910:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
912:   PetscInt       i,n_local = jac->n_local;

915:   for (i=0; i<n_local; i++) {
916:     KSPSetUp(jac->ksp[i]);
917:   }
918:   return(0);
919: }

921: /*
922:       Preconditioner for block Jacobi
923: */
926: PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
927: {
928:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
929:   PetscErrorCode        ierr;
930:   PetscInt              i,n_local = jac->n_local;
931:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
932:   PetscScalar           *xin,*yin;

935:   VecGetArray(x,&xin);
936:   VecGetArray(y,&yin);
937:   for (i=0; i<n_local; i++) {
938:     /*
939:        To avoid copying the subvector from x into a workspace we instead
940:        make the workspace vector array point to the subpart of the array of
941:        the global vector.
942:     */
943:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
944:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

946:     PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
947:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
948:     PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

950:     VecResetArray(bjac->x[i]);
951:     VecResetArray(bjac->y[i]);
952:   }
953:   VecRestoreArray(x,&xin);
954:   VecRestoreArray(y,&yin);
955:   return(0);
956: }

958: /*
959:       Preconditioner for block Jacobi
960: */
963: PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
964: {
965:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
966:   PetscErrorCode        ierr;
967:   PetscInt              i,n_local = jac->n_local;
968:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
969:   PetscScalar           *xin,*yin;

972:   VecGetArray(x,&xin);
973:   VecGetArray(y,&yin);
974:   for (i=0; i<n_local; i++) {
975:     /*
976:        To avoid copying the subvector from x into a workspace we instead
977:        make the workspace vector array point to the subpart of the array of
978:        the global vector.
979:     */
980:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
981:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

983:     PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
984:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
985:     PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);

987:     VecResetArray(bjac->x[i]);
988:     VecResetArray(bjac->y[i]);
989:   }
990:   VecRestoreArray(x,&xin);
991:   VecRestoreArray(y,&yin);
992:   return(0);
993: }

997: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
998: {
999:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
1000:   PetscErrorCode        ierr;
1001:   PetscInt              m,n_local,N,M,start,i;
1002:   const char            *prefix,*pprefix,*mprefix;
1003:   KSP                   ksp;
1004:   Vec                   x,y;
1005:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
1006:   PC                    subpc;
1007:   IS                    is;
1008:   MatReuse              scall;

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

1013:   n_local = jac->n_local;

1015:   if (pc->useAmat) {
1016:     PetscBool same;
1017:     PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1018:     if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1019:   }

1021:   if (!pc->setupcalled) {
1022:     scall = MAT_INITIAL_MATRIX;

1024:     if (!jac->ksp) {
1025:       pc->ops->reset         = PCReset_BJacobi_Multiblock;
1026:       pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
1027:       pc->ops->apply         = PCApply_BJacobi_Multiblock;
1028:       pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1029:       pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

1031:       PetscNewLog(pc,PC_BJacobi_Multiblock,&bjac);
1032:       PetscMalloc(n_local*sizeof(KSP),&jac->ksp);
1033:       PetscLogObjectMemory(pc,sizeof(n_local*sizeof(KSP)));
1034:       PetscMalloc2(n_local,Vec,&bjac->x,n_local,Vec,&bjac->y);
1035:       PetscMalloc(n_local*sizeof(PetscScalar),&bjac->starts);
1036:       PetscLogObjectMemory(pc,sizeof(n_local*sizeof(PetscScalar)));

1038:       jac->data = (void*)bjac;
1039:       PetscMalloc(n_local*sizeof(IS),&bjac->is);
1040:       PetscLogObjectMemory(pc,sizeof(n_local*sizeof(IS)));

1042:       for (i=0; i<n_local; i++) {
1043:         KSPCreate(PETSC_COMM_SELF,&ksp);
1044:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1045:         PetscLogObjectParent(pc,ksp);
1046:         KSPSetType(ksp,KSPPREONLY);
1047:         KSPGetPC(ksp,&subpc);
1048:         PCGetOptionsPrefix(pc,&prefix);
1049:         KSPSetOptionsPrefix(ksp,prefix);
1050:         KSPAppendOptionsPrefix(ksp,"sub_");

1052:         jac->ksp[i] = ksp;
1053:       }
1054:     } else {
1055:       bjac = (PC_BJacobi_Multiblock*)jac->data;
1056:     }

1058:     start = 0;
1059:     for (i=0; i<n_local; i++) {
1060:       m = jac->l_lens[i];
1061:       /*
1062:       The reason we need to generate these vectors is to serve
1063:       as the right-hand side and solution vector for the solve on the
1064:       block. We do not need to allocate space for the vectors since
1065:       that is provided via VecPlaceArray() just before the call to
1066:       KSPSolve() on the block.

1068:       */
1069:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
1070:       VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1071:       PetscLogObjectParent(pc,x);
1072:       PetscLogObjectParent(pc,y);

1074:       bjac->x[i]      = x;
1075:       bjac->y[i]      = y;
1076:       bjac->starts[i] = start;

1078:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1079:       bjac->is[i] = is;
1080:       PetscLogObjectParent(pc,is);

1082:       start += m;
1083:     }
1084:   } else {
1085:     bjac = (PC_BJacobi_Multiblock*)jac->data;
1086:     /*
1087:        Destroy the blocks from the previous iteration
1088:     */
1089:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1090:       MatDestroyMatrices(n_local,&bjac->pmat);
1091:       if (pc->useAmat) {
1092:         MatDestroyMatrices(n_local,&bjac->mat);
1093:       }
1094:       scall = MAT_INITIAL_MATRIX;
1095:     } else scall = MAT_REUSE_MATRIX;
1096:   }

1098:   MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1099:   if (pc->useAmat) {
1100:     PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1101:     MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1102:   }
1103:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
1104:      different boundary conditions for the submatrices than for the global problem) */
1105:   PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);

1107:   PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1108:   for (i=0; i<n_local; i++) {
1109:     PetscLogObjectParent(pc,bjac->pmat[i]);
1110:     PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1111:     if (pc->useAmat) {
1112:       PetscLogObjectParent(pc,bjac->mat[i]);
1113:       PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1114:       KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i],pc->flag);
1115:     } else {
1116:       KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i],pc->flag);
1117:     }
1118:     if (pc->setfromoptionscalled) {
1119:       KSPSetFromOptions(jac->ksp[i]);
1120:     }
1121:   }
1122:   return(0);
1123: }

1125: /* ---------------------------------------------------------------------------------------------*/
1126: /*
1127:       These are for a single block with multiple processes;
1128: */
1131: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1132: {
1133:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1134:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1135:   PetscErrorCode       ierr;

1138:   VecDestroy(&mpjac->ysub);
1139:   VecDestroy(&mpjac->xsub);
1140:   MatDestroy(&mpjac->submats);
1141:   if (jac->ksp) {KSPReset(jac->ksp[0]);}
1142:   return(0);
1143: }

1147: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1148: {
1149:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1150:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1151:   PetscErrorCode       ierr;

1154:   PCReset_BJacobi_Multiproc(pc);
1155:   KSPDestroy(&jac->ksp[0]);
1156:   PetscFree(jac->ksp);
1157:   PetscSubcommDestroy(&mpjac->psubcomm);

1159:   PetscFree(mpjac);
1160:   PetscFree(pc->data);
1161:   return(0);
1162: }

1166: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1167: {
1168:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1169:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1170:   PetscErrorCode       ierr;
1171:   PetscScalar          *xarray,*yarray;

1174:   /* place x's and y's local arrays into xsub and ysub */
1175:   VecGetArray(x,&xarray);
1176:   VecGetArray(y,&yarray);
1177:   VecPlaceArray(mpjac->xsub,xarray);
1178:   VecPlaceArray(mpjac->ysub,yarray);

1180:   /* apply preconditioner on each matrix block */
1181:   PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1182:   KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1183:   PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);

1185:   VecResetArray(mpjac->xsub);
1186:   VecResetArray(mpjac->ysub);
1187:   VecRestoreArray(x,&xarray);
1188:   VecRestoreArray(y,&yarray);
1189:   return(0);
1190: }

1192: #include <petsc-private/matimpl.h>
1195: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1196: {
1197:   PC_BJacobi           *jac   = (PC_BJacobi*)pc->data;
1198:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1199:   PetscErrorCode       ierr;
1200:   PetscInt             m,n;
1201:   MPI_Comm             comm,subcomm=0;
1202:   const char           *prefix;
1203:   PetscBool            wasSetup = PETSC_TRUE;

1206:   PetscObjectGetComm((PetscObject)pc,&comm);
1207:   if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1208:   jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1209:   if (!pc->setupcalled) {
1210:     wasSetup  = PETSC_FALSE;
1211:     PetscNewLog(pc,PC_BJacobi_Multiproc,&mpjac);
1212:     jac->data = (void*)mpjac;

1214:     /* initialize datastructure mpjac */
1215:     if (!jac->psubcomm) {
1216:       /* Create default contiguous subcommunicatiors if user does not provide them */
1217:       PetscSubcommCreate(comm,&jac->psubcomm);
1218:       PetscSubcommSetNumber(jac->psubcomm,jac->n);
1219:       PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1220:       PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
1221:     }
1222:     mpjac->psubcomm = jac->psubcomm;
1223:     subcomm         = mpjac->psubcomm->comm;

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

1229:     /* create a new PC that processors in each subcomm have copy of */
1230:     PetscMalloc(sizeof(KSP),&jac->ksp);
1231:     KSPCreate(subcomm,&jac->ksp[0]);
1232:     PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1233:     PetscLogObjectParent(pc,jac->ksp[0]);
1234:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats,pc->flag);
1235:     KSPGetPC(jac->ksp[0],&mpjac->pc);

1237:     PCGetOptionsPrefix(pc,&prefix);
1238:     KSPSetOptionsPrefix(jac->ksp[0],prefix);
1239:     KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1240:     /*
1241:       PetscMPIInt rank,subsize,subrank;
1242:       MPI_Comm_rank(comm,&rank);
1243:       MPI_Comm_size(subcomm,&subsize);
1244:       MPI_Comm_rank(subcomm,&subrank);

1246:       MatGetLocalSize(mpjac->submats,&m,NULL);
1247:       MatGetSize(mpjac->submats,&n,NULL);
1248:       PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1249:       PetscSynchronizedFlush(comm);
1250:     */

1252:     /* create dummy vectors xsub and ysub */
1253:     MatGetLocalSize(mpjac->submats,&m,&n);
1254:     VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1255:     VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1256:     PetscLogObjectParent(pc,mpjac->xsub);
1257:     PetscLogObjectParent(pc,mpjac->ysub);

1259:     pc->ops->reset   = PCReset_BJacobi_Multiproc;
1260:     pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1261:     pc->ops->apply   = PCApply_BJacobi_Multiproc;
1262:   } else { /* pc->setupcalled */
1263:     subcomm = mpjac->psubcomm->comm;
1264:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1265:       /* destroy old matrix blocks, then get new matrix blocks */
1266:       if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1267:       (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1268:     } else {
1269:       (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1270:     }
1271:     KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats,pc->flag);
1272:   }

1274:   if (!wasSetup && pc->setfromoptionscalled) {
1275:     KSPSetFromOptions(jac->ksp[0]);
1276:   }
1277:   return(0);
1278: }