Actual source code: mpisell.c

petsc-master 2020-07-10
Report Typos and Errors
  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:   Developer Notes:
 21:     Subclasses include MATSELLCUSP, MATSELLCUSPARSE, MATSELLPERM, MATSELLCRL, and also automatically switches over to use inodes when
 22:    enough exist.

 24:   Level: beginner

 26: .seealso: MatCreateSELL(), MatCreateSeqSELL(), MATSEQSELL, MATMPISELL
 27: M*/

 29: PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is)
 30: {
 32:   Mat_MPISELL    *sell=(Mat_MPISELL*)Y->data;

 35:   if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) {
 36:     MatDiagonalSet(sell->A,D,is);
 37:   } else {
 38:     MatDiagonalSet_Default(Y,D,is);
 39:   }
 40:   return(0);
 41: }

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

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

 71: #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \
 72:   { \
 73:     if (col <= lastcol1) low1 = 0; \
 74:     else                high1 = nrow1; \
 75:     lastcol1 = col; \
 76:     while (high1-low1 > 5) { \
 77:       t = (low1+high1)/2; \
 78:       if (*(cp1+8*t) > col) high1 = t; \
 79:       else                   low1 = t; \
 80:     } \
 81:     for (_i=low1; _i<high1; _i++) { \
 82:       if (*(cp1+8*_i) > col) break; \
 83:       if (*(cp1+8*_i) == col) { \
 84:         if (addv == ADD_VALUES) *(vp1+8*_i) += value;   \
 85:         else                     *(vp1+8*_i) = value; \
 86:         goto a_noinsert; \
 87:       } \
 88:     }  \
 89:     if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \
 90:     if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \
 91:     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); \
 92:     MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \
 93:     /* shift up all the later entries in this row */ \
 94:     for (ii=nrow1-1; ii>=_i; ii--) { \
 95:       *(cp1+8*(ii+1)) = *(cp1+8*ii); \
 96:       *(vp1+8*(ii+1)) = *(vp1+8*ii); \
 97:     } \
 98:     *(cp1+8*_i) = col; \
 99:     *(vp1+8*_i) = value; \
100:     a->nz++; nrow1++; A->nonzerostate++; \
101:     a_noinsert: ; \
102:     a->rlen[row] = nrow1; \
103:   }

105: #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \
106:   { \
107:     if (col <= lastcol2) low2 = 0; \
108:     else                high2 = nrow2; \
109:     lastcol2 = col; \
110:     while (high2-low2 > 5) { \
111:       t = (low2+high2)/2; \
112:       if (*(cp2+8*t) > col) high2 = t; \
113:       else low2  = t; \
114:     } \
115:     for (_i=low2; _i<high2; _i++) { \
116:       if (*(cp2+8*_i) > col) break; \
117:       if (*(cp2+8*_i) == col) { \
118:         if (addv == ADD_VALUES) *(vp2+8*_i) += value; \
119:         else                     *(vp2+8*_i) = value; \
120:         goto b_noinsert; \
121:       } \
122:     } \
123:     if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
124:     if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
125:     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); \
126:     MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \
127:     /* shift up all the later entries in this row */ \
128:     for (ii=nrow2-1; ii>=_i; ii--) { \
129:       *(cp2+8*(ii+1)) = *(cp2+8*ii); \
130:       *(vp2+8*(ii+1)) = *(vp2+8*ii); \
131:     } \
132:     *(cp2+8*_i) = col; \
133:     *(vp2+8*_i) = value; \
134:     b->nz++; nrow2++; B->nonzerostate++; \
135:     b_noinsert: ; \
136:     b->rlen[row] = nrow2; \
137:   }

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

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

