Actual source code: bjacobi.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Defines a block Jacobi preconditioner.
      5: */
     6:  #include src/mat/matimpl.h
     7:  #include src/ksp/pc/pcimpl.h
     8:  #include src/ksp/pc/impls/bjacobi/bjacobi.h

     10: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC,Mat,Mat);
     11: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC,Mat,Mat);

     15: static PetscErrorCode PCSetUp_BJacobi(PC pc)
     16: {
     17:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
     18:   Mat            mat = pc->mat,pmat = pc->pmat;
     19:   PetscErrorCode ierr,(*f)(Mat,PetscTruth*,MatReuse,Mat*);
     20:   PetscInt       N,M,start,i,sum,end;
     21:   PetscInt       bs,i_start=-1,i_end=-1;
     22:   PetscMPIInt    rank,size;
     23:   const char     *pprefix,*mprefix;

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

     31:   /* ----------
     32:       Determines the number of blocks assigned to each processor
     33:   */

     35:   /*   local block count  given */
     36:   if (jac->n_local > 0 && jac->n < 0) {
     37:     MPI_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,pc->comm);
     38:     if (jac->l_lens) { /* check that user set these correctly */
     39:       sum = 0;
     40:       for (i=0; i<jac->n_local; i++) {
     41:         if (jac->l_lens[i]/bs*bs !=jac->l_lens[i]) {
     42:           SETERRQ(PETSC_ERR_ARG_SIZ,"Mat blocksize doesn't match block Jacobi layout");
     43:         }
     44:         sum += jac->l_lens[i];
     45:       }
     46:       if (sum != M) SETERRQ(PETSC_ERR_ARG_SIZ,"Local lens sent incorrectly");
     47:     } else {
     48:       PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
     49:       for (i=0; i<jac->n_local; i++) {
     50:         jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
     51:       }
     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_ERR_ARG_SIZ,"Zero block not allowed");
     59:         if (jac->g_lens[i]/bs*bs != jac->g_lens[i]) {
     60:           SETERRQ(PETSC_ERR_ARG_SIZ,"Mat blocksize doesn't match block Jacobi layout");
     61:         }
     62:       }
     63:       if (size == 1) {
     64:         jac->n_local = jac->n;
     65:         PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
     66:         PetscMemcpy(jac->l_lens,jac->g_lens,jac->n_local*sizeof(PetscInt));
     67:         /* check that user set these correctly */
     68:         sum = 0;
     69:         for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
     70:         if (sum != M) SETERRQ(PETSC_ERR_ARG_SIZ,"Global lens sent incorrectly");
     71:       } else {
     72:         MatGetOwnershipRange(pc->pmat,&start,&end);
     73:         /* loop over blocks determing first one owned by me */
     74:         sum = 0;
     75:         for (i=0; i<jac->n+1; i++) {
     76:           if (sum == start) { i_start = i; goto start_1;}
     77:           if (i < jac->n) sum += jac->g_lens[i];
     78:         }
     79:         SETERRQ(PETSC_ERR_ARG_SIZ,"Block sizes\n\
     80:                    used in PCBJacobiSetTotalBlocks()\n\
     81:                    are not compatible with parallel matrix layout");
     82:  start_1:
     83:         for (i=i_start; i<jac->n+1; i++) {
     84:           if (sum == end) { i_end = i; goto end_1; }
     85:           if (i < jac->n) sum += jac->g_lens[i];
     86:         }
     87:         SETERRQ(PETSC_ERR_ARG_SIZ,"Block sizes\n\
     88:                       used in PCBJacobiSetTotalBlocks()\n\
     89:                       are not compatible with parallel matrix layout");
     90:  end_1:
     91:         jac->n_local = i_end - i_start;
     92:         PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
     93:         PetscMemcpy(jac->l_lens,jac->g_lens+i_start,jac->n_local*sizeof(PetscInt));
     94:       }
     95:     } else { /* no global blocks given, determine then using default layout */
     96:       jac->n_local = jac->n/size + ((jac->n % size) > rank);
     97:       PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
     98:       for (i=0; i<jac->n_local; i++) {
     99:         jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
    100:         if (!jac->l_lens[i]) SETERRQ(PETSC_ERR_ARG_SIZ,"Too many blocks given");
    101:       }
    102:     }
    103:   } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
    104:     jac->n         = size;
    105:     jac->n_local   = 1;
    106:     PetscMalloc(sizeof(PetscInt),&jac->l_lens);
    107:     jac->l_lens[0] = M;
    108:   }

    110:   MPI_Comm_size(pc->comm,&size);
    111:   PetscObjectQueryFunction((PetscObject)pc->mat,"MatGetDiagonalBlock_C",(void (**)(void))&f);
    112:   if (size == 1 && !f) {
    113:     mat  = pc->mat;
    114:     pmat = pc->pmat;
    115:   } else {
    116:     PetscTruth iscopy;
    117:     MatReuse   scall;

    119:     if (jac->use_true_local) {
    120:       scall = MAT_INITIAL_MATRIX;
    121:       if (pc->setupcalled) {
    122:         if (pc->flag == SAME_NONZERO_PATTERN) {
    123:           if (jac->tp_mat) {
    124:             scall = MAT_REUSE_MATRIX;
    125:             mat   = jac->tp_mat;
    126:           }
    127:         } else {
    128:           if (jac->tp_mat)  {
    129:             MatDestroy(jac->tp_mat);
    130:           }
    131:         }
    132:       }
    133:       if (!f) {
    134:         SETERRQ(PETSC_ERR_SUP,"This matrix does not support getting diagonal block");
    135:       }
    136:       (*f)(pc->mat,&iscopy,scall,&mat);
    137:       /* make submatrix have same prefix as entire matrix */
    138:       PetscObjectGetOptionsPrefix((PetscObject)pc->mat,&mprefix);
    139:       PetscObjectSetOptionsPrefix((PetscObject)mat,mprefix);
    140:       if (iscopy) {
    141:         jac->tp_mat = mat;
    142:       }
    143:     }
    144:     if (pc->pmat != pc->mat || !jac->use_true_local) {
    145:       scall = MAT_INITIAL_MATRIX;
    146:       if (pc->setupcalled) {
    147:         if (pc->flag == SAME_NONZERO_PATTERN) {
    148:           if (jac->tp_pmat) {
    149:             scall = MAT_REUSE_MATRIX;
    150:             pmat   = jac->tp_pmat;
    151:           }
    152:         } else {
    153:           if (jac->tp_pmat)  {
    154:             MatDestroy(jac->tp_pmat);
    155:           }
    156:         }
    157:       }
    158:       PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",(void (**)(void))&f);
    159:       if (!f) {
    160:         SETERRQ(PETSC_ERR_SUP,"This matrix does not support getting diagonal block");
    161:       }
    162:       (*f)(pc->pmat,&iscopy,scall,&pmat);
    163:       /* make submatrix have same prefix as entire matrix */
    164:       PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
    165:       PetscObjectSetOptionsPrefix((PetscObject)pmat,pprefix);
    166:       if (iscopy) {
    167:         jac->tp_pmat = pmat;
    168:       }
    169:     } else {
    170:       pmat = mat;
    171:     }
    172:   }

    174:   /* ------
    175:      Setup code depends on the number of blocks
    176:   */
    177:   if (jac->n_local == 1) {
    178:     PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
    179:   } else {
    180:     PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
    181:   }
    182:   return(0);
    183: }

    185: /* Default destroy, if it has never been setup */
    188: static PetscErrorCode PCDestroy_BJacobi(PC pc)
    189: {
    190:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

    194:   if (jac->g_lens) {PetscFree(jac->g_lens);}
    195:   if (jac->l_lens) {PetscFree(jac->l_lens);}
    196:   PetscFree(jac);
    197:   return(0);
    198: }

    202: static PetscErrorCode PCSetFromOptions_BJacobi(PC pc)
    203: {
    204:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
    206:   PetscInt       blocks;
    207:   PetscTruth     flg;

    210:   PetscOptionsHead("Block Jacobi options");
    211:     PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
    212:     if (flg) {
    213:       PCBJacobiSetTotalBlocks(pc,blocks,PETSC_NULL);
    214:     }
    215:     PetscOptionsName("-pc_bjacobi_truelocal","Use the true matrix, not preconditioner matrix to define matrix vector product in sub-problems","PCBJacobiSetUseTrueLocal",&flg);
    216:     if (flg) {
    217:       PCBJacobiSetUseTrueLocal(pc);
    218:     }
    219:   PetscOptionsTail();
    220:   return(0);
    221: }

    225: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
    226: {
    227:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
    229:   PetscMPIInt    rank;
    230:   PetscInt       i;
    231:   PetscTruth     iascii,isstring;
    232:   PetscViewer    sviewer;

    235:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    236:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
    237:   if (iascii) {
    238:     if (jac->use_true_local) {
    239:       PetscViewerASCIIPrintf(viewer,"  block Jacobi: using true local matrix, number of blocks = %D\n",jac->n);
    240:     }
    241:     PetscViewerASCIIPrintf(viewer,"  block Jacobi: number of blocks = %D\n",jac->n);
    242:     MPI_Comm_rank(pc->comm,&rank);
    243:     if (jac->same_local_solves) {
    244:       PetscViewerASCIIPrintf(viewer,"  Local solve is same for all blocks, in the following KSP and PC objects:\n");
    245:       PetscViewerGetSingleton(viewer,&sviewer);
    246:       if (!rank && jac->ksp) {
    247:         PetscViewerASCIIPushTab(viewer);
    248:         KSPView(jac->ksp[0],sviewer);
    249:         PetscViewerASCIIPopTab(viewer);
    250:       }
    251:       PetscViewerRestoreSingleton(viewer,&sviewer);
    252:     } else {
    253:       PetscInt n_global;
    254:       MPI_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,pc->comm);
    255:       PetscViewerASCIIPrintf(viewer,"  Local solve info for each block is in the following KSP and PC objects:\n");
    256:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
    257:                    rank,jac->n_local,jac->first_local);
    258:       PetscViewerASCIIPushTab(viewer);
    259:       for (i=0; i<n_global; i++) {
    260:         PetscViewerGetSingleton(viewer,&sviewer);
    261:         if (i < jac->n_local) {
    262:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
    263:           KSPView(jac->ksp[i],sviewer);
    264:           PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
    265:         }
    266:         PetscViewerRestoreSingleton(viewer,&sviewer);
    267:       }
    268:       PetscViewerASCIIPopTab(viewer);
    269:       PetscViewerFlush(viewer);
    270:     }
    271:   } else if (isstring) {
    272:     PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
    273:     PetscViewerGetSingleton(viewer,&sviewer);
    274:     if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
    275:     PetscViewerRestoreSingleton(viewer,&sviewer);
    276:   } else {
    277:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for block Jacobi",((PetscObject)viewer)->type_name);
    278:   }
    279:   return(0);
    280: }

    282: /* -------------------------------------------------------------------------------------*/

    287: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiSetUseTrueLocal_BJacobi(PC pc)
    288: {
    289:   PC_BJacobi   *jac;

    292:   jac                 = (PC_BJacobi*)pc->data;
    293:   jac->use_true_local = PETSC_TRUE;
    294:   return(0);
    295: }

    301: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
    302: {
    303:   PC_BJacobi   *jac = (PC_BJacobi*)pc->data;;

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

    308:   if (n_local)     *n_local     = jac->n_local;
    309:   if (first_local) *first_local = jac->first_local;
    310:   *ksp                          = jac->ksp;
    311:   jac->same_local_solves        = PETSC_FALSE; /* Assume that local solves are now different;
    312:                                                   not necessarily true though!  This flag is
    313:                                                   used only for PCView_BJacobi() */
    314:   return(0);
    315: }

    321: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
    322: {
    323:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

    328:   if (pc->setupcalled > 0) SETERRQ(PETSC_ERR_ORDER,"Cannot set number of blocks after PCSetUp()/KSPSetUp() has been called");
    329:   jac->n = blocks;
    330:   if (!lens) {
    331:     jac->g_lens = 0;
    332:   } else {
    333:     PetscMalloc(blocks*sizeof(PetscInt),&jac->g_lens);
    334:     PetscLogObjectMemory(pc,blocks*sizeof(PetscInt));
    335:     PetscMemcpy(jac->g_lens,lens,blocks*sizeof(PetscInt));
    336:   }
    337:   return(0);
    338: }

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

    349:   *blocks = jac->n;
    350:   if (lens) *lens = jac->g_lens;
    351:   return(0);
    352: }

    358: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
    359: {
    360:   PC_BJacobi     *jac;

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

    366:   jac->n_local = blocks;
    367:   if (!lens) {
    368:     jac->l_lens = 0;
    369:   } else {
    370:     PetscMalloc(blocks*sizeof(PetscInt),&jac->l_lens);
    371:     PetscLogObjectMemory(pc,blocks*sizeof(PetscInt));
    372:     PetscMemcpy(jac->l_lens,lens,blocks*sizeof(PetscInt));
    373:   }
    374:   return(0);
    375: }

    381: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
    382: {
    383:   PC_BJacobi *jac = (PC_BJacobi*) pc->data;

    386:   *blocks = jac->n_local;
    387:   if (lens) *lens = jac->l_lens;
    388:   return(0);
    389: }

    392: /* -------------------------------------------------------------------------------------*/

    396: /*@
    397:    PCBJacobiSetUseTrueLocal - Sets a flag to indicate that the block
    398:    problem is associated with the linear system matrix instead of the
    399:    default (where it is associated with the preconditioning matrix).
    400:    That is, if the local system is solved iteratively then it iterates
    401:    on the block from the matrix using the block from the preconditioner
    402:    as the preconditioner for the local block.

    404:    Collective on PC

    406:    Input Parameters:
    407: .  pc - the preconditioner context

    409:    Options Database Key:
    410: .  -pc_bjacobi_truelocal - Activates PCBJacobiSetUseTrueLocal()

    412:    Notes:
    413:    For the common case in which the preconditioning and linear
    414:    system matrices are identical, this routine is unnecessary.

    416:    Level: intermediate

    418: .keywords:  block, Jacobi, set, true, local, flag

    420: .seealso: PCSetOperators(), PCBJacobiSetLocalBlocks()
    421: @*/
    422: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiSetUseTrueLocal(PC pc)
    423: {
    424:   PetscErrorCode ierr,(*f)(PC);

    428:   PetscObjectQueryFunction((PetscObject)pc,"PCBJacobiSetUseTrueLocal_C",(void (**)(void))&f);
    429:   if (f) {
    430:     (*f)(pc);
    431:   }

    433:   return(0);
    434: }

    438: /*@C
    439:    PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
    440:    this processor.
    441:
    442:    Note Collective

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

    447:    Output Parameters:
    448: +  n_local - the number of blocks on this processor, or PETSC_NULL
    449: .  first_local - the global number of the first block on this processor, or PETSC_NULL
    450: -  ksp - the array of KSP contexts

    452:    Notes:
    453:    After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
    454:
    455:    Currently for some matrix implementations only 1 block per processor
    456:    is supported.
    457:
    458:    You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().

    460:    Level: advanced

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

    464: .seealso: PCBJacobiGetSubKSP()
    465: @*/
    466: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
    467: {
    468:   PetscErrorCode ierr,(*f)(PC,PetscInt *,PetscInt *,KSP **);

    472:   PetscObjectQueryFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",(void (**)(void))&f);
    473:   if (f) {
    474:     (*f)(pc,n_local,first_local,ksp);
    475:   } else {
    476:     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subsolvers for this preconditioner");
    477:   }
    478:   return(0);
    479: }

    483: /*@
    484:    PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
    485:    Jacobi preconditioner.

    487:    Collective on PC

    489:    Input Parameters:
    490: +  pc - the preconditioner context
    491: .  blocks - the number of blocks
    492: -  lens - [optional] integer array containing the size of each block

    494:    Options Database Key:
    495: .  -pc_bjacobi_blocks <blocks> - Sets the number of global blocks

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

    501:    Level: intermediate

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

    505: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiSetLocalBlocks()
    506: @*/
    507: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
    508: {
    509:   PetscErrorCode ierr,(*f)(PC,PetscInt,const PetscInt[]);

    513:   if (blocks <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
    514:   PetscObjectQueryFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",(void (**)(void))&f);
    515:   if (f) {
    516:     (*f)(pc,blocks,lens);
    517:   }
    518:   return(0);
    519: }

    523: /*@C
    524:    PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
    525:    Jacobi preconditioner.

    527:    Collective on PC

    529:    Input Parameter:
    530: .  pc - the preconditioner context

    532:    Output parameters:
    533: +  blocks - the number of blocks
    534: -  lens - integer array containing the size of each block

    536:    Level: intermediate

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

    540: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiGetLocalBlocks()
    541: @*/
    542: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
    543: {
    544:   PetscErrorCode ierr,(*f)(PC,PetscInt*, const PetscInt *[]);

    549:   PetscObjectQueryFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",(void (**)(void))&f);
    550:   if (f) {
    551:     (*f)(pc,blocks,lens);
    552:   }
    553:   return(0);
    554: }
    555:
    558: /*@
    559:    PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
    560:    Jacobi preconditioner.

    562:    Not Collective

    564:    Input Parameters:
    565: +  pc - the preconditioner context
    566: .  blocks - the number of blocks
    567: -  lens - [optional] integer array containing size of each block

    569:    Note:
    570:    Currently only a limited number of blocking configurations are supported.

    572:    Level: intermediate

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

    576: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiSetTotalBlocks()
    577: @*/
    578: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
    579: {
    580:   PetscErrorCode ierr,(*f)(PC,PetscInt,const PetscInt []);

    584:   if (blocks < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
    585:   PetscObjectQueryFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",(void (**)(void))&f);
    586:   if (f) {
    587:     (*f)(pc,blocks,lens);
    588:   }
    589:   return(0);
    590: }
    591:
    594: /*@C
    595:    PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
    596:    Jacobi preconditioner.

    598:    Not Collective

    600:    Input Parameters:
    601: +  pc - the preconditioner context
    602: .  blocks - the number of blocks
    603: -  lens - [optional] integer array containing size of each block

    605:    Note:
    606:    Currently only a limited number of blocking configurations are supported.

    608:    Level: intermediate

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

    612: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiGetTotalBlocks()
    613: @*/
    614: PetscErrorCode PETSCKSP_DLLEXPORT PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
    615: {
    616:   PetscErrorCode ierr,(*f)(PC,PetscInt*, const PetscInt *[]);

    621:   PetscObjectQueryFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",(void (**)(void))&f);
    622:   if (f) {
    623:     (*f)(pc,blocks,lens);
    624:   }
    625:   return(0);
    626: }

    628: /* -----------------------------------------------------------------------------------*/

    630: /*MC
    631:    PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
    632:            its own KSP object.

    634:    Options Database Keys:
    635: .  -pc_bjacobi_truelocal - Activates PCBJacobiSetUseTrueLocal()

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

    640:      To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
    641:         options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1 -sub_ksp_type preonly
    642:
    643:      To set the options on the solvers seperate for each block call PCBJacobiGetSubKSP()
    644:          and set the options directly on the resulting KSP object (you can access its PC
    645:          KSPGetPC())

    647:    Level: beginner

    649:    Concepts: block Jacobi

    651: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    652:            PCASM, PCBJacobiSetUseTrueLocal(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
    653:            PCBJacobiSetLocalBlocks(), PCSetModifySubmatrices()
    654: M*/

    659: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_BJacobi(PC pc)
    660: {
    662:   PetscMPIInt    rank;
    663:   PC_BJacobi     *jac;

    666:   PetscNew(PC_BJacobi,&jac);
    667:   PetscLogObjectMemory(pc,sizeof(PC_BJacobi));
    668:   MPI_Comm_rank(pc->comm,&rank);
    669:   pc->ops->apply              = 0;
    670:   pc->ops->applytranspose     = 0;
    671:   pc->ops->setup              = PCSetUp_BJacobi;
    672:   pc->ops->destroy            = PCDestroy_BJacobi;
    673:   pc->ops->setfromoptions     = PCSetFromOptions_BJacobi;
    674:   pc->ops->view               = PCView_BJacobi;
    675:   pc->ops->applyrichardson    = 0;

    677:   pc->data               = (void*)jac;
    678:   jac->n                 = -1;
    679:   jac->n_local           = -1;
    680:   jac->first_local       = rank;
    681:   jac->ksp              = 0;
    682:   jac->use_true_local    = PETSC_FALSE;
    683:   jac->same_local_solves = PETSC_TRUE;
    684:   jac->g_lens            = 0;
    685:   jac->l_lens            = 0;
    686:   jac->tp_mat            = 0;
    687:   jac->tp_pmat           = 0;

    689:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiSetUseTrueLocal_C",
    690:                     "PCBJacobiSetUseTrueLocal_BJacobi",
    691:                     PCBJacobiSetUseTrueLocal_BJacobi);
    692:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiGetSubKSP_C","PCBJacobiGetSubKSP_BJacobi",
    693:                     PCBJacobiGetSubKSP_BJacobi);
    694:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiSetTotalBlocks_C","PCBJacobiSetTotalBlocks_BJacobi",
    695:                     PCBJacobiSetTotalBlocks_BJacobi);
    696:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiGetTotalBlocks_C","PCBJacobiGetTotalBlocks_BJacobi",
    697:                     PCBJacobiGetTotalBlocks_BJacobi);
    698:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiSetLocalBlocks_C","PCBJacobiSetLocalBlocks_BJacobi",
    699:                     PCBJacobiSetLocalBlocks_BJacobi);
    700:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiGetLocalBlocks_C","PCBJacobiGetLocalBlocks_BJacobi",
    701:                     PCBJacobiGetLocalBlocks_BJacobi);

    703:   return(0);
    704: }

    707: /* --------------------------------------------------------------------------------------------*/
    708: /*
    709:         These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
    710: */
    713: PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
    714: {
    715:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    716:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
    717:   PetscErrorCode         ierr;

    720:   /*
    721:         If the on processor block had to be generated via a MatGetDiagonalBlock()
    722:      that creates a copy (for example MPIBDiag matrices do), this frees the space
    723:   */
    724:   if (jac->tp_mat) {
    725:     MatDestroy(jac->tp_mat);
    726:   }
    727:   if (jac->tp_pmat) {
    728:     MatDestroy(jac->tp_pmat);
    729:   }

    731:   KSPDestroy(jac->ksp[0]);
    732:   PetscFree(jac->ksp);
    733:   VecDestroy(bjac->x);
    734:   VecDestroy(bjac->y);
    735:   if (jac->l_lens) {PetscFree(jac->l_lens);}
    736:   if (jac->g_lens) {PetscFree(jac->g_lens);}
    737:   PetscFree(bjac);
    738:   PetscFree(jac);
    739:   return(0);
    740: }

    744: PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
    745: {
    747:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;

    750:   KSPSetUp(jac->ksp[0]);
    751:   return(0);
    752: }

    756: PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
    757: {
    758:   PetscErrorCode         ierr;
    759:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    760:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
    761:   PetscScalar            *x_array,*y_array;

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

    784: PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
    785: {
    786:   PetscErrorCode         ierr;
    787:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    788:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
    789:   PetscScalar            *x_array,*y_array;
    790:   PC                     subpc;

    793:   /*
    794:       The VecPlaceArray() is to avoid having to copy the
    795:     y vector into the bjac->x vector. The reason for
    796:     the bjac->x vector is that we need a sequential vector
    797:     for the sequential solve.
    798:   */
    799:   VecGetArray(x,&x_array);
    800:   VecGetArray(y,&y_array);
    801:   VecPlaceArray(bjac->x,x_array);
    802:   VecPlaceArray(bjac->y,y_array);

    804:   /* apply the symmetric left portion of the inner PC operator */
    805:   /* note this by-passes the inner KSP and its options completely */

    807:   KSPGetPC(jac->ksp[0],&subpc);
    808:   PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
    809:   VecResetArray(bjac->x);
    810:   VecResetArray(bjac->y);

    812:   VecRestoreArray(x,&x_array);
    813:   VecRestoreArray(y,&y_array);
    814:   return(0);
    815: }

    819: PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
    820: {
    821:   PetscErrorCode         ierr;
    822:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    823:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
    824:   PetscScalar            *x_array,*y_array;
    825:   PC                     subpc;

    828:   /*
    829:       The VecPlaceArray() is to avoid having to copy the
    830:     y vector into the bjac->x vector. The reason for
    831:     the bjac->x vector is that we need a sequential vector
    832:     for the sequential solve.
    833:   */
    834:   VecGetArray(x,&x_array);
    835:   VecGetArray(y,&y_array);
    836:   VecPlaceArray(bjac->x,x_array);
    837:   VecPlaceArray(bjac->y,y_array);

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

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

    845:   VecRestoreArray(x,&x_array);
    846:   VecRestoreArray(y,&y_array);
    847:   return(0);
    848: }

    852: PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
    853: {
    854:   PetscErrorCode         ierr;
    855:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    856:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
    857:   PetscScalar            *x_array,*y_array;

    860:   /*
    861:       The VecPlaceArray() is to avoid having to copy the
    862:     y vector into the bjac->x vector. The reason for
    863:     the bjac->x vector is that we need a sequential vector
    864:     for the sequential solve.
    865:   */
    866:   VecGetArray(x,&x_array);
    867:   VecGetArray(y,&y_array);
    868:   VecPlaceArray(bjac->x,x_array);
    869:   VecPlaceArray(bjac->y,y_array);
    870:   KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
    871:   VecResetArray(bjac->x);
    872:   VecResetArray(bjac->y);
    873:   VecRestoreArray(x,&x_array);
    874:   VecRestoreArray(y,&y_array);
    875:   return(0);
    876: }

    880: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
    881: {
    882:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    883:   PetscErrorCode         ierr;
    884:   PetscInt               m;
    885:   KSP                    ksp;
    886:   Vec                    x,y;
    887:   PC_BJacobi_Singleblock *bjac;
    888:   PC                     subpc;
    889:   PetscTruth             wasSetup;

    893:   /* set default direct solver with no Krylov method */
    894:   if (!pc->setupcalled) {
    895:     const char *prefix;
    896:     wasSetup = PETSC_FALSE;
    897:     KSPCreate(PETSC_COMM_SELF,&ksp);
    898:     PetscLogObjectParent(pc,ksp);
    899:     KSPSetType(ksp,KSPPREONLY);
    900:     KSPGetPC(ksp,&subpc);
    901:     PCGetOptionsPrefix(pc,&prefix);
    902:     KSPSetOptionsPrefix(ksp,prefix);
    903:     KSPAppendOptionsPrefix(ksp,"sub_");
    904:     /*
    905:       The reason we need to generate these vectors is to serve
    906:       as the right-hand side and solution vector for the solve on the
    907:       block. We do not need to allocate space for the vectors since
    908:       that is provided via VecPlaceArray() just before the call to
    909:       KSPSolve() on the block.
    910:     */
    911:     MatGetSize(pmat,&m,&m);
    912:     VecCreateSeqWithArray(PETSC_COMM_SELF,m,PETSC_NULL,&x);
    913:     VecCreateSeqWithArray(PETSC_COMM_SELF,m,PETSC_NULL,&y);
    914:     PetscLogObjectParent(pc,x);
    915:     PetscLogObjectParent(pc,y);

    917:     pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
    918:     pc->ops->apply               = PCApply_BJacobi_Singleblock;
    919:     pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
    920:     pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
    921:     pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
    922:     pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

    924:     PetscMalloc(sizeof(PC_BJacobi_Singleblock),&bjac);
    925:     PetscLogObjectMemory(pc,sizeof(PC_BJacobi_Singleblock));
    926:     bjac->x      = x;
    927:     bjac->y      = y;

    929:     PetscMalloc(sizeof(KSP),&jac->ksp);
    930:     jac->ksp[0] = ksp;
    931:     jac->data    = (void*)bjac;
    932:   } else {
    933:     wasSetup = PETSC_TRUE;
    934:     ksp = jac->ksp[0];
    935:     bjac = (PC_BJacobi_Singleblock *)jac->data;
    936:   }
    937:   if (jac->use_true_local) {
    938:     KSPSetOperators(ksp,mat,pmat,pc->flag);
    939:   }  else {
    940:     KSPSetOperators(ksp,pmat,pmat,pc->flag);
    941:   }
    942:   if (!wasSetup) {
    943:     KSPSetFromOptions(ksp);
    944:   }
    945:   return(0);
    946: }

    948: /* ---------------------------------------------------------------------------------------------*/

    952: PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
    953: {
    954:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
    955:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
    956:   PetscErrorCode        ierr;
    957:   PetscInt              i;

    960:   MatDestroyMatrices(jac->n_local,&bjac->pmat);
    961:   if (jac->use_true_local) {
    962:     MatDestroyMatrices(jac->n_local,&bjac->mat);
    963:   }

    965:   /*
    966:         If the on processor block had to be generated via a MatGetDiagonalBlock()
    967:      that creates a copy (for example MPIBDiag matrices do), this frees the space
    968:   */
    969:   if (jac->tp_mat) {
    970:     MatDestroy(jac->tp_mat);
    971:   }
    972:   if (jac->tp_pmat) {
    973:     MatDestroy(jac->tp_pmat);
    974:   }

    976:   for (i=0; i<jac->n_local; i++) {
    977:     KSPDestroy(jac->ksp[i]);
    978:     VecDestroy(bjac->x[i]);
    979:     VecDestroy(bjac->y[i]);
    980:     ISDestroy(bjac->is[i]);
    981:   }
    982:   PetscFree(jac->ksp);
    983:   PetscFree(bjac->x);
    984:   PetscFree(bjac->starts);
    985:   PetscFree(bjac->is);
    986:   PetscFree(bjac);
    987:   if (jac->l_lens) {PetscFree(jac->l_lens);}
    988:   if (jac->g_lens) {PetscFree(jac->g_lens);}
    989:   PetscFree(jac);
    990:   return(0);
    991: }

    995: PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
    996: {
    997:   PC_BJacobi     *jac = (PC_BJacobi*)pc->data;
    999:   PetscInt       i,n_local = jac->n_local;

    1002:   for (i=0; i<n_local; i++) {
    1003:     KSPSetUp(jac->ksp[i]);
    1004:   }
    1005:   return(0);
    1006: }

    1008: /*
    1009:       Preconditioner for block Jacobi
    1010: */
    1013: PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
    1014: {
    1015:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
    1016:   PetscErrorCode        ierr;
    1017:   PetscInt              i,n_local = jac->n_local;
    1018:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
    1019:   PetscScalar           *xin,*yin;
    1020:   static PetscTruth     flag = PETSC_TRUE;
    1021:   static PetscEvent     SUBKspSolve;

    1024:   if (flag) {
    1025:     PetscLogEventRegister(&SUBKspSolve,"SubKspSolve",KSP_COOKIE);
    1026:     flag = PETSC_FALSE;
    1027:   }
    1028:   VecGetArray(x,&xin);
    1029:   VecGetArray(y,&yin);
    1030:   for (i=0; i<n_local; i++) {
    1031:     /*
    1032:        To avoid copying the subvector from x into a workspace we instead
    1033:        make the workspace vector array point to the subpart of the array of
    1034:        the global vector.
    1035:     */
    1036:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
    1037:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

    1039:     PetscLogEventBegin(SUBKspSolve,jac->ksp[i],bjac->x[i],bjac->y[i],0);
    1040:     KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
    1041:     PetscLogEventEnd(SUBKspSolve,jac->ksp[i],bjac->x[i],bjac->y[i],0);

    1043:     VecResetArray(bjac->x[i]);
    1044:     VecResetArray(bjac->y[i]);
    1045:   }
    1046:   VecRestoreArray(x,&xin);
    1047:   VecRestoreArray(y,&yin);
    1048:   return(0);
    1049: }

    1051: /*
    1052:       Preconditioner for block Jacobi
    1053: */
    1056: PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
    1057: {
    1058:   PC_BJacobi            *jac = (PC_BJacobi*)pc->data;
    1059:   PetscErrorCode        ierr;
    1060:   PetscInt              i,n_local = jac->n_local;
    1061:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
    1062:   PetscScalar           *xin,*yin;
    1063:   static PetscTruth     flag = PETSC_TRUE;
    1064:   static PetscEvent     SUBKspSolve;

    1067:   if (flag) {
    1068:     PetscLogEventRegister(&SUBKspSolve,"SubKspSolveTranspose",KSP_COOKIE);
    1069:     flag = PETSC_FALSE;
    1070:   }
    1071:   VecGetArray(x,&xin);
    1072:   VecGetArray(y,&yin);
    1073:   for (i=0; i<n_local; i++) {
    1074:     /*
    1075:        To avoid copying the subvector from x into a workspace we instead
    1076:        make the workspace vector array point to the subpart of the array of
    1077:        the global vector.
    1078:     */
    1079:     VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
    1080:     VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);

    1082:     PetscLogEventBegin(SUBKspSolve,jac->ksp[i],bjac->x[i],bjac->y[i],0);
    1083:     KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
    1084:     PetscLogEventEnd(SUBKspSolve,jac->ksp[i],bjac->x[i],bjac->y[i],0);
    1085:   }
    1086:   VecRestoreArray(x,&xin);
    1087:   VecRestoreArray(y,&yin);
    1088:   return(0);
    1089: }

    1093: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
    1094: {
    1095:   PC_BJacobi             *jac = (PC_BJacobi*)pc->data;
    1096:   PetscErrorCode         ierr;
    1097:   PetscInt               m,n_local,N,M,start,i;
    1098:   const char             *prefix,*pprefix,*mprefix;
    1099:   KSP                    ksp;
    1100:   Vec                    x,y;
    1101:   PC_BJacobi_Multiblock  *bjac = (PC_BJacobi_Multiblock*)jac->data;
    1102:   PC                     subpc;
    1103:   IS                     is;
    1104:   MatReuse               scall = MAT_REUSE_MATRIX;

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

    1109:   n_local = jac->n_local;

    1111:   if (jac->use_true_local) {
    1112:     if (mat->type != pmat->type) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
    1113:   }

    1115:   /* set default direct solver with no Krylov method */
    1116:   if (!pc->setupcalled) {
    1117:     scall                  = MAT_INITIAL_MATRIX;
    1118:     pc->ops->destroy       = PCDestroy_BJacobi_Multiblock;
    1119:     pc->ops->apply         = PCApply_BJacobi_Multiblock;
    1120:     pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
    1121:     pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;

    1123:     PetscMalloc(sizeof(PC_BJacobi_Multiblock),&bjac);
    1124:     PetscLogObjectMemory(pc,sizeof(PC_BJacobi_Multiblock));
    1125:     PetscMalloc(n_local*sizeof(KSP),&jac->ksp);
    1126:     PetscLogObjectMemory(pc,sizeof(n_local*sizeof(KSP)));
    1127:     PetscMalloc(2*n_local*sizeof(Vec),&bjac->x);
    1128:     PetscLogObjectMemory(pc,sizeof(2*n_local*sizeof(Vec)));
    1129:     bjac->y      = bjac->x + n_local;
    1130:     PetscMalloc(n_local*sizeof(PetscScalar),&bjac->starts);
    1131:     PetscLogObjectMemory(pc,sizeof(n_local*sizeof(PetscScalar)));
    1132:
    1133:     jac->data    = (void*)bjac;
    1134:     PetscMalloc(n_local*sizeof(IS),&bjac->is);
    1135:     PetscLogObjectMemory(pc,sizeof(n_local*sizeof(IS)));

    1137:     start = 0;
    1138:     for (i=0; i<n_local; i++) {
    1139:       KSPCreate(PETSC_COMM_SELF,&ksp);
    1140:       PetscLogObjectParent(pc,ksp);
    1141:       KSPSetType(ksp,KSPPREONLY);
    1142:       KSPGetPC(ksp,&subpc);
    1143:       PCGetOptionsPrefix(pc,&prefix);
    1144:       KSPSetOptionsPrefix(ksp,prefix);
    1145:       KSPAppendOptionsPrefix(ksp,"sub_");

    1147:       m = jac->l_lens[i];

    1149:       /*
    1150:       The reason we need to generate these vectors is to serve
    1151:       as the right-hand side and solution vector for the solve on the
    1152:       block. We do not need to allocate space for the vectors since
    1153:       that is provided via VecPlaceArray() just before the call to
    1154:       KSPSolve() on the block.

    1156:       */
    1157:       VecCreateSeq(PETSC_COMM_SELF,m,&x);
    1158:       VecCreateSeqWithArray(PETSC_COMM_SELF,m,PETSC_NULL,&y);
    1159:       PetscLogObjectParent(pc,x);
    1160:       PetscLogObjectParent(pc,y);
    1161:       bjac->x[i]      = x;
    1162:       bjac->y[i]      = y;
    1163:       bjac->starts[i] = start;
    1164:       jac->ksp[i]    = ksp;

    1166:       ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
    1167:       bjac->is[i] = is;
    1168:       PetscLogObjectParent(pc,is);

    1170:       start += m;
    1171:     }
    1172:   } else {
    1173:     bjac = (PC_BJacobi_Multiblock*)jac->data;
    1174:     /*
    1175:        Destroy the blocks from the previous iteration
    1176:     */
    1177:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    1178:       MatDestroyMatrices(n_local,&bjac->pmat);
    1179:       if (jac->use_true_local) {
    1180:         MatDestroyMatrices(n_local,&bjac->mat);
    1181:       }
    1182:       scall = MAT_INITIAL_MATRIX;
    1183:     }
    1184:   }

    1186:   MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
    1187:   if (jac->use_true_local) {
    1188:     PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
    1189:     MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
    1190:   }
    1191:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
    1192:      different boundary conditions for the submatrices than for the global problem) */
    1193:   PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);

    1195:   PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
    1196:   for (i=0; i<n_local; i++) {
    1197:     PetscLogObjectParent(pc,bjac->pmat[i]);
    1198:     PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
    1199:     if (jac->use_true_local) {
    1200:       PetscLogObjectParent(pc,bjac->mat[i]);
    1201:       PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
    1202:       KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i],pc->flag);
    1203:     } else {
    1204:       KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i],pc->flag);
    1205:     }
    1206:     KSPSetFromOptions(jac->ksp[i]);
    1207:   }

    1209:   return(0);
    1210: }