Actual source code: mpisell.c

      1: #include <../src/mat/impls/aij/mpi/mpiaij.h>
      2: #include <../src/mat/impls/sell/mpi/mpisell.h>
      3: #include <petsc/private/vecimpl.h>
      4: #include <petsc/private/isimpl.h>
      5: #include <petscblaslapack.h>
      6: #include <petscsf.h>

      8: /*MC
      9:    MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices.

     11:    This matrix type is identical to MATSEQSELL when constructed with a single process communicator,
     12:    and MATMPISELL otherwise.  As a result, for single process communicators,
     13:   MatSeqSELLSetPreallocation is supported, and similarly MatMPISELLSetPreallocation is supported
     14:   for communicators controlling multiple processes.  It is recommended that you call both of
     15:   the above preallocation routines for simplicity.

     17:    Options Database Keys:
     18: . -mat_type sell - sets the matrix type to "sell" during a call to MatSetFromOptions()

     20:   Level: beginner

     22: .seealso: MatCreateSELL(), MatCreateSeqSELL(), MATSEQSELL, MATMPISELL
     23: M*/

     25: PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is)
     26: {
     28:   Mat_MPISELL    *sell=(Mat_MPISELL*)Y->data;

     31:   if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) {
     32:     MatDiagonalSet(sell->A,D,is);
     33:   } else {
     34:     MatDiagonalSet_Default(Y,D,is);
     35:   }
     36:   return(0);
     37: }

     39: /*
     40:   Local utility routine that creates a mapping from the global column
     41: number to the local number in the off-diagonal part of the local
     42: storage of the matrix.  When PETSC_USE_CTABLE is used this is scalable at
     43: a slightly higher hash table cost; without it it is not scalable (each processor
     44: has an order N integer array but is fast to acess.
     45: */
     46: PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat)
     47: {
     48:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
     50:   PetscInt       n=sell->B->cmap->n,i;

     53:   if (!sell->garray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MPISELL Matrix was assembled but is missing garray");
     54: #if defined(PETSC_USE_CTABLE)
     55:   PetscTableCreate(n,mat->cmap->N+1,&sell->colmap);
     56:   for (i=0; i<n; i++) {
     57:     PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES);
     58:   }
     59: #else
     60:   PetscCalloc1(mat->cmap->N+1,&sell->colmap);
     61:   PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt));
     62:   for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1;
     63: #endif
     64:   return(0);
     65: }

     67: #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \
     68:   { \
     69:     if (col <= lastcol1) low1 = 0; \
     70:     else                high1 = nrow1; \
     71:     lastcol1 = col; \
     72:     while (high1-low1 > 5) { \
     73:       t = (low1+high1)/2; \
     74:       if (*(cp1+8*t) > col) high1 = t; \
     75:       else                   low1 = t; \
     76:     } \
     77:     for (_i=low1; _i<high1; _i++) { \
     78:       if (*(cp1+8*_i) > col) break; \
     79:       if (*(cp1+8*_i) == col) { \
     80:         if (addv == ADD_VALUES) *(vp1+8*_i) += value;   \
     81:         else                     *(vp1+8*_i) = value; \
     82:         goto a_noinsert; \
     83:       } \
     84:     }  \
     85:     if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \
     86:     if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \
     87:     if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \
     88:     MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \
     89:     /* shift up all the later entries in this row */ \
     90:     for (ii=nrow1-1; ii>=_i; ii--) { \
     91:       *(cp1+8*(ii+1)) = *(cp1+8*ii); \
     92:       *(vp1+8*(ii+1)) = *(vp1+8*ii); \
     93:     } \
     94:     *(cp1+8*_i) = col; \
     95:     *(vp1+8*_i) = value; \
     96:     a->nz++; nrow1++; A->nonzerostate++; \
     97:     a_noinsert: ; \
     98:     a->rlen[row] = nrow1; \
     99:   }

    101: #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \
    102:   { \
    103:     if (col <= lastcol2) low2 = 0; \
    104:     else                high2 = nrow2; \
    105:     lastcol2 = col; \
    106:     while (high2-low2 > 5) { \
    107:       t = (low2+high2)/2; \
    108:       if (*(cp2+8*t) > col) high2 = t; \
    109:       else low2  = t; \
    110:     } \
    111:     for (_i=low2; _i<high2; _i++) { \
    112:       if (*(cp2+8*_i) > col) break; \
    113:       if (*(cp2+8*_i) == col) { \
    114:         if (addv == ADD_VALUES) *(vp2+8*_i) += value; \
    115:         else                     *(vp2+8*_i) = value; \
    116:         goto b_noinsert; \
    117:       } \
    118:     } \
    119:     if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
    120:     if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
    121:     if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \
    122:     MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \
    123:     /* shift up all the later entries in this row */ \
    124:     for (ii=nrow2-1; ii>=_i; ii--) { \
    125:       *(cp2+8*(ii+1)) = *(cp2+8*ii); \
    126:       *(vp2+8*(ii+1)) = *(vp2+8*ii); \
    127:     } \
    128:     *(cp2+8*_i) = col; \
    129:     *(vp2+8*_i) = value; \
    130:     b->nz++; nrow2++; B->nonzerostate++; \
    131:     b_noinsert: ; \
    132:     b->rlen[row] = nrow2; \
    133:   }

    135: PetscErrorCode MatSetValues_MPISELL(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
    136: {
    137:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
    138:   PetscScalar    value;
    140:   PetscInt       i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend,shift1,shift2;
    141:   PetscInt       cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;
    142:   PetscBool      roworiented=sell->roworiented;

    144:   /* Some Variables required in the macro */
    145:   Mat            A=sell->A;
    146:   Mat_SeqSELL    *a=(Mat_SeqSELL*)A->data;
    147:   PetscBool      ignorezeroentries=a->ignorezeroentries,found;
    148:   Mat            B=sell->B;
    149:   Mat_SeqSELL    *b=(Mat_SeqSELL*)B->data;
    150:   PetscInt       *cp1,*cp2,ii,_i,nrow1,nrow2,low1,high1,low2,high2,t,lastcol1,lastcol2;
    151:   MatScalar      *vp1,*vp2;

    154:   for (i=0; i<m; i++) {
    155:     if (im[i] < 0) continue;
    156:     if (PetscUnlikelyDebug(im[i] >= mat->rmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap->N-1);
    157:     if (im[i] >= rstart && im[i] < rend) {
    158:       row      = im[i] - rstart;
    159:       lastcol1 = -1;
    160:       shift1   = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */
    161:       cp1      = a->colidx+shift1;
    162:       vp1      = a->val+shift1;
    163:       nrow1    = a->rlen[row];
    164:       low1     = 0;
    165:       high1    = nrow1;
    166:       lastcol2 = -1;
    167:       shift2   = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
    168:       cp2      = b->colidx+shift2;
    169:       vp2      = b->val+shift2;
    170:       nrow2    = b->rlen[row];
    171:       low2     = 0;
    172:       high2    = nrow2;

    174:       for (j=0; j<n; j++) {
    175:         if (roworiented) value = v[i*n+j];
    176:         else             value = v[i+j*m];
    177:         if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue;
    178:         if (in[j] >= cstart && in[j] < cend) {
    179:           col   = in[j] - cstart;
    180:           MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */
    181:         } else if (in[j] < 0) continue;
    182:         else if (PetscUnlikelyDebug(in[j] >= mat->cmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[j],mat->cmap->N-1);
    183:         else {
    184:           if (mat->was_assembled) {
    185:             if (!sell->colmap) {
    186:               MatCreateColmap_MPISELL_Private(mat);
    187:             }
    188: #if defined(PETSC_USE_CTABLE)
    189:             PetscTableFind(sell->colmap,in[j]+1,&col);
    190:             col--;
    191: #else
    192:             col = sell->colmap[in[j]] - 1;
    193: #endif
    194:             if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) {
    195:               MatDisAssemble_MPISELL(mat);
    196:               col    = in[j];
    197:               /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */
    198:               B      = sell->B;
    199:               b      = (Mat_SeqSELL*)B->data;
    200:               shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
    201:               cp2    = b->colidx+shift2;
    202:               vp2    = b->val+shift2;
    203:               nrow2  = b->rlen[row];
    204:               low2   = 0;
    205:               high2  = nrow2;
    206:             } else if (col < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", im[i], in[j]);
    207:           } else col = in[j];
    208:           MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */
    209:         }
    210:       }
    211:     } else {
    212:       if (mat->nooffprocentries) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process row %D even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set",im[i]);
    213:       if (!sell->donotstash) {
    214:         mat->assembled = PETSC_FALSE;
    215:         if (roworiented) {
    216:           MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
    217:         } else {
    218:           MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
    219:         }
    220:       }
    221:     }
    222:   }
    223:   return(0);
    224: }

    226: PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
    227: {
    228:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
    230:   PetscInt       i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend;
    231:   PetscInt       cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;

    234:   for (i=0; i<m; i++) {
    235:     if (idxm[i] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",idxm[i]);*/
    236:     if (idxm[i] >= mat->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm[i],mat->rmap->N-1);
    237:     if (idxm[i] >= rstart && idxm[i] < rend) {
    238:       row = idxm[i] - rstart;
    239:       for (j=0; j<n; j++) {
    240:         if (idxn[j] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",idxn[j]); */
    241:         if (idxn[j] >= mat->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",idxn[j],mat->cmap->N-1);
    242:         if (idxn[j] >= cstart && idxn[j] < cend) {
    243:           col  = idxn[j] - cstart;
    244:           MatGetValues(sell->A,1,&row,1,&col,v+i*n+j);
    245:         } else {
    246:           if (!sell->colmap) {
    247:             MatCreateColmap_MPISELL_Private(mat);
    248:           }
    249: #if defined(PETSC_USE_CTABLE)
    250:           PetscTableFind(sell->colmap,idxn[j]+1,&col);
    251:           col--;
    252: #else
    253:           col = sell->colmap[idxn[j]] - 1;
    254: #endif
    255:           if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0;
    256:           else {
    257:             MatGetValues(sell->B,1,&row,1,&col,v+i*n+j);
    258:           }
    259:         }
    260:       }
    261:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported");
    262:   }
    263:   return(0);
    264: }

    266: extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec);

    268: PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode)
    269: {
    270:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
    272:   PetscInt       nstash,reallocs;

    275:   if (sell->donotstash || mat->nooffprocentries) return(0);

    277:   MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
    278:   MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);
    279:   PetscInfo2(sell->A,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);
    280:   return(0);
    281: }

    283: PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode)
    284: {
    285:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
    287:   PetscMPIInt    n;
    288:   PetscInt       i,flg;
    289:   PetscInt       *row,*col;
    290:   PetscScalar    *val;
    291:   PetscBool      other_disassembled;
    292:   /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */
    294:   if (!sell->donotstash && !mat->nooffprocentries) {
    295:     while (1) {
    296:       MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);
    297:       if (!flg) break;

    299:       for (i=0; i<n; i++) { /* assemble one by one */
    300:         MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode);
    301:       }
    302:     }
    303:     MatStashScatterEnd_Private(&mat->stash);
    304:   }
    305:   MatAssemblyBegin(sell->A,mode);
    306:   MatAssemblyEnd(sell->A,mode);

    308:   /*
    309:      determine if any processor has disassembled, if so we must
    310:      also disassemble ourselfs, in order that we may reassemble.
    311:   */
    312:   /*
    313:      if nonzero structure of submatrix B cannot change then we know that
    314:      no processor disassembled thus we can skip this stuff
    315:   */
    316:   if (!((Mat_SeqSELL*)sell->B->data)->nonew) {
    317:     MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat));
    318:     if (mat->was_assembled && !other_disassembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet\n");
    319:   }
    320:   if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) {
    321:     MatSetUpMultiply_MPISELL(mat);
    322:   }
    323:   /*
    324:   MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE);
    325:   */
    326:   MatAssemblyBegin(sell->B,mode);
    327:   MatAssemblyEnd(sell->B,mode);
    328:   PetscFree2(sell->rowvalues,sell->rowindices);
    329:   sell->rowvalues = NULL;
    330:   VecDestroy(&sell->diag);

    332:   /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
    333:   if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) {
    334:     PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate;
    335:     MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));
    336:   }
    337:   return(0);
    338: }

    340: PetscErrorCode MatZeroEntries_MPISELL(Mat A)
    341: {
    342:   Mat_MPISELL    *l=(Mat_MPISELL*)A->data;

    346:   MatZeroEntries(l->A);
    347:   MatZeroEntries(l->B);
    348:   return(0);
    349: }

    351: PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy)
    352: {
    353:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
    355:   PetscInt       nt;

    358:   VecGetLocalSize(xx,&nt);
    359:   if (nt != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%D) and xx (%D)",A->cmap->n,nt);
    360:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
    361:   (*a->A->ops->mult)(a->A,xx,yy);
    362:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
    363:   (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);
    364:   return(0);
    365: }

    367: PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx)
    368: {
    369:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    373:   MatMultDiagonalBlock(a->A,bb,xx);
    374:   return(0);
    375: }

    377: PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
    378: {
    379:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    383:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
    384:   (*a->A->ops->multadd)(a->A,xx,yy,zz);
    385:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
    386:   (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
    387:   return(0);
    388: }

    390: PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy)
    391: {
    392:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    396:   /* do nondiagonal part */
    397:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
    398:   /* do local part */
    399:   (*a->A->ops->multtranspose)(a->A,xx,yy);
    400:   /* add partial results together */
    401:   VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
    402:   VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
    403:   return(0);
    404: }

    406: PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f)
    407: {
    408:   MPI_Comm       comm;
    409:   Mat_MPISELL    *Asell=(Mat_MPISELL*)Amat->data,*Bsell;
    410:   Mat            Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs;
    411:   IS             Me,Notme;
    413:   PetscInt       M,N,first,last,*notme,i;
    414:   PetscMPIInt    size;

    417:   /* Easy test: symmetric diagonal block */
    418:   Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A;
    419:   MatIsTranspose(Adia,Bdia,tol,f);
    420:   if (!*f) return(0);
    421:   PetscObjectGetComm((PetscObject)Amat,&comm);
    422:   MPI_Comm_size(comm,&size);
    423:   if (size == 1) return(0);

    425:   /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */
    426:   MatGetSize(Amat,&M,&N);
    427:   MatGetOwnershipRange(Amat,&first,&last);
    428:   PetscMalloc1(N-last+first,&notme);
    429:   for (i=0; i<first; i++) notme[i] = i;
    430:   for (i=last; i<M; i++) notme[i-last+first] = i;
    431:   ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme);
    432:   ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me);
    433:   MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs);
    434:   Aoff = Aoffs[0];
    435:   MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs);
    436:   Boff = Boffs[0];
    437:   MatIsTranspose(Aoff,Boff,tol,f);
    438:   MatDestroyMatrices(1,&Aoffs);
    439:   MatDestroyMatrices(1,&Boffs);
    440:   ISDestroy(&Me);
    441:   ISDestroy(&Notme);
    442:   PetscFree(notme);
    443:   return(0);
    444: }

    446: PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
    447: {
    448:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    452:   /* do nondiagonal part */
    453:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
    454:   /* do local part */
    455:   (*a->A->ops->multtransposeadd)(a->A,xx,yy,zz);
    456:   /* add partial results together */
    457:   VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
    458:   VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
    459:   return(0);
    460: }

    462: /*
    463:   This only works correctly for square matrices where the subblock A->A is the
    464:    diagonal block
    465: */
    466: PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v)
    467: {
    469:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    472:   if (A->rmap->N != A->cmap->N) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block");
    473:   if (A->rmap->rstart != A->cmap->rstart || A->rmap->rend != A->cmap->rend) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"row partition must equal col partition");
    474:   MatGetDiagonal(a->A,v);
    475:   return(0);
    476: }

    478: PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa)
    479: {
    480:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    484:   MatScale(a->A,aa);
    485:   MatScale(a->B,aa);
    486:   return(0);
    487: }

    489: PetscErrorCode MatDestroy_MPISELL(Mat mat)
    490: {
    491:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

    495: #if defined(PETSC_USE_LOG)
    496:   PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D",mat->rmap->N,mat->cmap->N);
    497: #endif
    498:   MatStashDestroy_Private(&mat->stash);
    499:   VecDestroy(&sell->diag);
    500:   MatDestroy(&sell->A);
    501:   MatDestroy(&sell->B);
    502: #if defined(PETSC_USE_CTABLE)
    503:   PetscTableDestroy(&sell->colmap);
    504: #else
    505:   PetscFree(sell->colmap);
    506: #endif
    507:   PetscFree(sell->garray);
    508:   VecDestroy(&sell->lvec);
    509:   VecScatterDestroy(&sell->Mvctx);
    510:   PetscFree2(sell->rowvalues,sell->rowindices);
    511:   PetscFree(sell->ld);
    512:   PetscFree(mat->data);

    514:   PetscObjectChangeTypeName((PetscObject)mat,NULL);
    515:   PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);
    516:   PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);
    517:   PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL);
    518:   PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL);
    519:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL);
    520:   PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);
    521:   return(0);
    522: }

    524: #include <petscdraw.h>
    525: PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer)
    526: {
    527:   Mat_MPISELL       *sell=(Mat_MPISELL*)mat->data;
    528:   PetscErrorCode    ierr;
    529:   PetscMPIInt       rank=sell->rank,size=sell->size;
    530:   PetscBool         isdraw,iascii,isbinary;
    531:   PetscViewer       sviewer;
    532:   PetscViewerFormat format;

    535:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
    536:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    537:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
    538:   if (iascii) {
    539:     PetscViewerGetFormat(viewer,&format);
    540:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    541:       MatInfo   info;
    542:       PetscInt *inodes;

    544:       MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank);
    545:       MatGetInfo(mat,MAT_LOCAL,&info);
    546:       MatInodeGetInodeSizes(sell->A,NULL,&inodes,NULL);
    547:       PetscViewerASCIIPushSynchronized(viewer);
    548:       if (!inodes) {
    549:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D mem %D, not using I-node routines\n",
    550:                                                   rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
    551:       } else {
    552:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D mem %D, using I-node routines\n",
    553:                                                   rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
    554:       }
    555:       MatGetInfo(sell->A,MAT_LOCAL,&info);
    556:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);
    557:       MatGetInfo(sell->B,MAT_LOCAL,&info);
    558:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);
    559:       PetscViewerFlush(viewer);
    560:       PetscViewerASCIIPopSynchronized(viewer);
    561:       PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n");
    562:       VecScatterView(sell->Mvctx,viewer);
    563:       return(0);
    564:     } else if (format == PETSC_VIEWER_ASCII_INFO) {
    565:       PetscInt inodecount,inodelimit,*inodes;
    566:       MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit);
    567:       if (inodes) {
    568:         PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %D nodes, limit used is %D\n",inodecount,inodelimit);
    569:       } else {
    570:         PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n");
    571:       }
    572:       return(0);
    573:     } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
    574:       return(0);
    575:     }
    576:   } else if (isbinary) {
    577:     if (size == 1) {
    578:       PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name);
    579:       MatView(sell->A,viewer);
    580:     } else {
    581:       /* MatView_MPISELL_Binary(mat,viewer); */
    582:     }
    583:     return(0);
    584:   } else if (isdraw) {
    585:     PetscDraw draw;
    586:     PetscBool isnull;
    587:     PetscViewerDrawGetDraw(viewer,0,&draw);
    588:     PetscDrawIsNull(draw,&isnull);
    589:     if (isnull) return(0);
    590:   }

    592:   {
    593:     /* assemble the entire matrix onto first processor. */
    594:     Mat         A;
    595:     Mat_SeqSELL *Aloc;
    596:     PetscInt    M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j;
    597:     MatScalar   *aval;
    598:     PetscBool   isnonzero;

    600:     MatCreate(PetscObjectComm((PetscObject)mat),&A);
    601:     if (!rank) {
    602:       MatSetSizes(A,M,N,M,N);
    603:     } else {
    604:       MatSetSizes(A,0,0,M,N);
    605:     }
    606:     /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */
    607:     MatSetType(A,MATMPISELL);
    608:     MatMPISELLSetPreallocation(A,0,NULL,0,NULL);
    609:     MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
    610:     PetscLogObjectParent((PetscObject)mat,(PetscObject)A);

    612:     /* copy over the A part */
    613:     Aloc = (Mat_SeqSELL*)sell->A->data;
    614:     acolidx = Aloc->colidx; aval = Aloc->val;
    615:     for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */
    616:       for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
    617:         isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
    618:         if (isnonzero) { /* check the mask bit */
    619:           row  = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */
    620:           col  = *acolidx + mat->rmap->rstart;
    621:           MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
    622:         }
    623:         aval++; acolidx++;
    624:       }
    625:     }

    627:     /* copy over the B part */
    628:     Aloc = (Mat_SeqSELL*)sell->B->data;
    629:     acolidx = Aloc->colidx; aval = Aloc->val;
    630:     for (i=0; i<Aloc->totalslices; i++) {
    631:       for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
    632:         isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
    633:         if (isnonzero) {
    634:           row  = (i<<3)+(j&0x07) + mat->rmap->rstart;
    635:           col  = sell->garray[*acolidx];
    636:           MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
    637:         }
    638:         aval++; acolidx++;
    639:       }
    640:     }

    642:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    643:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    644:     /*
    645:        Everyone has to call to draw the matrix since the graphics waits are
    646:        synchronized across all processors that share the PetscDraw object
    647:     */
    648:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    649:     if (!rank) {
    650:       PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name);
    651:       MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer);
    652:     }
    653:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    654:     PetscViewerFlush(viewer);
    655:     MatDestroy(&A);
    656:   }
    657:   return(0);
    658: }

    660: PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer)
    661: {
    663:   PetscBool      iascii,isdraw,issocket,isbinary;

    666:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    667:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
    668:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
    669:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
    670:   if (iascii || isdraw || isbinary || issocket) {
    671:     MatView_MPISELL_ASCIIorDraworSocket(mat,viewer);
    672:   }
    673:   return(0);
    674: }

    676: PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
    677: {
    678:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

    682:   MatGetSize(sell->B,NULL,nghosts);
    683:   if (ghosts) *ghosts = sell->garray;
    684:   return(0);
    685: }

    687: PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info)
    688: {
    689:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
    690:   Mat            A=mat->A,B=mat->B;
    692:   PetscLogDouble isend[5],irecv[5];

    695:   info->block_size = 1.0;
    696:   MatGetInfo(A,MAT_LOCAL,info);

    698:   isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded;
    699:   isend[3] = info->memory;  isend[4] = info->mallocs;

    701:   MatGetInfo(B,MAT_LOCAL,info);

    703:   isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded;
    704:   isend[3] += info->memory;  isend[4] += info->mallocs;
    705:   if (flag == MAT_LOCAL) {
    706:     info->nz_used      = isend[0];
    707:     info->nz_allocated = isend[1];
    708:     info->nz_unneeded  = isend[2];
    709:     info->memory       = isend[3];
    710:     info->mallocs      = isend[4];
    711:   } else if (flag == MAT_GLOBAL_MAX) {
    712:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin));

    714:     info->nz_used      = irecv[0];
    715:     info->nz_allocated = irecv[1];
    716:     info->nz_unneeded  = irecv[2];
    717:     info->memory       = irecv[3];
    718:     info->mallocs      = irecv[4];
    719:   } else if (flag == MAT_GLOBAL_SUM) {
    720:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin));

    722:     info->nz_used      = irecv[0];
    723:     info->nz_allocated = irecv[1];
    724:     info->nz_unneeded  = irecv[2];
    725:     info->memory       = irecv[3];
    726:     info->mallocs      = irecv[4];
    727:   }
    728:   info->fill_ratio_given  = 0; /* no parallel LU/ILU/Cholesky */
    729:   info->fill_ratio_needed = 0;
    730:   info->factor_mallocs    = 0;
    731:   return(0);
    732: }

    734: PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg)
    735: {
    736:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    740:   switch (op) {
    741:   case MAT_NEW_NONZERO_LOCATIONS:
    742:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
    743:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
    744:   case MAT_KEEP_NONZERO_PATTERN:
    745:   case MAT_NEW_NONZERO_LOCATION_ERR:
    746:   case MAT_USE_INODES:
    747:   case MAT_IGNORE_ZERO_ENTRIES:
    748:     MatCheckPreallocated(A,1);
    749:     MatSetOption(a->A,op,flg);
    750:     MatSetOption(a->B,op,flg);
    751:     break;
    752:   case MAT_ROW_ORIENTED:
    753:     MatCheckPreallocated(A,1);
    754:     a->roworiented = flg;

    756:     MatSetOption(a->A,op,flg);
    757:     MatSetOption(a->B,op,flg);
    758:     break;
    759:   case MAT_FORCE_DIAGONAL_ENTRIES:
    760:   case MAT_SORTED_FULL:
    761:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
    762:     break;
    763:   case MAT_IGNORE_OFF_PROC_ENTRIES:
    764:     a->donotstash = flg;
    765:     break;
    766:   case MAT_SPD:
    767:     A->spd_set = PETSC_TRUE;
    768:     A->spd     = flg;
    769:     if (flg) {
    770:       A->symmetric                  = PETSC_TRUE;
    771:       A->structurally_symmetric     = PETSC_TRUE;
    772:       A->symmetric_set              = PETSC_TRUE;
    773:       A->structurally_symmetric_set = PETSC_TRUE;
    774:     }
    775:     break;
    776:   case MAT_SYMMETRIC:
    777:     MatCheckPreallocated(A,1);
    778:     MatSetOption(a->A,op,flg);
    779:     break;
    780:   case MAT_STRUCTURALLY_SYMMETRIC:
    781:     MatCheckPreallocated(A,1);
    782:     MatSetOption(a->A,op,flg);
    783:     break;
    784:   case MAT_HERMITIAN:
    785:     MatCheckPreallocated(A,1);
    786:     MatSetOption(a->A,op,flg);
    787:     break;
    788:   case MAT_SYMMETRY_ETERNAL:
    789:     MatCheckPreallocated(A,1);
    790:     MatSetOption(a->A,op,flg);
    791:     break;
    792:   default:
    793:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
    794:   }
    795:   return(0);
    796: }

    799: PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr)
    800: {
    801:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
    802:   Mat            a=sell->A,b=sell->B;
    804:   PetscInt       s1,s2,s3;

    807:   MatGetLocalSize(mat,&s2,&s3);
    808:   if (rr) {
    809:     VecGetLocalSize(rr,&s1);
    810:     if (s1!=s3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size");
    811:     /* Overlap communication with computation. */
    812:     VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
    813:   }
    814:   if (ll) {
    815:     VecGetLocalSize(ll,&s1);
    816:     if (s1!=s2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size");
    817:     (*b->ops->diagonalscale)(b,ll,NULL);
    818:   }
    819:   /* scale  the diagonal block */
    820:   (*a->ops->diagonalscale)(a,ll,rr);

    822:   if (rr) {
    823:     /* Do a scatter end and then right scale the off-diagonal block */
    824:     VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
    825:     (*b->ops->diagonalscale)(b,NULL,sell->lvec);
    826:   }
    827:   return(0);
    828: }

    830: PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
    831: {
    832:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    836:   MatSetUnfactored(a->A);
    837:   return(0);
    838: }

    840: PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool  *flag)
    841: {
    842:   Mat_MPISELL    *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data;
    843:   Mat            a,b,c,d;
    844:   PetscBool      flg;

    848:   a = matA->A; b = matA->B;
    849:   c = matB->A; d = matB->B;

    851:   MatEqual(a,c,&flg);
    852:   if (flg) {
    853:     MatEqual(b,d,&flg);
    854:   }
    855:   MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
    856:   return(0);
    857: }

    859: PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str)
    860: {
    862:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
    863:   Mat_MPISELL    *b=(Mat_MPISELL*)B->data;

    866:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
    867:   if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
    868:     /* because of the column compression in the off-processor part of the matrix a->B,
    869:        the number of columns in a->B and b->B may be different, hence we cannot call
    870:        the MatCopy() directly on the two parts. If need be, we can provide a more
    871:        efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices
    872:        then copying the submatrices */
    873:     MatCopy_Basic(A,B,str);
    874:   } else {
    875:     MatCopy(a->A,b->A,str);
    876:     MatCopy(a->B,b->B,str);
    877:   }
    878:   return(0);
    879: }

    881: PetscErrorCode MatSetUp_MPISELL(Mat A)
    882: {

    886:    MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
    887:   return(0);
    888: }

    891: extern PetscErrorCode MatConjugate_SeqSELL(Mat);

    893: PetscErrorCode MatConjugate_MPISELL(Mat mat)
    894: {
    895: #if defined(PETSC_USE_COMPLEX)
    897:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

    900:   MatConjugate_SeqSELL(sell->A);
    901:   MatConjugate_SeqSELL(sell->B);
    902: #else
    904: #endif
    905:   return(0);
    906: }

    908: PetscErrorCode MatRealPart_MPISELL(Mat A)
    909: {
    910:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    914:   MatRealPart(a->A);
    915:   MatRealPart(a->B);
    916:   return(0);
    917: }

    919: PetscErrorCode MatImaginaryPart_MPISELL(Mat A)
    920: {
    921:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    925:   MatImaginaryPart(a->A);
    926:   MatImaginaryPart(a->B);
    927:   return(0);
    928: }

    930: PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values)
    931: {
    932:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    936:   MatInvertBlockDiagonal(a->A,values);
    937:   A->factorerrortype = a->A->factorerrortype;
    938:   return(0);
    939: }

    941: static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx)
    942: {
    944:   Mat_MPISELL    *sell=(Mat_MPISELL*)x->data;

    947:   MatSetRandom(sell->A,rctx);
    948:   MatSetRandom(sell->B,rctx);
    949:   MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
    950:   MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
    951:   return(0);
    952: }

    954: PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A)
    955: {

    959:   PetscOptionsHead(PetscOptionsObject,"MPISELL options");
    960:   PetscOptionsTail();
    961:   return(0);
    962: }

    964: PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a)
    965: {
    967:   Mat_MPISELL    *msell=(Mat_MPISELL*)Y->data;
    968:   Mat_SeqSELL    *sell=(Mat_SeqSELL*)msell->A->data;

    971:   if (!Y->preallocated) {
    972:     MatMPISELLSetPreallocation(Y,1,NULL,0,NULL);
    973:   } else if (!sell->nz) {
    974:     PetscInt nonew = sell->nonew;
    975:     MatSeqSELLSetPreallocation(msell->A,1,NULL);
    976:     sell->nonew = nonew;
    977:   }
    978:   MatShift_Basic(Y,a);
    979:   return(0);
    980: }

    982: PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool  *missing,PetscInt *d)
    983: {
    984:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

    988:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices");
    989:   MatMissingDiagonal(a->A,missing,d);
    990:   if (d) {
    991:     PetscInt rstart;
    992:     MatGetOwnershipRange(A,&rstart,NULL);
    993:     *d += rstart;

    995:   }
    996:   return(0);
    997: }

    999: PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a)
    1000: {
    1002:   *a = ((Mat_MPISELL*)A->data)->A;
    1003:   return(0);
    1004: }

    1006: /* -------------------------------------------------------------------*/
    1007: static struct _MatOps MatOps_Values = {MatSetValues_MPISELL,
    1008:                                        NULL,
    1009:                                        NULL,
    1010:                                        MatMult_MPISELL,
    1011:                                 /* 4*/ MatMultAdd_MPISELL,
    1012:                                        MatMultTranspose_MPISELL,
    1013:                                        MatMultTransposeAdd_MPISELL,
    1014:                                        NULL,
    1015:                                        NULL,
    1016:                                        NULL,
    1017:                                 /*10*/ NULL,
    1018:                                        NULL,
    1019:                                        NULL,
    1020:                                        MatSOR_MPISELL,
    1021:                                        NULL,
    1022:                                 /*15*/ MatGetInfo_MPISELL,
    1023:                                        MatEqual_MPISELL,
    1024:                                        MatGetDiagonal_MPISELL,
    1025:                                        MatDiagonalScale_MPISELL,
    1026:                                        NULL,
    1027:                                 /*20*/ MatAssemblyBegin_MPISELL,
    1028:                                        MatAssemblyEnd_MPISELL,
    1029:                                        MatSetOption_MPISELL,
    1030:                                        MatZeroEntries_MPISELL,
    1031:                                 /*24*/ NULL,
    1032:                                        NULL,
    1033:                                        NULL,
    1034:                                        NULL,
    1035:                                        NULL,
    1036:                                 /*29*/ MatSetUp_MPISELL,
    1037:                                        NULL,
    1038:                                        NULL,
    1039:                                        MatGetDiagonalBlock_MPISELL,
    1040:                                        NULL,
    1041:                                 /*34*/ MatDuplicate_MPISELL,
    1042:                                        NULL,
    1043:                                        NULL,
    1044:                                        NULL,
    1045:                                        NULL,
    1046:                                 /*39*/ NULL,
    1047:                                        NULL,
    1048:                                        NULL,
    1049:                                        MatGetValues_MPISELL,
    1050:                                        MatCopy_MPISELL,
    1051:                                 /*44*/ NULL,
    1052:                                        MatScale_MPISELL,
    1053:                                        MatShift_MPISELL,
    1054:                                        MatDiagonalSet_MPISELL,
    1055:                                        NULL,
    1056:                                 /*49*/ MatSetRandom_MPISELL,
    1057:                                        NULL,
    1058:                                        NULL,
    1059:                                        NULL,
    1060:                                        NULL,
    1061:                                 /*54*/ MatFDColoringCreate_MPIXAIJ,
    1062:                                        NULL,
    1063:                                        MatSetUnfactored_MPISELL,
    1064:                                        NULL,
    1065:                                        NULL,
    1066:                                 /*59*/ NULL,
    1067:                                        MatDestroy_MPISELL,
    1068:                                        MatView_MPISELL,
    1069:                                        NULL,
    1070:                                        NULL,
    1071:                                 /*64*/ NULL,
    1072:                                        NULL,
    1073:                                        NULL,
    1074:                                        NULL,
    1075:                                        NULL,
    1076:                                 /*69*/ NULL,
    1077:                                        NULL,
    1078:                                        NULL,
    1079:                                        NULL,
    1080:                                        NULL,
    1081:                                        NULL,
    1082:                                 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */
    1083:                                        MatSetFromOptions_MPISELL,
    1084:                                        NULL,
    1085:                                        NULL,
    1086:                                        NULL,
    1087:                                 /*80*/ NULL,
    1088:                                        NULL,
    1089:                                        NULL,
    1090:                                 /*83*/ NULL,
    1091:                                        NULL,
    1092:                                        NULL,
    1093:                                        NULL,
    1094:                                        NULL,
    1095:                                        NULL,
    1096:                                 /*89*/ NULL,
    1097:                                        NULL,
    1098:                                        NULL,
    1099:                                        NULL,
    1100:                                        NULL,
    1101:                                 /*94*/ NULL,
    1102:                                        NULL,
    1103:                                        NULL,
    1104:                                        NULL,
    1105:                                        NULL,
    1106:                                 /*99*/ NULL,
    1107:                                        NULL,
    1108:                                        NULL,
    1109:                                        MatConjugate_MPISELL,
    1110:                                        NULL,
    1111:                                 /*104*/NULL,
    1112:                                        MatRealPart_MPISELL,
    1113:                                        MatImaginaryPart_MPISELL,
    1114:                                        NULL,
    1115:                                        NULL,
    1116:                                 /*109*/NULL,
    1117:                                        NULL,
    1118:                                        NULL,
    1119:                                        NULL,
    1120:                                        MatMissingDiagonal_MPISELL,
    1121:                                 /*114*/NULL,
    1122:                                        NULL,
    1123:                                        MatGetGhosts_MPISELL,
    1124:                                        NULL,
    1125:                                        NULL,
    1126:                                 /*119*/NULL,
    1127:                                        NULL,
    1128:                                        NULL,
    1129:                                        NULL,
    1130:                                        NULL,
    1131:                                 /*124*/NULL,
    1132:                                        NULL,
    1133:                                        MatInvertBlockDiagonal_MPISELL,
    1134:                                        NULL,
    1135:                                        NULL,
    1136:                                 /*129*/NULL,
    1137:                                        NULL,
    1138:                                        NULL,
    1139:                                        NULL,
    1140:                                        NULL,
    1141:                                 /*134*/NULL,
    1142:                                        NULL,
    1143:                                        NULL,
    1144:                                        NULL,
    1145:                                        NULL,
    1146:                                 /*139*/NULL,
    1147:                                        NULL,
    1148:                                        NULL,
    1149:                                        MatFDColoringSetUp_MPIXAIJ,
    1150:                                        NULL,
    1151:                                 /*144*/NULL
    1152: };

    1154: /* ----------------------------------------------------------------------------------------*/

    1156: PetscErrorCode MatStoreValues_MPISELL(Mat mat)
    1157: {
    1158:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

    1162:   MatStoreValues(sell->A);
    1163:   MatStoreValues(sell->B);
    1164:   return(0);
    1165: }

    1167: PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
    1168: {
    1169:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

    1173:   MatRetrieveValues(sell->A);
    1174:   MatRetrieveValues(sell->B);
    1175:   return(0);
    1176: }

    1178: PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[])
    1179: {
    1180:   Mat_MPISELL    *b;

    1184:   PetscLayoutSetUp(B->rmap);
    1185:   PetscLayoutSetUp(B->cmap);
    1186:   b = (Mat_MPISELL*)B->data;

    1188:   if (!B->preallocated) {
    1189:     /* Explicitly create 2 MATSEQSELL matrices. */
    1190:     MatCreate(PETSC_COMM_SELF,&b->A);
    1191:     MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
    1192:     MatSetBlockSizesFromMats(b->A,B,B);
    1193:     MatSetType(b->A,MATSEQSELL);
    1194:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
    1195:     MatCreate(PETSC_COMM_SELF,&b->B);
    1196:     MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
    1197:     MatSetBlockSizesFromMats(b->B,B,B);
    1198:     MatSetType(b->B,MATSEQSELL);
    1199:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
    1200:   }

    1202:   MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);
    1203:   MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);
    1204:   B->preallocated  = PETSC_TRUE;
    1205:   B->was_assembled = PETSC_FALSE;
    1206:   /*
    1207:     critical for MatAssemblyEnd to work.
    1208:     MatAssemblyBegin checks it to set up was_assembled
    1209:     and MatAssemblyEnd checks was_assembled to determine whether to build garray
    1210:   */
    1211:   B->assembled     = PETSC_FALSE;
    1212:   return(0);
    1213: }

    1215: PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat)
    1216: {
    1217:   Mat            mat;
    1218:   Mat_MPISELL    *a,*oldmat=(Mat_MPISELL*)matin->data;

    1222:   *newmat = NULL;
    1223:   MatCreate(PetscObjectComm((PetscObject)matin),&mat);
    1224:   MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);
    1225:   MatSetBlockSizesFromMats(mat,matin,matin);
    1226:   MatSetType(mat,((PetscObject)matin)->type_name);
    1227:   a       = (Mat_MPISELL*)mat->data;

    1229:   mat->factortype   = matin->factortype;
    1230:   mat->assembled    = PETSC_TRUE;
    1231:   mat->insertmode   = NOT_SET_VALUES;
    1232:   mat->preallocated = PETSC_TRUE;

    1234:   a->size         = oldmat->size;
    1235:   a->rank         = oldmat->rank;
    1236:   a->donotstash   = oldmat->donotstash;
    1237:   a->roworiented  = oldmat->roworiented;
    1238:   a->rowindices   = NULL;
    1239:   a->rowvalues    = NULL;
    1240:   a->getrowactive = PETSC_FALSE;

    1242:   PetscLayoutReference(matin->rmap,&mat->rmap);
    1243:   PetscLayoutReference(matin->cmap,&mat->cmap);

    1245:   if (oldmat->colmap) {
    1246: #if defined(PETSC_USE_CTABLE)
    1247:     PetscTableCreateCopy(oldmat->colmap,&a->colmap);
    1248: #else
    1249:     PetscMalloc1(mat->cmap->N,&a->colmap);
    1250:     PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));
    1251:     PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N);
    1252: #endif
    1253:   } else a->colmap = NULL;
    1254:   if (oldmat->garray) {
    1255:     PetscInt len;
    1256:     len  = oldmat->B->cmap->n;
    1257:     PetscMalloc1(len+1,&a->garray);
    1258:     PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));
    1259:     if (len) { PetscArraycpy(a->garray,oldmat->garray,len); }
    1260:   } else a->garray = NULL;

    1262:   VecDuplicate(oldmat->lvec,&a->lvec);
    1263:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);
    1264:   VecScatterCopy(oldmat->Mvctx,&a->Mvctx);
    1265:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);
    1266:   MatDuplicate(oldmat->A,cpvalues,&a->A);
    1267:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);
    1268:   MatDuplicate(oldmat->B,cpvalues,&a->B);
    1269:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);
    1270:   PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);
    1271:   *newmat = mat;
    1272:   return(0);
    1273: }

    1275: /*@C
    1276:    MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format.
    1277:    For good matrix assembly performance the user should preallocate the matrix storage by
    1278:    setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz).

    1280:    Collective

    1282:    Input Parameters:
    1283: +  B - the matrix
    1284: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
    1285:            (same value is used for all local rows)
    1286: .  d_nnz - array containing the number of nonzeros in the various rows of the
    1287:            DIAGONAL portion of the local submatrix (possibly different for each row)
    1288:            or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
    1289:            The size of this array is equal to the number of local rows, i.e 'm'.
    1290:            For matrices that will be factored, you must leave room for (and set)
    1291:            the diagonal entry even if it is zero.
    1292: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
    1293:            submatrix (same value is used for all local rows).
    1294: -  o_nnz - array containing the number of nonzeros in the various rows of the
    1295:            OFF-DIAGONAL portion of the local submatrix (possibly different for
    1296:            each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
    1297:            structure. The size of this array is equal to the number
    1298:            of local rows, i.e 'm'.

    1300:    If the *_nnz parameter is given then the *_nz parameter is ignored

    1302:    The stored row and column indices begin with zero.

    1304:    The parallel matrix is partitioned such that the first m0 rows belong to
    1305:    process 0, the next m1 rows belong to process 1, the next m2 rows belong
    1306:    to process 2 etc.. where m0,m1,m2... are the input parameter 'm'.

    1308:    The DIAGONAL portion of the local submatrix of a processor can be defined
    1309:    as the submatrix which is obtained by extraction the part corresponding to
    1310:    the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
    1311:    first row that belongs to the processor, r2 is the last row belonging to
    1312:    the this processor, and c1-c2 is range of indices of the local part of a
    1313:    vector suitable for applying the matrix to.  This is an mxn matrix.  In the
    1314:    common case of a square matrix, the row and column ranges are the same and
    1315:    the DIAGONAL part is also square. The remaining portion of the local
    1316:    submatrix (mxN) constitute the OFF-DIAGONAL portion.

    1318:    If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored.

    1320:    You can call MatGetInfo() to get information on how effective the preallocation was;
    1321:    for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
    1322:    You can also run with the option -info and look for messages with the string
    1323:    malloc in them to see if additional memory allocation was needed.

    1325:    Example usage:

    1327:    Consider the following 8x8 matrix with 34 non-zero values, that is
    1328:    assembled across 3 processors. Lets assume that proc0 owns 3 rows,
    1329:    proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
    1330:    as follows:

    1332: .vb
    1333:             1  2  0  |  0  3  0  |  0  4
    1334:     Proc0   0  5  6  |  7  0  0  |  8  0
    1335:             9  0 10  | 11  0  0  | 12  0
    1336:     -------------------------------------
    1337:            13  0 14  | 15 16 17  |  0  0
    1338:     Proc1   0 18  0  | 19 20 21  |  0  0
    1339:             0  0  0  | 22 23  0  | 24  0
    1340:     -------------------------------------
    1341:     Proc2  25 26 27  |  0  0 28  | 29  0
    1342:            30  0  0  | 31 32 33  |  0 34
    1343: .ve

    1345:    This can be represented as a collection of submatrices as:

    1347: .vb
    1348:       A B C
    1349:       D E F
    1350:       G H I
    1351: .ve

    1353:    Where the submatrices A,B,C are owned by proc0, D,E,F are
    1354:    owned by proc1, G,H,I are owned by proc2.

    1356:    The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
    1357:    The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
    1358:    The 'M','N' parameters are 8,8, and have the same values on all procs.

    1360:    The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
    1361:    submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
    1362:    corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
    1363:    Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
    1364:    part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
    1365:    matrix, ans [DF] as another SeqSELL matrix.

    1367:    When d_nz, o_nz parameters are specified, d_nz storage elements are
    1368:    allocated for every row of the local diagonal submatrix, and o_nz
    1369:    storage locations are allocated for every row of the OFF-DIAGONAL submat.
    1370:    One way to choose d_nz and o_nz is to use the max nonzerors per local
    1371:    rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
    1372:    In this case, the values of d_nz,o_nz are:
    1373: .vb
    1374:      proc0 : dnz = 2, o_nz = 2
    1375:      proc1 : dnz = 3, o_nz = 2
    1376:      proc2 : dnz = 1, o_nz = 4
    1377: .ve
    1378:    We are allocating m*(d_nz+o_nz) storage locations for every proc. This
    1379:    translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
    1380:    for proc3. i.e we are using 12+15+10=37 storage locations to store
    1381:    34 values.

    1383:    When d_nnz, o_nnz parameters are specified, the storage is specified
    1384:    for every row, coresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
    1385:    In the above case the values for d_nnz,o_nnz are:
    1386: .vb
    1387:      proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
    1388:      proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
    1389:      proc2: d_nnz = [1,1]   and o_nnz = [4,4]
    1390: .ve
    1391:    Here the space allocated is according to nz (or maximum values in the nnz
    1392:    if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37

    1394:    Level: intermediate

    1396: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(),
    1397:           MATMPISELL, MatGetInfo(), PetscSplitOwnership()
    1398: @*/
    1399: PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
    1400: {

    1406:   PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
    1407:   return(0);
    1408: }

    1410: /*@C
    1411:    MatCreateSELL - Creates a sparse parallel matrix in SELL format.

    1413:    Collective

    1415:    Input Parameters:
    1416: +  comm - MPI communicator
    1417: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
    1418:            This value should be the same as the local size used in creating the
    1419:            y vector for the matrix-vector product y = Ax.
    1420: .  n - This value should be the same as the local size used in creating the
    1421:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
    1422:        calculated if N is given) For square matrices n is almost always m.
    1423: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
    1424: .  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
    1425: .  d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix
    1426:                (same value is used for all local rows)
    1427: .  d_rlen - array containing the number of nonzeros in the various rows of the
    1428:             DIAGONAL portion of the local submatrix (possibly different for each row)
    1429:             or NULL, if d_rlenmax is used to specify the nonzero structure.
    1430:             The size of this array is equal to the number of local rows, i.e 'm'.
    1431: .  o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local
    1432:                submatrix (same value is used for all local rows).
    1433: -  o_rlen - array containing the number of nonzeros in the various rows of the
    1434:             OFF-DIAGONAL portion of the local submatrix (possibly different for
    1435:             each row) or NULL, if o_rlenmax is used to specify the nonzero
    1436:             structure. The size of this array is equal to the number
    1437:             of local rows, i.e 'm'.

    1439:    Output Parameter:
    1440: .  A - the matrix

    1442:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
    1443:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
    1444:    [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation]

    1446:    Notes:
    1447:    If the *_rlen parameter is given then the *_rlenmax parameter is ignored

    1449:    m,n,M,N parameters specify the size of the matrix, and its partitioning across
    1450:    processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate
    1451:    storage requirements for this matrix.

    1453:    If PETSC_DECIDE or  PETSC_DETERMINE is used for a particular argument on one
    1454:    processor than it must be used on all processors that share the object for
    1455:    that argument.

    1457:    The user MUST specify either the local or global matrix dimensions
    1458:    (possibly both).

    1460:    The parallel matrix is partitioned across processors such that the
    1461:    first m0 rows belong to process 0, the next m1 rows belong to
    1462:    process 1, the next m2 rows belong to process 2 etc.. where
    1463:    m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores
    1464:    values corresponding to [m x N] submatrix.

    1466:    The columns are logically partitioned with the n0 columns belonging
    1467:    to 0th partition, the next n1 columns belonging to the next
    1468:    partition etc.. where n0,n1,n2... are the input parameter 'n'.

    1470:    The DIAGONAL portion of the local submatrix on any given processor
    1471:    is the submatrix corresponding to the rows and columns m,n
    1472:    corresponding to the given processor. i.e diagonal matrix on
    1473:    process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
    1474:    etc. The remaining portion of the local submatrix [m x (N-n)]
    1475:    constitute the OFF-DIAGONAL portion. The example below better
    1476:    illustrates this concept.

    1478:    For a square global matrix we define each processor's diagonal portion
    1479:    to be its local rows and the corresponding columns (a square submatrix);
    1480:    each processor's off-diagonal portion encompasses the remainder of the
    1481:    local matrix (a rectangular submatrix).

    1483:    If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored.

    1485:    When calling this routine with a single process communicator, a matrix of
    1486:    type SEQSELL is returned.  If a matrix of type MATMPISELL is desired for this
    1487:    type of communicator, use the construction mechanism:
    1488:      MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...);

    1490:    Options Database Keys:
    1491: -  -mat_sell_oneindex - Internally use indexing starting at 1
    1492:         rather than 0.  Note that when calling MatSetValues(),
    1493:         the user still MUST index entries starting at 0!

    1496:    Example usage:

    1498:    Consider the following 8x8 matrix with 34 non-zero values, that is
    1499:    assembled across 3 processors. Lets assume that proc0 owns 3 rows,
    1500:    proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
    1501:    as follows:

    1503: .vb
    1504:             1  2  0  |  0  3  0  |  0  4
    1505:     Proc0   0  5  6  |  7  0  0  |  8  0
    1506:             9  0 10  | 11  0  0  | 12  0
    1507:     -------------------------------------
    1508:            13  0 14  | 15 16 17  |  0  0
    1509:     Proc1   0 18  0  | 19 20 21  |  0  0
    1510:             0  0  0  | 22 23  0  | 24  0
    1511:     -------------------------------------
    1512:     Proc2  25 26 27  |  0  0 28  | 29  0
    1513:            30  0  0  | 31 32 33  |  0 34
    1514: .ve

    1516:    This can be represented as a collection of submatrices as:

    1518: .vb
    1519:       A B C
    1520:       D E F
    1521:       G H I
    1522: .ve

    1524:    Where the submatrices A,B,C are owned by proc0, D,E,F are
    1525:    owned by proc1, G,H,I are owned by proc2.

    1527:    The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
    1528:    The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
    1529:    The 'M','N' parameters are 8,8, and have the same values on all procs.

    1531:    The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
    1532:    submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
    1533:    corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
    1534:    Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
    1535:    part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
    1536:    matrix, ans [DF] as another SeqSELL matrix.

    1538:    When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are
    1539:    allocated for every row of the local diagonal submatrix, and o_rlenmax
    1540:    storage locations are allocated for every row of the OFF-DIAGONAL submat.
    1541:    One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local
    1542:    rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
    1543:    In this case, the values of d_rlenmax,o_rlenmax are:
    1544: .vb
    1545:      proc0 : d_rlenmax = 2, o_rlenmax = 2
    1546:      proc1 : d_rlenmax = 3, o_rlenmax = 2
    1547:      proc2 : d_rlenmax = 1, o_rlenmax = 4
    1548: .ve
    1549:    We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This
    1550:    translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
    1551:    for proc3. i.e we are using 12+15+10=37 storage locations to store
    1552:    34 values.

    1554:    When d_rlen, o_rlen parameters are specified, the storage is specified
    1555:    for every row, coresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
    1556:    In the above case the values for d_nnz,o_nnz are:
    1557: .vb
    1558:      proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
    1559:      proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
    1560:      proc2: d_nnz = [1,1]   and o_nnz = [4,4]
    1561: .ve
    1562:    Here the space allocated is still 37 though there are 34 nonzeros because
    1563:    the allocation is always done according to rlenmax.

    1565:    Level: intermediate

    1567: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(),
    1568:           MATMPISELL, MatCreateMPISELLWithArrays()
    1569: @*/
    1570: PetscErrorCode MatCreateSELL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[],Mat *A)
    1571: {
    1573:   PetscMPIInt    size;

    1576:   MatCreate(comm,A);
    1577:   MatSetSizes(*A,m,n,M,N);
    1578:   MPI_Comm_size(comm,&size);
    1579:   if (size > 1) {
    1580:     MatSetType(*A,MATMPISELL);
    1581:     MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);
    1582:   } else {
    1583:     MatSetType(*A,MATSEQSELL);
    1584:     MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);
    1585:   }
    1586:   return(0);
    1587: }

    1589: PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[])
    1590: {
    1591:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
    1592:   PetscBool      flg;

    1596:   PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);
    1597:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input");
    1598:   if (Ad)     *Ad     = a->A;
    1599:   if (Ao)     *Ao     = a->B;
    1600:   if (colmap) *colmap = a->garray;
    1601:   return(0);
    1602: }

    1604: /*@C
    1605:      MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns

    1607:     Not Collective

    1609:    Input Parameters:
    1610: +    A - the matrix
    1611: .    scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
    1612: -    row, col - index sets of rows and columns to extract (or NULL)

    1614:    Output Parameter:
    1615: .    A_loc - the local sequential matrix generated

    1617:     Level: developer

    1619: .seealso: MatGetOwnershipRange(), MatMPISELLGetLocalMat()

    1621: @*/
    1622: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc)
    1623: {
    1624:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
    1626:   PetscInt       i,start,end,ncols,nzA,nzB,*cmap,imark,*idx;
    1627:   IS             isrowa,iscola;
    1628:   Mat            *aloc;
    1629:   PetscBool      match;

    1632:   PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match);
    1633:   if (!match) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input");
    1634:   PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0);
    1635:   if (!row) {
    1636:     start = A->rmap->rstart; end = A->rmap->rend;
    1637:     ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa);
    1638:   } else {
    1639:     isrowa = *row;
    1640:   }
    1641:   if (!col) {
    1642:     start = A->cmap->rstart;
    1643:     cmap  = a->garray;
    1644:     nzA   = a->A->cmap->n;
    1645:     nzB   = a->B->cmap->n;
    1646:     PetscMalloc1(nzA+nzB, &idx);
    1647:     ncols = 0;
    1648:     for (i=0; i<nzB; i++) {
    1649:       if (cmap[i] < start) idx[ncols++] = cmap[i];
    1650:       else break;
    1651:     }
    1652:     imark = i;
    1653:     for (i=0; i<nzA; i++) idx[ncols++] = start + i;
    1654:     for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i];
    1655:     ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola);
    1656:   } else {
    1657:     iscola = *col;
    1658:   }
    1659:   if (scall != MAT_INITIAL_MATRIX) {
    1660:     PetscMalloc1(1,&aloc);
    1661:     aloc[0] = *A_loc;
    1662:   }
    1663:   MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc);
    1664:   *A_loc = aloc[0];
    1665:   PetscFree(aloc);
    1666:   if (!row) {
    1667:     ISDestroy(&isrowa);
    1668:   }
    1669:   if (!col) {
    1670:     ISDestroy(&iscola);
    1671:   }
    1672:   PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0);
    1673:   return(0);
    1674: }

    1676: #include <../src/mat/impls/aij/mpi/mpiaij.h>

    1678: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
    1679: {
    1681:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
    1682:   Mat            B;
    1683:   Mat_MPIAIJ     *b;

    1686:   if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled");

    1688:   if (reuse == MAT_REUSE_MATRIX) {
    1689:     B = *newmat;
    1690:   } else {
    1691:     MatCreate(PetscObjectComm((PetscObject)A),&B);
    1692:     MatSetType(B,MATMPIAIJ);
    1693:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
    1694:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
    1695:     MatSeqAIJSetPreallocation(B,0,NULL);
    1696:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
    1697:   }
    1698:   b    = (Mat_MPIAIJ*) B->data;

    1700:   if (reuse == MAT_REUSE_MATRIX) {
    1701:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);
    1702:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);
    1703:   } else {
    1704:     MatDestroy(&b->A);
    1705:     MatDestroy(&b->B);
    1706:     MatDisAssemble_MPISELL(A);
    1707:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);
    1708:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);
    1709:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    1710:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    1711:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    1712:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    1713:   }

    1715:   if (reuse == MAT_INPLACE_MATRIX) {
    1716:     MatHeaderReplace(A,&B);
    1717:   } else {
    1718:     *newmat = B;
    1719:   }
    1720:   return(0);
    1721: }

    1723: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
    1724: {
    1726:   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)A->data;
    1727:   Mat            B;
    1728:   Mat_MPISELL    *b;

    1731:   if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled");

    1733:   if (reuse == MAT_REUSE_MATRIX) {
    1734:     B = *newmat;
    1735:   } else {
    1736:     MatCreate(PetscObjectComm((PetscObject)A),&B);
    1737:     MatSetType(B,MATMPISELL);
    1738:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
    1739:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
    1740:     MatSeqAIJSetPreallocation(B,0,NULL);
    1741:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
    1742:   }
    1743:   b    = (Mat_MPISELL*) B->data;

    1745:   if (reuse == MAT_REUSE_MATRIX) {
    1746:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);
    1747:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);
    1748:   } else {
    1749:     MatDestroy(&b->A);
    1750:     MatDestroy(&b->B);
    1751:     MatDisAssemble_MPIAIJ(A);
    1752:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);
    1753:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);
    1754:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    1755:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    1756:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    1757:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    1758:   }

    1760:   if (reuse == MAT_INPLACE_MATRIX) {
    1761:     MatHeaderReplace(A,&B);
    1762:   } else {
    1763:     *newmat = B;
    1764:   }
    1765:   return(0);
    1766: }

    1768: PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
    1769: {
    1770:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
    1772:   Vec            bb1=NULL;

    1775:   if (flag == SOR_APPLY_UPPER) {
    1776:     (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
    1777:     return(0);
    1778:   }

    1780:   if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) {
    1781:     VecDuplicate(bb,&bb1);
    1782:   }

    1784:   if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
    1785:     if (flag & SOR_ZERO_INITIAL_GUESS) {
    1786:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
    1787:       its--;
    1788:     }

    1790:     while (its--) {
    1791:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
    1792:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

    1794:       /* update rhs: bb1 = bb - B*x */
    1795:       VecScale(mat->lvec,-1.0);
    1796:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

    1798:       /* local sweep */
    1799:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);
    1800:     }
    1801:   } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
    1802:     if (flag & SOR_ZERO_INITIAL_GUESS) {
    1803:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
    1804:       its--;
    1805:     }
    1806:     while (its--) {
    1807:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
    1808:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

    1810:       /* update rhs: bb1 = bb - B*x */
    1811:       VecScale(mat->lvec,-1.0);
    1812:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

    1814:       /* local sweep */
    1815:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);
    1816:     }
    1817:   } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
    1818:     if (flag & SOR_ZERO_INITIAL_GUESS) {
    1819:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
    1820:       its--;
    1821:     }
    1822:     while (its--) {
    1823:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
    1824:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

    1826:       /* update rhs: bb1 = bb - B*x */
    1827:       VecScale(mat->lvec,-1.0);
    1828:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

    1830:       /* local sweep */
    1831:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx);
    1832:     }
    1833:   } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported");

    1835:   VecDestroy(&bb1);

    1837:   matin->factorerrortype = mat->A->factorerrortype;
    1838:   return(0);
    1839: }

    1841: /*MC
    1842:    MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices.

    1844:    Options Database Keys:
    1845: . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions()

    1847:   Level: beginner

    1849: .seealso: MatCreateSELL()
    1850: M*/
    1851: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
    1852: {
    1853:   Mat_MPISELL    *b;
    1855:   PetscMPIInt    size;

    1858:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
    1859:   PetscNewLog(B,&b);
    1860:   B->data       = (void*)b;
    1861:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
    1862:   B->assembled  = PETSC_FALSE;
    1863:   B->insertmode = NOT_SET_VALUES;
    1864:   b->size       = size;
    1865:   MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);
    1866:   /* build cache for off array entries formed */
    1867:   MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);

    1869:   b->donotstash  = PETSC_FALSE;
    1870:   b->colmap      = NULL;
    1871:   b->garray      = NULL;
    1872:   b->roworiented = PETSC_TRUE;

    1874:   /* stuff used for matrix vector multiply */
    1875:   b->lvec  = NULL;
    1876:   b->Mvctx = NULL;

    1878:   /* stuff for MatGetRow() */
    1879:   b->rowindices   = NULL;
    1880:   b->rowvalues    = NULL;
    1881:   b->getrowactive = PETSC_FALSE;

    1883:   PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);
    1884:   PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);
    1885:   PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);
    1886:   PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);
    1887:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);
    1888:   PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);
    1889:   PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);
    1890:   return(0);
    1891: }