Actual source code: mpisell.c

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

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

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

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

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

276: PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode)
277: {
278:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
280:   PetscInt       nstash,reallocs;

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

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

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

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

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

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

351: PetscErrorCode MatZeroEntries_MPISELL(Mat A)
352: {
353:   Mat_MPISELL    *l=(Mat_MPISELL*)A->data;

357:   MatZeroEntries(l->A);
358:   MatZeroEntries(l->B);
359:   return(0);
360: }

362: PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy)
363: {
364:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
366:   PetscInt       nt;

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

378: PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx)
379: {
380:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

384:   MatMultDiagonalBlock(a->A,bb,xx);
385:   return(0);
386: }

388: PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
389: {
390:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

394:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
395:   (*a->A->ops->multadd)(a->A,xx,yy,zz);
396:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
397:   (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
398:   return(0);
399: }

401: PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy)
402: {
403:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

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

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

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

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

457: PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
458: {
459:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

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

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

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

489: PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa)
490: {
491:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

495:   MatScale(a->A,aa);
496:   MatScale(a->B,aa);
497:   return(0);
498: }

500: PetscErrorCode MatDestroy_MPISELL(Mat mat)
501: {
502:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

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

525:   PetscObjectChangeTypeName((PetscObject)mat,0);
526:   PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);
527:   PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);
528:   PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL);
529:   PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL);
530:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL);
531:   PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);
532:   return(0);
533: }

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

546:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
547:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
548:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
549:   if (iascii) {
550:     PetscViewerGetFormat(viewer,&format);
551:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
552:       MatInfo   info;
553:       PetscBool inodes;

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

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

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

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

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

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

671: PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer)
672: {
674:   PetscBool      iascii,isdraw,issocket,isbinary;

677:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
678:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
679:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
680:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
681:   if (iascii || isdraw || isbinary || issocket) {
682:     MatView_MPISELL_ASCIIorDraworSocket(mat,viewer);
683:   }
684:   return(0);
685: }

687: PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
688: {
689:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

693:   MatGetSize(sell->B,NULL,nghosts);
694:   if (ghosts) *ghosts = sell->garray;
695:   return(0);
696: }

698: PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info)
699: {
700:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
701:   Mat            A=mat->A,B=mat->B;
703:   PetscReal      isend[5],irecv[5];

706:   info->block_size = 1.0;
707:   MatGetInfo(A,MAT_LOCAL,info);

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

712:   MatGetInfo(B,MAT_LOCAL,info);

714:   isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded;
715:   isend[3] += info->memory;  isend[4] += info->mallocs;
716:   if (flag == MAT_LOCAL) {
717:     info->nz_used      = isend[0];
718:     info->nz_allocated = isend[1];
719:     info->nz_unneeded  = isend[2];
720:     info->memory       = isend[3];
721:     info->mallocs      = isend[4];
722:   } else if (flag == MAT_GLOBAL_MAX) {
723:     MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin));

725:     info->nz_used      = irecv[0];
726:     info->nz_allocated = irecv[1];
727:     info->nz_unneeded  = irecv[2];
728:     info->memory       = irecv[3];
729:     info->mallocs      = irecv[4];
730:   } else if (flag == MAT_GLOBAL_SUM) {
731:     MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)matin));

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

745: PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg)
746: {
747:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

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

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


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

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

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

840: PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
841: {
842:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

846:   MatSetUnfactored(a->A);
847:   return(0);
848: }

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

858:   a = matA->A; b = matA->B;
859:   c = matB->A; d = matB->B;

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

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

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

891: PetscErrorCode MatSetUp_MPISELL(Mat A)
892: {

896:    MatMPISELLSetPreallocation(A,PETSC_DEFAULT,0,PETSC_DEFAULT,0);
897:   return(0);
898: }


901: extern PetscErrorCode MatConjugate_SeqSELL(Mat);

903: PetscErrorCode MatConjugate_MPISELL(Mat mat)
904: {
905: #if defined(PETSC_USE_COMPLEX)
907:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

910:   MatConjugate_SeqSELL(sell->A);
911:   MatConjugate_SeqSELL(sell->B);
912: #else
914: #endif
915:   return(0);
916: }

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

924:   MatRealPart(a->A);
925:   MatRealPart(a->B);
926:   return(0);
927: }

929: PetscErrorCode MatImaginaryPart_MPISELL(Mat A)
930: {
931:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

935:   MatImaginaryPart(a->A);
936:   MatImaginaryPart(a->B);
937:   return(0);
938: }

940: PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values)
941: {
942:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

946:   MatInvertBlockDiagonal(a->A,values);
947:   A->factorerrortype = a->A->factorerrortype;
948:   return(0);
949: }

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