158:   for (i=0; i<m; i++) {
159:     if (im[i] < 0) continue;
160:     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);
161:     if (im[i] >= rstart && im[i] < rend) {
162:       row      = im[i] - rstart;
163:       lastcol1 = -1;
164:       shift1   = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */
165:       cp1      = a->colidx+shift1;
166:       vp1      = a->val+shift1;
167:       nrow1    = a->rlen[row];
168:       low1     = 0;
169:       high1    = nrow1;
170:       lastcol2 = -1;
171:       shift2   = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
172:       cp2      = b->colidx+shift2;
173:       vp2      = b->val+shift2;
174:       nrow2    = b->rlen[row];
175:       low2     = 0;
176:       high2    = nrow2;

178:       for (j=0; j<n; j++) {
179:         if (roworiented) value = v[i*n+j];
180:         else             value = v[i+j*m];
181:         if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue;
182:         if (in[j] >= cstart && in[j] < cend) {
183:           col   = in[j] - cstart;
184:           MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */
185:         } else if (in[j] < 0) continue;
186:         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);
187:         else {
188:           if (mat->was_assembled) {
189:             if (!sell->colmap) {
190:               MatCreateColmap_MPISELL_Private(mat);
191:             }
192: #if defined(PETSC_USE_CTABLE)
193:             PetscTableFind(sell->colmap,in[j]+1,&col);
194:             col--;
195: #else
196:             col = sell->colmap[in[j]] - 1;
197: #endif
198:             if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) {
199:               MatDisAssemble_MPISELL(mat);
200:               col    = in[j];
201:               /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */
202:               B      = sell->B;
203:               b      = (Mat_SeqSELL*)B->data;
204:               shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
205:               cp2    = b->colidx+shift2;
206:               vp2    = b->val+shift2;
207:               nrow2  = b->rlen[row];
208:               low2   = 0;
209:               high2  = nrow2;
210:             } 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]);
211:           } else col = in[j];
212:           MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */
213:         }
214:       }
215:     } else {
216:       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]);
217:       if (!sell->donotstash) {
218:         mat->assembled = PETSC_FALSE;
219:         if (roworiented) {
220:           MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
221:         } else {
222:           MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
223:         }
224:       }
225:     }
226:   }
227:   return(0);
228: }

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

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

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

272: PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode)
273: {
274:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
276:   PetscInt       nstash,reallocs;

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

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

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

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

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

339:   /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
340:   if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) {
341:     PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate;
342:     MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));
343:   }
344:   return(0);
345: }

347: PetscErrorCode MatZeroEntries_MPISELL(Mat A)
348: {
349:   Mat_MPISELL    *l=(Mat_MPISELL*)A->data;

353:   MatZeroEntries(l->A);
354:   MatZeroEntries(l->B);
355:   return(0);
356: }

358: PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy)
359: {
360:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
362:   PetscInt       nt;

365:   VecGetLocalSize(xx,&nt);
366:   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);
367:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
368:   (*a->A->ops->mult)(a->A,xx,yy);
369:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
370:   (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);
371:   return(0);
372: }

374: PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx)
375: {
376:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

380:   MatMultDiagonalBlock(a->A,bb,xx);
381:   return(0);
382: }

384: PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
385: {
386:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

390:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
391:   (*a->A->ops->multadd)(a->A,xx,yy,zz);
392:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
393:   (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
394:   return(0);
395: }

397: PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy)
398: {
399:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

403:   /* do nondiagonal part */
404:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
405:   /* do local part */
406:   (*a->A->ops->multtranspose)(a->A,xx,yy);
407:   /* add partial results together */
408:   VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
409:   VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
410:   return(0);
411: }

