Actual source code: bjacobi.c

      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: }