957:   MatSetRandom(sell->A,rctx);
958:   MatSetRandom(sell->B,rctx);
959:   MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
960:   MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
961:   return(0);
962: }

964: PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A)
965: {

969:   PetscOptionsHead(PetscOptionsObject,"MPISELL options");
970:   PetscOptionsTail();
971:   return(0);
972: }

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

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

992: PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool  *missing,PetscInt *d)
993: {
994:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

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

1005:   }
1006:   return(0);
1007: }

1009: PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a)
1010: {
1012:   *a = ((Mat_MPISELL*)A->data)->A;
1013:   return(0);
1014: }

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

1164: /* ----------------------------------------------------------------------------------------*/

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

1172:   MatStoreValues(sell->A);
1173:   MatStoreValues(sell->B);
1174:   return(0);
1175: }

1177: PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
1178: {
1179:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

1183:   MatRetrieveValues(sell->A);
1184:   MatRetrieveValues(sell->B);
1185:   return(0);
1186: }

1188: PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[])
1189: {
1190:   Mat_MPISELL    *b;

1194:   PetscLayoutSetUp(B->rmap);
1195:   PetscLayoutSetUp(B->cmap);
1196:   b = (Mat_MPISELL*)B->data;

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

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

1225: PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat)
1226: {
1227:   Mat            mat;
1228:   Mat_MPISELL    *a,*oldmat=(Mat_MPISELL*)matin->data;

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

1239:   mat->factortype   = matin->factortype;
1240:   mat->assembled    = PETSC_TRUE;
1241:   mat->insertmode   = NOT_SET_VALUES;
1242:   mat->preallocated = PETSC_TRUE;

1244:   a->size         = oldmat->size;
1245:   a->rank         = oldmat->rank;
1246:   a->donotstash   = oldmat->donotstash;
1247:   a->roworiented  = oldmat->roworiented;
1248:   a->rowindices   = 0;
1249:   a->rowvalues    = 0;
1250:   a->getrowactive = PETSC_FALSE;

1252:   PetscLayoutReference(matin->rmap,&mat->rmap);
1253:   PetscLayoutReference(matin->cmap,&mat->cmap);

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

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

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

1290:    Collective on MPI_Comm

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

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

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

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

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

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

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

1335:    Example usage:

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

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

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

1357: .vb
1358:       A B C
1359:       D E F
1360:       G H I
1361: .ve

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

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

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

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

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

1404:    Level: intermediate

1406: .keywords: matrix, sell, sparse, parallel

1408: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(),
1409:           MATMPISELL, MatGetInfo(), PetscSplitOwnership()
1410: @*/
1411: PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1412: {

1418:   PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1419:   return(0);
1420: }

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

1425:    Collective on MPI_Comm

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

1451:    Output Parameter:
1452: .  A - the matrix

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

1458:    Notes:
1459:    If the *_rlen parameter is given then the *_rlenmax parameter is ignored

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

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

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

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

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

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

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

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

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

1502:    Options Database Keys:
1503: -  -mat_sell_oneindex - Internally use indexing starting at 1
1504:         rather than 0.  Note that when calling MatSetValues(),
1505:         the user still MUST index entries starting at 0!


1508:    Example usage:

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

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

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

1530: .vb
1531:       A B C
1532:       D E F
1533:       G H I
1534: .ve

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

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

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

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

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

1577:    Level: intermediate

1579: .keywords: matrix, sell, sparse, parallel

1581: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(),
1582:           MATMPISELL, MatCreateMPISELLWithArrays()
1583: @*/
1584: 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)
1585: {
1587:   PetscMPIInt    size;

1590:   MatCreate(comm,A);
1591:   MatSetSizes(*A,m,n,M,N);
1592:   MPI_Comm_size(comm,&size);
1593:   if (size > 1) {
1594:     MatSetType(*A,MATMPISELL);
1595:     MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);
1596:   } else {
1597:     MatSetType(*A,MATSEQSELL);
1598:     MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);
1599:   }
1600:   return(0);
1601: }

1603: PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[])
1604: {
1605:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1606:   PetscBool      flg;

1610:   PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);
1611:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input");
1612:   if (Ad)     *Ad     = a->A;
1613:   if (Ao)     *Ao     = a->B;
1614:   if (colmap) *colmap = a->garray;
1615:   return(0);
1616: }

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

1621:     Not Collective

1623:    Input Parameters:
1624: +    A - the matrix
1625: .    scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
1626: -    row, col - index sets of rows and columns to extract (or NULL)

1628:    Output Parameter:
1629: .    A_loc - the local sequential matrix generated