413: PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f)
414: {
415:   MPI_Comm       comm;
416:   Mat_MPISELL    *Asell=(Mat_MPISELL*)Amat->data,*Bsell;
417:   Mat            Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs;
418:   IS             Me,Notme;
420:   PetscInt       M,N,first,last,*notme,i;
421:   PetscMPIInt    size;

424:   /* Easy test: symmetric diagonal block */
425:   Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A;
426:   MatIsTranspose(Adia,Bdia,tol,f);
427:   if (!*f) return(0);
428:   PetscObjectGetComm((PetscObject)Amat,&comm);
429:   MPI_Comm_size(comm,&size);
430:   if (size == 1) return(0);

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

453: PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
454: {
455:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

459:   /* do nondiagonal part */
460:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
461:   /* do local part */
462:   (*a->A->ops->multtransposeadd)(a->A,xx,yy,zz);
463:   /* add partial results together */
464:   VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
465:   VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
466:   return(0);
467: }

469: /*
470:   This only works correctly for square matrices where the subblock A->A is the
471:    diagonal block
472: */
473: PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v)
474: {
476:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

479:   if (A->rmap->N != A->cmap->N) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block");
480:   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");
481:   MatGetDiagonal(a->A,v);
482:   return(0);
483: }

485: PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa)
486: {
487:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

491:   MatScale(a->A,aa);
492:   MatScale(a->B,aa);
493:   return(0);
494: }

496: PetscErrorCode MatDestroy_MPISELL(Mat mat)
497: {
498:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

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

521:   PetscObjectChangeTypeName((PetscObject)mat,0);
522:   PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);
523:   PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);
524:   PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL);
525:   PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL);
526:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL);
527:   PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);
528:   return(0);
529: }

531:  #include <petscdraw.h>
532: PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer)
533: {
534:   Mat_MPISELL       *sell=(Mat_MPISELL*)mat->data;
535:   PetscErrorCode    ierr;
536:   PetscMPIInt       rank=sell->rank,size=sell->size;
537:   PetscBool         isdraw,iascii,isbinary;
538:   PetscViewer       sviewer;
539:   PetscViewerFormat format;

542:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
543:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
544:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
545:   if (iascii) {
546:     PetscViewerGetFormat(viewer,&format);
547:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
548:       MatInfo   info;
549:       PetscBool inodes;

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

599:   {
600:     /* assemble the entire matrix onto first processor. */
601:     Mat         A;
602:     Mat_SeqSELL *Aloc;
603:     PetscInt    M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j;
604:     MatScalar   *aval;
605:     PetscBool   isnonzero;

607:     MatCreate(PetscObjectComm((PetscObject)mat),&A);
608:     if (!rank) {
609:       MatSetSizes(A,M,N,M,N);
610:     } else {
611:       MatSetSizes(A,0,0,M,N);
612:     }
613:     /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */
614:     MatSetType(A,MATMPISELL);
615:     MatMPISELLSetPreallocation(A,0,NULL,0,NULL);
616:     MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
617:     PetscLogObjectParent((PetscObject)mat,(PetscObject)A);

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

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

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

667: PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer)
668: {
670:   PetscBool      iascii,isdraw,issocket,isbinary;

673:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
674:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
675:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
676:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
677:   if (iascii || isdraw || isbinary || issocket) {
678:     MatView_MPISELL_ASCIIorDraworSocket(mat,viewer);
679:   }
680:   return(0);
681: }

683: PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
684: {
685:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

689:   MatGetSize(sell->B,NULL,nghosts);
690:   if (ghosts) *ghosts = sell->garray;
691:   return(0);
692: }

694: PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info)
695: {
696:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
697:   Mat            A=mat->A,B=mat->B;
699:   PetscLogDouble isend[5],irecv[5];

702:   info->block_size = 1.0;
703:   MatGetInfo(A,MAT_LOCAL,info);

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

708:   MatGetInfo(B,MAT_LOCAL,info);

710:   isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded;
711:   isend[3] += info->memory;  isend[4] += info->mallocs;
712:   if (flag == MAT_LOCAL) {
713:     info->nz_used      = isend[0];
714:     info->nz_allocated = isend[1];
715:     info->nz_unneeded  = isend[2];
716:     info->memory       = isend[3];
717:     info->mallocs      = isend[4];
718:   } else if (flag == MAT_GLOBAL_MAX) {
719:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin));

721:     info->nz_used      = irecv[0];
722:     info->nz_allocated = irecv[1];
723:     info->nz_unneeded  = irecv[2];
724:     info->memory       = irecv[3];
725:     info->mallocs      = irecv[4];
726:   } else if (flag == MAT_GLOBAL_SUM) {
727:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin));

729:     info->nz_used      = irecv[0];
730:     info->nz_allocated = irecv[1];
731:     info->nz_unneeded  = irecv[2];
732:     info->memory       = irecv[3];
733:     info->mallocs      = irecv[4];
734:   }
735:   info->fill_ratio_given  = 0; /* no parallel LU/ILU/Cholesky */
736:   info->fill_ratio_needed = 0;
737:   info->factor_mallocs    = 0;
738:   return(0);
739: }

741: PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg)
742: {
743:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

747:   switch (op) {
748:   case MAT_NEW_NONZERO_LOCATIONS:
749:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
750:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
751:   case MAT_KEEP_NONZERO_PATTERN:
752:   case MAT_NEW_NONZERO_LOCATION_ERR:
753:   case MAT_USE_INODES:
754:   case MAT_IGNORE_ZERO_ENTRIES:
755:     MatCheckPreallocated(A,1);
756:     MatSetOption(a->A,op,flg);
757:     MatSetOption(a->B,op,flg);
758:     break;
759:   case MAT_ROW_ORIENTED:
760:     MatCheckPreallocated(A,1);
761:     a->roworiented = flg;

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


806: PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr)
807: {
808:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
809:   Mat            a=sell->A,b=sell->B;
811:   PetscInt       s1,s2,s3;

814:   MatGetLocalSize(mat,&s2,&s3);
815:   if (rr) {
816:     VecGetLocalSize(rr,&s1);
817:     if (s1!=s3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size");
818:     /* Overlap communication with computation. */
819:     VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
820:   }
821:   if (ll) {
822:     VecGetLocalSize(ll,&s1);
823:     if (s1!=s2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size");
824:     (*b->ops->diagonalscale)(b,ll,0);
825:   }
826:   /* scale  the diagonal block */
827:   (*a->ops->diagonalscale)(a,ll,rr);

829:   if (rr) {
830:     /* Do a scatter end and then right scale the off-diagonal block */
831:     VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
832:     (*b->ops->diagonalscale)(b,0,sell->lvec);
833:   }
834:   return(0);
835: }

837: PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
838: {
839:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

843:   MatSetUnfactored(a->A);
844:   return(0);
845: }

847: PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool  *flag)
848: {
849:   Mat_MPISELL    *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data;
850:   Mat            a,b,c,d;
851:   PetscBool      flg;

855:   a = matA->A; b = matA->B;
856:   c = matB->A; d = matB->B;

858:   MatEqual(a,c,&flg);
859:   if (flg) {
860:     MatEqual(b,d,&flg);
861:   }
862:   MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
863:   return(0);
864: }

866: PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str)
867: {
869:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
870:   Mat_MPISELL    *b=(Mat_MPISELL*)B->data;

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

888: PetscErrorCode MatSetUp_MPISELL(Mat A)
889: {

893:    MatMPISELLSetPreallocation(A,PETSC_DEFAULT,0,PETSC_DEFAULT,0);
894:   return(0);
895: }


898: extern PetscErrorCode MatConjugate_SeqSELL(Mat);

900: PetscErrorCode MatConjugate_MPISELL(Mat mat)
901: {
902: #if defined(PETSC_USE_COMPLEX)
904:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

907:   MatConjugate_SeqSELL(sell->A);
908:   MatConjugate_SeqSELL(sell->B);
909: #else
911: #endif
912:   return(0);
913: }

915: PetscErrorCode MatRealPart_MPISELL(Mat A)
916: {
917:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

921:   MatRealPart(a->A);
922:   MatRealPart(a->B);
923:   return(0);
924: }

926: PetscErrorCode MatImaginaryPart_MPISELL(Mat A)
927: {
928:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

932:   MatImaginaryPart(a->A);
933:   MatImaginaryPart(a->B);
934:   return(0);
935: }

937: PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values)
938: {
939:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

943:   MatInvertBlockDiagonal(a->A,values);
944:   A->factorerrortype = a->A->factorerrortype;
945:   return(0);
946: }

948: static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx)
949: {
951:   Mat_MPISELL    *sell=(Mat_MPISELL*)x->data;

954:   MatSetRandom(sell->A,rctx);
955:   MatSetRandom(sell->B,rctx);
956:   MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
957:   MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
958:   return(0);
959: }