1631:     Level: developer

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

1635: @*/
1636: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc)
1637: {
1638:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1640:   PetscInt       i,start,end,ncols,nzA,nzB,*cmap,imark,*idx;
1641:   IS             isrowa,iscola;
1642:   Mat            *aloc;
1643:   PetscBool      match;

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

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

1692: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1693: {
1695:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1696:   Mat            B;
1697:   Mat_MPIAIJ     *b;

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

1702:   if (reuse == MAT_REUSE_MATRIX) {
1703:     B = *newmat;
1704:   } else {
1705:     MatCreate(PetscObjectComm((PetscObject)A),&B);
1706:     MatSetType(B,MATMPIAIJ);
1707:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1708:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1709:     MatSeqAIJSetPreallocation(B,0,NULL);
1710:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1711:   }
1712:   b    = (Mat_MPIAIJ*) B->data;

1714:   if (reuse == MAT_REUSE_MATRIX) {
1715:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);
1716:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);
1717:   } else {
1718:     MatDestroy(&b->A);
1719:     MatDestroy(&b->B);
1720:     MatDisAssemble_MPISELL(A);
1721:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);
1722:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);
1723:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1724:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1725:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1726:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1727:   }

1729:   if (reuse == MAT_INPLACE_MATRIX) {
1730:     MatHeaderReplace(A,&B);
1731:   } else {
1732:     *newmat = B;
1733:   }
1734:   return(0);
1735: }

1737: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1738: {
1740:   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)A->data;
1741:   Mat            B;
1742:   Mat_MPISELL    *b;

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

1747:   if (reuse == MAT_REUSE_MATRIX) {
1748:     B = *newmat;
1749:   } else {
1750:     MatCreate(PetscObjectComm((PetscObject)A),&B);
1751:     MatSetType(B,MATMPISELL);
1752:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1753:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1754:     MatSeqAIJSetPreallocation(B,0,NULL);
1755:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1756:   }
1757:   b    = (Mat_MPISELL*) B->data;

1759:   if (reuse == MAT_REUSE_MATRIX) {
1760:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);
1761:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);
1762:   } else {
1763:     MatDestroy(&b->A);
1764:     MatDestroy(&b->B);
1765:     MatDisAssemble_MPIAIJ(A);
1766:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);
1767:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);
1768:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1769:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1770:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1771:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1772:   }

1774:   if (reuse == MAT_INPLACE_MATRIX) {
1775:     MatHeaderReplace(A,&B);
1776:   } else {
1777:     *newmat = B;
1778:   }
1779:   return(0);
1780: }

1782: PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1783: {
1784:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
1786:   Vec            bb1=0;

1789:   if (flag == SOR_APPLY_UPPER) {
1790:     (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1791:     return(0);
1792:   }

1794:   if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) {
1795:     VecDuplicate(bb,&bb1);
1796:   }

1798:   if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1799:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1800:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1801:       its--;
1802:     }

1804:     while (its--) {
1805:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1806:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

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

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

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

1828:       /* local sweep */
1829:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);
1830:     }
1831:   } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1832:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1833:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1834:       its--;
1835:     }
1836:     while (its--) {
1837:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1838:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1840:       /* update rhs: bb1 = bb - B*x */
1841:       VecScale(mat->lvec,-1.0);
1842:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

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

1849:   VecDestroy(&bb1);

1851:   matin->factorerrortype = mat->A->factorerrortype;
1852:   return(0);
1853: }

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

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

1861:   Level: beginner

1863: .seealso: MatCreateSELL()
1864: M*/
1865: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1866: {
1867:   Mat_MPISELL    *b;
1869:   PetscMPIInt    size;

1872:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1873:   PetscNewLog(B,&b);
1874:   B->data       = (void*)b;
1875:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1876:   B->assembled  = PETSC_FALSE;
1877:   B->insertmode = NOT_SET_VALUES;
1878:   b->size       = size;
1879:   MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);
1880:   /* build cache for off array entries formed */
1881:   MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);

1883:   b->donotstash  = PETSC_FALSE;
1884:   b->colmap      = 0;
1885:   b->garray      = 0;
1886:   b->roworiented = PETSC_TRUE;

1888:   /* stuff used for matrix vector multiply */
1889:   b->lvec  = NULL;
1890:   b->Mvctx = NULL;

1892:   /* stuff for MatGetRow() */
1893:   b->rowindices   = 0;
1894:   b->rowvalues    = 0;
1895:   b->getrowactive = PETSC_FALSE;

1897:   PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);
1898:   PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);
1899:   PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);
1900:   PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);
1901:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);
1902:   PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);
1903:   PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);
1904:   return(0);
1905: }