961: PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A)
962: {

966:   PetscOptionsHead(PetscOptionsObject,"MPISELL options");
967:   PetscOptionsTail();
968:   return(0);
969: }

971: PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a)
972: {
974:   Mat_MPISELL    *msell=(Mat_MPISELL*)Y->data;
975:   Mat_SeqSELL    *sell=(Mat_SeqSELL*)msell->A->data;

978:   if (!Y->preallocated) {
979:     MatMPISELLSetPreallocation(Y,1,NULL,0,NULL);
980:   } else if (!sell->nz) {
981:     PetscInt nonew = sell->nonew;
982:     MatSeqSELLSetPreallocation(msell->A,1,NULL);
983:     sell->nonew = nonew;
984:   }
985:   MatShift_Basic(Y,a);
986:   return(0);
987: }

989: PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool  *missing,PetscInt *d)
990: {
991:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

995:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices");
996:   MatMissingDiagonal(a->A,missing,d);
997:   if (d) {
998:     PetscInt rstart;
999:     MatGetOwnershipRange(A,&rstart,NULL);
1000:     *d += rstart;

1002:   }
1003:   return(0);
1004: }

1006: PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a)
1007: {
1009:   *a = ((Mat_MPISELL*)A->data)->A;
1010:   return(0);
1011: }

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

1161: /* ----------------------------------------------------------------------------------------*/

1163: PetscErrorCode MatStoreValues_MPISELL(Mat mat)
1164: {
1165:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

1169:   MatStoreValues(sell->A);
1170:   MatStoreValues(sell->B);
1171:   return(0);
1172: }

1174: PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
1175: {
1176:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

1180:   MatRetrieveValues(sell->A);
1181:   MatRetrieveValues(sell->B);
1182:   return(0);
1183: }

1185: PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[])
1186: {
1187:   Mat_MPISELL    *b;

1191:   PetscLayoutSetUp(B->rmap);
1192:   PetscLayoutSetUp(B->cmap);
1193:   b = (Mat_MPISELL*)B->data;

1195:   if (!B->preallocated) {
1196:     /* Explicitly create 2 MATSEQSELL matrices. */
1197:     MatCreate(PETSC_COMM_SELF,&b->A);
1198:     MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
1199:     MatSetBlockSizesFromMats(b->A,B,B);
1200:     MatSetType(b->A,MATSEQSELL);
1201:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
1202:     MatCreate(PETSC_COMM_SELF,&b->B);
1203:     MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
1204:     MatSetBlockSizesFromMats(b->B,B,B);
1205:     MatSetType(b->B,MATSEQSELL);
1206:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
1207:   }

1209:   MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);
1210:   MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);
1211:   B->preallocated  = PETSC_TRUE;
1212:   B->was_assembled = PETSC_FALSE;
1213:   /*
1214:     critical for MatAssemblyEnd to work.
1215:     MatAssemblyBegin checks it to set up was_assembled
1216:     and MatAssemblyEnd checks was_assembled to determine whether to build garray
1217:   */
1218:   B->assembled     = PETSC_FALSE;
1219:   return(0);
1220: }

1222: PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat)
1223: {
1224:   Mat            mat;
1225:   Mat_MPISELL    *a,*oldmat=(Mat_MPISELL*)matin->data;

1229:   *newmat = 0;
1230:   MatCreate(PetscObjectComm((PetscObject)matin),&mat);
1231:   MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);
1232:   MatSetBlockSizesFromMats(mat,matin,matin);
1233:   MatSetType(mat,((PetscObject)matin)->type_name);
1234:   a       = (Mat_MPISELL*)mat->data;

1236:   mat->factortype   = matin->factortype;
1237:   mat->assembled    = PETSC_TRUE;
1238:   mat->insertmode   = NOT_SET_VALUES;
1239:   mat->preallocated = PETSC_TRUE;

1241:   a->size         = oldmat->size;
1242:   a->rank         = oldmat->rank;
1243:   a->donotstash   = oldmat->donotstash;
1244:   a->roworiented  = oldmat->roworiented;
1245:   a->rowindices   = 0;
1246:   a->rowvalues    = 0;
1247:   a->getrowactive = PETSC_FALSE;

1249:   PetscLayoutReference(matin->rmap,&mat->rmap);
1250:   PetscLayoutReference(matin->cmap,&mat->cmap);

1252:   if (oldmat->colmap) {
1253: #if defined(PETSC_USE_CTABLE)
1254:     PetscTableCreateCopy(oldmat->colmap,&a->colmap);
1255: #else
1256:     PetscMalloc1(mat->cmap->N,&a->colmap);
1257:     PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));
1258:     PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N);
1259: #endif
1260:   } else a->colmap = 0;
1261:   if (oldmat->garray) {
1262:     PetscInt len;
1263:     len  = oldmat->B->cmap->n;
1264:     PetscMalloc1(len+1,&a->garray);
1265:     PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));
1266:     if (len) { PetscArraycpy(a->garray,oldmat->garray,len); }
1267:   } else a->garray = 0;

1269:   VecDuplicate(oldmat->lvec,&a->lvec);
1270:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);
1271:   VecScatterCopy(oldmat->Mvctx,&a->Mvctx);
1272:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);
1273:   MatDuplicate(oldmat->A,cpvalues,&a->A);
1274:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);
1275:   MatDuplicate(oldmat->B,cpvalues,&a->B);
1276:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);
1277:   PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);
1278:   *newmat = mat;
1279:   return(0);
1280: }

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

1287:    Collective

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

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

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

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

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

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

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

1332:    Example usage:

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

1339: .vb
1340:             1  2  0  |  0  3  0  |  0  4
1341:     Proc0   0  5  6  |  7  0  0  |  8  0
1342:             9  0 10  | 11  0  0  | 12  0
1343:     -------------------------------------
1344:            13  0 14  | 15 16 17  |  0  0
1345:     Proc1   0 18  0  | 19 20 21  |  0  0
1346:             0  0  0  | 22 23  0  | 24  0
1347:     -------------------------------------
1348:     Proc2  25 26 27  |  0  0 28  | 29  0
1349:            30  0  0  | 31 32 33  |  0 34
1350: .ve

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

1354: .vb
1355:       A B C
1356:       D E F
1357:       G H I
1358: .ve

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

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

1367:    The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1368:    submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1369:    corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1370:    Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1371:    part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1372:    matrix, ans [DF] as another SeqSELL matrix.

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

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

1401:    Level: intermediate

1403: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(),
1404:           MATMPISELL, MatGetInfo(), PetscSplitOwnership()
1405: @*/
1406: PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1407: {

1413:   PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1414:   return(0);
1415: }

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

1420:    Collective

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

1446:    Output Parameter:
1447: .  A - the matrix

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

1453:    Notes:
1454:    If the *_rlen parameter is given then the *_rlenmax parameter is ignored

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

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

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

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

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

1477:    The DIAGONAL portion of the local submatrix on any given processor
1478:    is the submatrix corresponding to the rows and columns m,n
1479:    corresponding to the given processor. i.e diagonal matrix on
1480:    process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
1481:    etc. The remaining portion of the local submatrix [m x (N-n)]
1482:    constitute the OFF-DIAGONAL portion. The example below better
1483:    illustrates this concept.

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

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

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

1497:    Options Database Keys:
1498: -  -mat_sell_oneindex - Internally use indexing starting at 1
1499:         rather than 0.  Note that when calling MatSetValues(),
1500:         the user still MUST index entries starting at 0!


1503:    Example usage:

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

1510: .vb
1511:             1  2  0  |  0  3  0  |  0  4
1512:     Proc0   0  5  6  |  7  0  0  |  8  0
1513:             9  0 10  | 11  0  0  | 12  0
1514:     -------------------------------------
1515:            13  0 14  | 15 16 17  |  0  0
1516:     Proc1   0 18  0  | 19 20 21  |  0  0
1517:             0  0  0  | 22 23  0  | 24  0
1518:     -------------------------------------
1519:     Proc2  25 26 27  |  0  0 28  | 29  0
1520:            30  0  0  | 31 32 33  |  0 34
1521: .ve

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

1525: .vb
1526:       A B C
1527:       D E F
1528:       G H I
1529: .ve

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

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

1538:    The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1539:    submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1540:    corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1541:    Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1542:    part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1543:    matrix, ans [DF] as another SeqSELL matrix.

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

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

1572:    Level: intermediate

1574: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(),
1575:           MATMPISELL, MatCreateMPISELLWithArrays()
1576: @*/
1577: 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)
1578: {
1580:   PetscMPIInt    size;

1583:   MatCreate(comm,A);
1584:   MatSetSizes(*A,m,n,M,N);
1585:   MPI_Comm_size(comm,&size);
1586:   if (size > 1) {
1587:     MatSetType(*A,MATMPISELL);
1588:     MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);
1589:   } else {
1590:     MatSetType(*A,MATSEQSELL);
1591:     MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);
1592:   }
1593:   return(0);
1594: }

1596: PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[])
1597: {
1598:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1599:   PetscBool      flg;

1603:   PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);
1604:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input");
1605:   if (Ad)     *Ad     = a->A;
1606:   if (Ao)     *Ao     = a->B;
1607:   if (colmap) *colmap = a->garray;
1608:   return(0);
1609: }

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

1614:     Not Collective

1616:    Input Parameters:
1617: +    A - the matrix
1618: .    scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
1619: -    row, col - index sets of rows and columns to extract (or NULL)

1621:    Output Parameter:
1622: .    A_loc - the local sequential matrix generated

1624:     Level: developer

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

1628: @*/
1629: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc)
1630: {
1631:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1633:   PetscInt       i,start,end,ncols,nzA,nzB,*cmap,imark,*idx;
1634:   IS             isrowa,iscola;
1635:   Mat            *aloc;
1636:   PetscBool      match;

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

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

1685: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1686: {
1688:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1689:   Mat            B;
1690:   Mat_MPIAIJ     *b;

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

1695:   if (reuse == MAT_REUSE_MATRIX) {
1696:     B = *newmat;
1697:   } else {
1698:     MatCreate(PetscObjectComm((PetscObject)A),&B);
1699:     MatSetType(B,MATMPIAIJ);
1700:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1701:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1702:     MatSeqAIJSetPreallocation(B,0,NULL);
1703:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1704:   }
1705:   b    = (Mat_MPIAIJ*) B->data;

1707:   if (reuse == MAT_REUSE_MATRIX) {
1708:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);
1709:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);
1710:   } else {
1711:     MatDestroy(&b->A);
1712:     MatDestroy(&b->B);
1713:     MatDisAssemble_MPISELL(A);
1714:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);
1715:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);
1716:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1717:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1718:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1719:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1720:   }

1722:   if (reuse == MAT_INPLACE_MATRIX) {
1723:     MatHeaderReplace(A,&B);
1724:   } else {
1725:     *newmat = B;
1726:   }
1727:   return(0);
1728: }

1730: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1731: {
1733:   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)A->data;
1734:   Mat            B;
1735:   Mat_MPISELL    *b;

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

1740:   if (reuse == MAT_REUSE_MATRIX) {
1741:     B = *newmat;
1742:   } else {
1743:     MatCreate(PetscObjectComm((PetscObject)A),&B);
1744:     MatSetType(B,MATMPISELL);
1745:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1746:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1747:     MatSeqAIJSetPreallocation(B,0,NULL);
1748:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1749:   }
1750:   b    = (Mat_MPISELL*) B->data;

1752:   if (reuse == MAT_REUSE_MATRIX) {
1753:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);
1754:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);
1755:   } else {
1756:     MatDestroy(&b->A);
1757:     MatDestroy(&b->B);
1758:     MatDisAssemble_MPIAIJ(A);
1759:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);
1760:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);
1761:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1762:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1763:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1764:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1765:   }

1767:   if (reuse == MAT_INPLACE_MATRIX) {
1768:     MatHeaderReplace(A,&B);
1769:   } else {
1770:     *newmat = B;
1771:   }
1772:   return(0);
1773: }

1775: PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1776: {
1777:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
1779:   Vec            bb1=0;

1782:   if (flag == SOR_APPLY_UPPER) {
1783:     (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1784:     return(0);
1785:   }

1787:   if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) {
1788:     VecDuplicate(bb,&bb1);
1789:   }

1791:   if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1792:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1793:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1794:       its--;
1795:     }

1797:     while (its--) {
1798:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1799:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1801:       /* update rhs: bb1 = bb - B*x */
1802:       VecScale(mat->lvec,-1.0);
1803:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

1805:       /* local sweep */
1806:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);
1807:     }
1808:   } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
1809:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1810:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1811:       its--;
1812:     }
1813:     while (its--) {
1814:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1815:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1817:       /* update rhs: bb1 = bb - B*x */
1818:       VecScale(mat->lvec,-1.0);
1819:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

1821:       /* local sweep */
1822:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);
1823:     }
1824:   } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1825:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1826:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1827:       its--;
1828:     }
1829:     while (its--) {
1830:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1831:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1833:       /* update rhs: bb1 = bb - B*x */
1834:       VecScale(mat->lvec,-1.0);
1835:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

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

1842:   VecDestroy(&bb1);

1844:   matin->factorerrortype = mat->A->factorerrortype;
1845:   return(0);
1846: }

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

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

1854:   Level: beginner

1856: .seealso: MatCreateSELL()
1857: M*/
1858: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1859: {
1860:   Mat_MPISELL    *b;
1862:   PetscMPIInt    size;

1865:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1866:   PetscNewLog(B,&b);
1867:   B->data       = (void*)b;
1868:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1869:   B->assembled  = PETSC_FALSE;
1870:   B->insertmode = NOT_SET_VALUES;
1871:   b->size       = size;
1872:   MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);
1873:   /* build cache for off array entries formed */
1874:   MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);

1876:   b->donotstash  = PETSC_FALSE;
1877:   b->colmap      = 0;
1878:   b->garray      = 0;
1879:   b->roworiented = PETSC_TRUE;

1881:   /* stuff used for matrix vector multiply */
1882:   b->lvec  = NULL;
1883:   b->Mvctx = NULL;

1885:   /* stuff for MatGetRow() */
1886:   b->rowindices   = 0;
1887:   b->rowvalues    = 0;
1888:   b->getrowactive = PETSC_FALSE;

1890:   PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);
1891:   PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);
1892:   PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);
1893:   PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);
1894:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);
1895:   PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);
1896:   PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);
1897:   return(0);
1898: }