Actual source code: mpisell.c

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

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

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

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

 20:   Level: beginner

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

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

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

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

 50: #if defined(PETSC_USE_CTABLE)
 51:   PetscTableCreate(n,mat->cmap->N+1,&sell->colmap);
 52:   for (i=0; i<n; i++) {
 53:     PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES);
 54:   }
 55: #else
 56:   PetscCalloc1(mat->cmap->N+1,&sell->colmap);
 57:   PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt));
 58:   for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1;
 59: #endif
 60:   return 0;
 61: }

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

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

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

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

148:   for (i=0; i<m; i++) {
149:     if (im[i] < 0) continue;
151:     if (im[i] >= rstart && im[i] < rend) {
152:       row      = im[i] - rstart;
153:       lastcol1 = -1;
154:       shift1   = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */
155:       cp1      = a->colidx+shift1;
156:       vp1      = a->val+shift1;
157:       nrow1    = a->rlen[row];
158:       low1     = 0;
159:       high1    = nrow1;
160:       lastcol2 = -1;
161:       shift2   = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
162:       cp2      = b->colidx+shift2;
163:       vp2      = b->val+shift2;
164:       nrow2    = b->rlen[row];
165:       low2     = 0;
166:       high2    = nrow2;

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

220: PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
221: {
222:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
223:   PetscInt       i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend;
224:   PetscInt       cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;

226:   for (i=0; i<m; i++) {
227:     if (idxm[i] < 0) continue; /* negative row */
229:     if (idxm[i] >= rstart && idxm[i] < rend) {
230:       row = idxm[i] - rstart;
231:       for (j=0; j<n; j++) {
232:         if (idxn[j] < 0) continue; /* negative column */
234:         if (idxn[j] >= cstart && idxn[j] < cend) {
235:           col  = idxn[j] - cstart;
236:           MatGetValues(sell->A,1,&row,1,&col,v+i*n+j);
237:         } else {
238:           if (!sell->colmap) {
239:             MatCreateColmap_MPISELL_Private(mat);
240:           }
241: #if defined(PETSC_USE_CTABLE)
242:           PetscTableFind(sell->colmap,idxn[j]+1,&col);
243:           col--;
244: #else
245:           col = sell->colmap[idxn[j]] - 1;
246: #endif
247:           if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0;
248:           else {
249:             MatGetValues(sell->B,1,&row,1,&col,v+i*n+j);
250:           }
251:         }
252:       }
253:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported");
254:   }
255:   return 0;
256: }

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

260: PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode)
261: {
262:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
263:   PetscInt       nstash,reallocs;

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

267:   MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
268:   MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);
269:   PetscInfo(sell->A,"Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n",nstash,reallocs);
270:   return 0;
271: }

273: PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode)
274: {
275:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
276:   PetscMPIInt    n;
277:   PetscInt       i,flg;
278:   PetscInt       *row,*col;
279:   PetscScalar    *val;
280:   PetscBool      other_disassembled;
281:   /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */
282:   if (!sell->donotstash && !mat->nooffprocentries) {
283:     while (1) {
284:       MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);
285:       if (!flg) break;

287:       for (i=0; i<n; i++) { /* assemble one by one */
288:         MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode);
289:       }
290:     }
291:     MatStashScatterEnd_Private(&mat->stash);
292:   }
293:   MatAssemblyBegin(sell->A,mode);
294:   MatAssemblyEnd(sell->A,mode);

296:   /*
297:      determine if any processor has disassembled, if so we must
298:      also disassemble ourselfs, in order that we may reassemble.
299:   */
300:   /*
301:      if nonzero structure of submatrix B cannot change then we know that
302:      no processor disassembled thus we can skip this stuff
303:   */
304:   if (!((Mat_SeqSELL*)sell->B->data)->nonew) {
305:     MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat));
307:   }
308:   if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) {
309:     MatSetUpMultiply_MPISELL(mat);
310:   }
311:   /*
312:   MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE);
313:   */
314:   MatAssemblyBegin(sell->B,mode);
315:   MatAssemblyEnd(sell->B,mode);
316:   PetscFree2(sell->rowvalues,sell->rowindices);
317:   sell->rowvalues = NULL;
318:   VecDestroy(&sell->diag);

320:   /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
321:   if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) {
322:     PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate;
323:     MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));
324:   }
325:   return 0;
326: }

328: PetscErrorCode MatZeroEntries_MPISELL(Mat A)
329: {
330:   Mat_MPISELL    *l=(Mat_MPISELL*)A->data;

332:   MatZeroEntries(l->A);
333:   MatZeroEntries(l->B);
334:   return 0;
335: }

337: PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy)
338: {
339:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
340:   PetscInt       nt;

342:   VecGetLocalSize(xx,&nt);
344:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
345:   (*a->A->ops->mult)(a->A,xx,yy);
346:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
347:   (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);
348:   return 0;
349: }

351: PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx)
352: {
353:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

355:   MatMultDiagonalBlock(a->A,bb,xx);
356:   return 0;
357: }

359: PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
360: {
361:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

363:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
364:   (*a->A->ops->multadd)(a->A,xx,yy,zz);
365:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
366:   (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
367:   return 0;
368: }

370: PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy)
371: {
372:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

374:   /* do nondiagonal part */
375:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
376:   /* do local part */
377:   (*a->A->ops->multtranspose)(a->A,xx,yy);
378:   /* add partial results together */
379:   VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
380:   VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
381:   return 0;
382: }

384: PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f)
385: {
386:   MPI_Comm       comm;
387:   Mat_MPISELL    *Asell=(Mat_MPISELL*)Amat->data,*Bsell;
388:   Mat            Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs;
389:   IS             Me,Notme;
390:   PetscInt       M,N,first,last,*notme,i;
391:   PetscMPIInt    size;

393:   /* Easy test: symmetric diagonal block */
394:   Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A;
395:   MatIsTranspose(Adia,Bdia,tol,f);
396:   if (!*f) return 0;
397:   PetscObjectGetComm((PetscObject)Amat,&comm);
398:   MPI_Comm_size(comm,&size);
399:   if (size == 1) return 0;

401:   /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */
402:   MatGetSize(Amat,&M,&N);
403:   MatGetOwnershipRange(Amat,&first,&last);
404:   PetscMalloc1(N-last+first,&notme);
405:   for (i=0; i<first; i++) notme[i] = i;
406:   for (i=last; i<M; i++) notme[i-last+first] = i;
407:   ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme);
408:   ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me);
409:   MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs);
410:   Aoff = Aoffs[0];
411:   MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs);
412:   Boff = Boffs[0];
413:   MatIsTranspose(Aoff,Boff,tol,f);
414:   MatDestroyMatrices(1,&Aoffs);
415:   MatDestroyMatrices(1,&Boffs);
416:   ISDestroy(&Me);
417:   ISDestroy(&Notme);
418:   PetscFree(notme);
419:   return 0;
420: }

422: PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
423: {
424:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

426:   /* do nondiagonal part */
427:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
428:   /* do local part */
429:   (*a->A->ops->multtransposeadd)(a->A,xx,yy,zz);
430:   /* add partial results together */
431:   VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
432:   VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
433:   return 0;
434: }

436: /*
437:   This only works correctly for square matrices where the subblock A->A is the
438:    diagonal block
439: */
440: PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v)
441: {
442:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

446:   MatGetDiagonal(a->A,v);
447:   return 0;
448: }

450: PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa)
451: {
452:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

454:   MatScale(a->A,aa);
455:   MatScale(a->B,aa);
456:   return 0;
457: }

459: PetscErrorCode MatDestroy_MPISELL(Mat mat)
460: {
461:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

463: #if defined(PETSC_USE_LOG)
464:   PetscLogObjectState((PetscObject)mat,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT,mat->rmap->N,mat->cmap->N);
465: #endif
466:   MatStashDestroy_Private(&mat->stash);
467:   VecDestroy(&sell->diag);
468:   MatDestroy(&sell->A);
469:   MatDestroy(&sell->B);
470: #if defined(PETSC_USE_CTABLE)
471:   PetscTableDestroy(&sell->colmap);
472: #else
473:   PetscFree(sell->colmap);
474: #endif
475:   PetscFree(sell->garray);
476:   VecDestroy(&sell->lvec);
477:   VecScatterDestroy(&sell->Mvctx);
478:   PetscFree2(sell->rowvalues,sell->rowindices);
479:   PetscFree(sell->ld);
480:   PetscFree(mat->data);

482:   PetscObjectChangeTypeName((PetscObject)mat,NULL);
483:   PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);
484:   PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);
485:   PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL);
486:   PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL);
487:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL);
488:   PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);
489:   return 0;
490: }

492: #include <petscdraw.h>
493: PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer)
494: {
495:   Mat_MPISELL       *sell=(Mat_MPISELL*)mat->data;
496:   PetscErrorCode    ierr;
497:   PetscMPIInt       rank=sell->rank,size=sell->size;
498:   PetscBool         isdraw,iascii,isbinary;
499:   PetscViewer       sviewer;
500:   PetscViewerFormat format;

502:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
503:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
504:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
505:   if (iascii) {
506:     PetscViewerGetFormat(viewer,&format);
507:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
508:       MatInfo   info;
509:       PetscInt *inodes;

511:       MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank);
512:       MatGetInfo(mat,MAT_LOCAL,&info);
513:       MatInodeGetInodeSizes(sell->A,NULL,&inodes,NULL);
514:       PetscViewerASCIIPushSynchronized(viewer);
515:       if (!inodes) {
516:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", not using I-node routines\n",
517:                                                   rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
518:       } else {
519:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", using I-node routines\n",
520:                                                   rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
521:       }
522:       MatGetInfo(sell->A,MAT_LOCAL,&info);
523:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used);
524:       MatGetInfo(sell->B,MAT_LOCAL,&info);
525:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used);
526:       PetscViewerFlush(viewer);
527:       PetscViewerASCIIPopSynchronized(viewer);
528:       PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n");
529:       VecScatterView(sell->Mvctx,viewer);
530:       return 0;
531:     } else if (format == PETSC_VIEWER_ASCII_INFO) {
532:       PetscInt inodecount,inodelimit,*inodes;
533:       MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit);
534:       if (inodes) {
535:         PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n",inodecount,inodelimit);
536:       } else {
537:         PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n");
538:       }
539:       return 0;
540:     } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
541:       return 0;
542:     }
543:   } else if (isbinary) {
544:     if (size == 1) {
545:       PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name);
546:       MatView(sell->A,viewer);
547:     } else {
548:       /* MatView_MPISELL_Binary(mat,viewer); */
549:     }
550:     return 0;
551:   } else if (isdraw) {
552:     PetscDraw draw;
553:     PetscBool isnull;
554:     PetscViewerDrawGetDraw(viewer,0,&draw);
555:     PetscDrawIsNull(draw,&isnull);
556:     if (isnull) return 0;
557:   }

559:   {
560:     /* assemble the entire matrix onto first processor. */
561:     Mat         A;
562:     Mat_SeqSELL *Aloc;
563:     PetscInt    M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j;
564:     MatScalar   *aval;
565:     PetscBool   isnonzero;

567:     MatCreate(PetscObjectComm((PetscObject)mat),&A);
568:     if (rank == 0) {
569:       MatSetSizes(A,M,N,M,N);
570:     } else {
571:       MatSetSizes(A,0,0,M,N);
572:     }
573:     /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */
574:     MatSetType(A,MATMPISELL);
575:     MatMPISELLSetPreallocation(A,0,NULL,0,NULL);
576:     MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
577:     PetscLogObjectParent((PetscObject)mat,(PetscObject)A);

579:     /* copy over the A part */
580:     Aloc = (Mat_SeqSELL*)sell->A->data;
581:     acolidx = Aloc->colidx; aval = Aloc->val;
582:     for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */
583:       for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
584:         isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
585:         if (isnonzero) { /* check the mask bit */
586:           row  = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */
587:           col  = *acolidx + mat->rmap->rstart;
588:           MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
589:         }
590:         aval++; acolidx++;
591:       }
592:     }

594:     /* copy over the B part */
595:     Aloc = (Mat_SeqSELL*)sell->B->data;
596:     acolidx = Aloc->colidx; aval = Aloc->val;
597:     for (i=0; i<Aloc->totalslices; i++) {
598:       for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
599:         isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
600:         if (isnonzero) {
601:           row  = (i<<3)+(j&0x07) + mat->rmap->rstart;
602:           col  = sell->garray[*acolidx];
603:           MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
604:         }
605:         aval++; acolidx++;
606:       }
607:     }

609:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
610:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
611:     /*
612:        Everyone has to call to draw the matrix since the graphics waits are
613:        synchronized across all processors that share the PetscDraw object
614:     */
615:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
616:     if (rank == 0) {
617:       PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name);
618:       MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer);
619:     }
620:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
621:     PetscViewerFlush(viewer);
622:     MatDestroy(&A);
623:   }
624:   return 0;
625: }

627: PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer)
628: {
629:   PetscBool      iascii,isdraw,issocket,isbinary;

631:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
632:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
633:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
634:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
635:   if (iascii || isdraw || isbinary || issocket) {
636:     MatView_MPISELL_ASCIIorDraworSocket(mat,viewer);
637:   }
638:   return 0;
639: }

641: PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
642: {
643:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

645:   MatGetSize(sell->B,NULL,nghosts);
646:   if (ghosts) *ghosts = sell->garray;
647:   return 0;
648: }

650: PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info)
651: {
652:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
653:   Mat            A=mat->A,B=mat->B;
654:   PetscLogDouble isend[5],irecv[5];

656:   info->block_size = 1.0;
657:   MatGetInfo(A,MAT_LOCAL,info);

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

662:   MatGetInfo(B,MAT_LOCAL,info);

664:   isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded;
665:   isend[3] += info->memory;  isend[4] += info->mallocs;
666:   if (flag == MAT_LOCAL) {
667:     info->nz_used      = isend[0];
668:     info->nz_allocated = isend[1];
669:     info->nz_unneeded  = isend[2];
670:     info->memory       = isend[3];
671:     info->mallocs      = isend[4];
672:   } else if (flag == MAT_GLOBAL_MAX) {
673:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin));

675:     info->nz_used      = irecv[0];
676:     info->nz_allocated = irecv[1];
677:     info->nz_unneeded  = irecv[2];
678:     info->memory       = irecv[3];
679:     info->mallocs      = irecv[4];
680:   } else if (flag == MAT_GLOBAL_SUM) {
681:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin));

683:     info->nz_used      = irecv[0];
684:     info->nz_allocated = irecv[1];
685:     info->nz_unneeded  = irecv[2];
686:     info->memory       = irecv[3];
687:     info->mallocs      = irecv[4];
688:   }
689:   info->fill_ratio_given  = 0; /* no parallel LU/ILU/Cholesky */
690:   info->fill_ratio_needed = 0;
691:   info->factor_mallocs    = 0;
692:   return 0;
693: }

695: PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg)
696: {
697:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

699:   switch (op) {
700:   case MAT_NEW_NONZERO_LOCATIONS:
701:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
702:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
703:   case MAT_KEEP_NONZERO_PATTERN:
704:   case MAT_NEW_NONZERO_LOCATION_ERR:
705:   case MAT_USE_INODES:
706:   case MAT_IGNORE_ZERO_ENTRIES:
707:     MatCheckPreallocated(A,1);
708:     MatSetOption(a->A,op,flg);
709:     MatSetOption(a->B,op,flg);
710:     break;
711:   case MAT_ROW_ORIENTED:
712:     MatCheckPreallocated(A,1);
713:     a->roworiented = flg;

715:     MatSetOption(a->A,op,flg);
716:     MatSetOption(a->B,op,flg);
717:     break;
718:   case MAT_FORCE_DIAGONAL_ENTRIES:
719:   case MAT_SORTED_FULL:
720:     PetscInfo(A,"Option %s ignored\n",MatOptions[op]);
721:     break;
722:   case MAT_IGNORE_OFF_PROC_ENTRIES:
723:     a->donotstash = flg;
724:     break;
725:   case MAT_SPD:
726:     A->spd_set = PETSC_TRUE;
727:     A->spd     = flg;
728:     if (flg) {
729:       A->symmetric                  = PETSC_TRUE;
730:       A->structurally_symmetric     = PETSC_TRUE;
731:       A->symmetric_set              = PETSC_TRUE;
732:       A->structurally_symmetric_set = PETSC_TRUE;
733:     }
734:     break;
735:   case MAT_SYMMETRIC:
736:     MatCheckPreallocated(A,1);
737:     MatSetOption(a->A,op,flg);
738:     break;
739:   case MAT_STRUCTURALLY_SYMMETRIC:
740:     MatCheckPreallocated(A,1);
741:     MatSetOption(a->A,op,flg);
742:     break;
743:   case MAT_HERMITIAN:
744:     MatCheckPreallocated(A,1);
745:     MatSetOption(a->A,op,flg);
746:     break;
747:   case MAT_SYMMETRY_ETERNAL:
748:     MatCheckPreallocated(A,1);
749:     MatSetOption(a->A,op,flg);
750:     break;
751:   default:
752:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
753:   }
754:   return 0;
755: }

757: PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr)
758: {
759:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;
760:   Mat            a=sell->A,b=sell->B;
761:   PetscInt       s1,s2,s3;

763:   MatGetLocalSize(mat,&s2,&s3);
764:   if (rr) {
765:     VecGetLocalSize(rr,&s1);
767:     /* Overlap communication with computation. */
768:     VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
769:   }
770:   if (ll) {
771:     VecGetLocalSize(ll,&s1);
773:     (*b->ops->diagonalscale)(b,ll,NULL);
774:   }
775:   /* scale  the diagonal block */
776:   (*a->ops->diagonalscale)(a,ll,rr);

778:   if (rr) {
779:     /* Do a scatter end and then right scale the off-diagonal block */
780:     VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
781:     (*b->ops->diagonalscale)(b,NULL,sell->lvec);
782:   }
783:   return 0;
784: }

786: PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
787: {
788:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

790:   MatSetUnfactored(a->A);
791:   return 0;
792: }

794: PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool  *flag)
795: {
796:   Mat_MPISELL    *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data;
797:   Mat            a,b,c,d;
798:   PetscBool      flg;

800:   a = matA->A; b = matA->B;
801:   c = matB->A; d = matB->B;

803:   MatEqual(a,c,&flg);
804:   if (flg) {
805:     MatEqual(b,d,&flg);
806:   }
807:   MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
808:   return 0;
809: }

811: PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str)
812: {
813:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
814:   Mat_MPISELL    *b=(Mat_MPISELL*)B->data;

816:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
817:   if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
818:     /* because of the column compression in the off-processor part of the matrix a->B,
819:        the number of columns in a->B and b->B may be different, hence we cannot call
820:        the MatCopy() directly on the two parts. If need be, we can provide a more
821:        efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices
822:        then copying the submatrices */
823:     MatCopy_Basic(A,B,str);
824:   } else {
825:     MatCopy(a->A,b->A,str);
826:     MatCopy(a->B,b->B,str);
827:   }
828:   return 0;
829: }

831: PetscErrorCode MatSetUp_MPISELL(Mat A)
832: {
833:   MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
834:   return 0;
835: }

837: extern PetscErrorCode MatConjugate_SeqSELL(Mat);

839: PetscErrorCode MatConjugate_MPISELL(Mat mat)
840: {
841:   if (PetscDefined(USE_COMPLEX)) {
842:     Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;

844:     MatConjugate_SeqSELL(sell->A);
845:     MatConjugate_SeqSELL(sell->B);
846:   }
847:   return 0;
848: }

850: PetscErrorCode MatRealPart_MPISELL(Mat A)
851: {
852:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

854:   MatRealPart(a->A);
855:   MatRealPart(a->B);
856:   return 0;
857: }

859: PetscErrorCode MatImaginaryPart_MPISELL(Mat A)
860: {
861:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

863:   MatImaginaryPart(a->A);
864:   MatImaginaryPart(a->B);
865:   return 0;
866: }

868: PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values)
869: {
870:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

872:   MatInvertBlockDiagonal(a->A,values);
873:   A->factorerrortype = a->A->factorerrortype;
874:   return 0;
875: }

877: static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx)
878: {
879:   Mat_MPISELL    *sell=(Mat_MPISELL*)x->data;

881:   MatSetRandom(sell->A,rctx);
882:   MatSetRandom(sell->B,rctx);
883:   MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
884:   MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
885:   return 0;
886: }

888: PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A)
889: {
890:   PetscOptionsHead(PetscOptionsObject,"MPISELL options");
891:   PetscOptionsTail();
892:   return 0;
893: }

895: PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a)
896: {
897:   Mat_MPISELL    *msell=(Mat_MPISELL*)Y->data;
898:   Mat_SeqSELL    *sell=(Mat_SeqSELL*)msell->A->data;

900:   if (!Y->preallocated) {
901:     MatMPISELLSetPreallocation(Y,1,NULL,0,NULL);
902:   } else if (!sell->nz) {
903:     PetscInt nonew = sell->nonew;
904:     MatSeqSELLSetPreallocation(msell->A,1,NULL);
905:     sell->nonew = nonew;
906:   }
907:   MatShift_Basic(Y,a);
908:   return 0;
909: }

911: PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool  *missing,PetscInt *d)
912: {
913:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;

916:   MatMissingDiagonal(a->A,missing,d);
917:   if (d) {
918:     PetscInt rstart;
919:     MatGetOwnershipRange(A,&rstart,NULL);
920:     *d += rstart;

922:   }
923:   return 0;
924: }

926: PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a)
927: {
928:   *a = ((Mat_MPISELL*)A->data)->A;
929:   return 0;
930: }

932: /* -------------------------------------------------------------------*/
933: static struct _MatOps MatOps_Values = {MatSetValues_MPISELL,
934:                                        NULL,
935:                                        NULL,
936:                                        MatMult_MPISELL,
937:                                 /* 4*/ MatMultAdd_MPISELL,
938:                                        MatMultTranspose_MPISELL,
939:                                        MatMultTransposeAdd_MPISELL,
940:                                        NULL,
941:                                        NULL,
942:                                        NULL,
943:                                 /*10*/ NULL,
944:                                        NULL,
945:                                        NULL,
946:                                        MatSOR_MPISELL,
947:                                        NULL,
948:                                 /*15*/ MatGetInfo_MPISELL,
949:                                        MatEqual_MPISELL,
950:                                        MatGetDiagonal_MPISELL,
951:                                        MatDiagonalScale_MPISELL,
952:                                        NULL,
953:                                 /*20*/ MatAssemblyBegin_MPISELL,
954:                                        MatAssemblyEnd_MPISELL,
955:                                        MatSetOption_MPISELL,
956:                                        MatZeroEntries_MPISELL,
957:                                 /*24*/ NULL,
958:                                        NULL,
959:                                        NULL,
960:                                        NULL,
961:                                        NULL,
962:                                 /*29*/ MatSetUp_MPISELL,
963:                                        NULL,
964:                                        NULL,
965:                                        MatGetDiagonalBlock_MPISELL,
966:                                        NULL,
967:                                 /*34*/ MatDuplicate_MPISELL,
968:                                        NULL,
969:                                        NULL,
970:                                        NULL,
971:                                        NULL,
972:                                 /*39*/ NULL,
973:                                        NULL,
974:                                        NULL,
975:                                        MatGetValues_MPISELL,
976:                                        MatCopy_MPISELL,
977:                                 /*44*/ NULL,
978:                                        MatScale_MPISELL,
979:                                        MatShift_MPISELL,
980:                                        MatDiagonalSet_MPISELL,
981:                                        NULL,
982:                                 /*49*/ MatSetRandom_MPISELL,
983:                                        NULL,
984:                                        NULL,
985:                                        NULL,
986:                                        NULL,
987:                                 /*54*/ MatFDColoringCreate_MPIXAIJ,
988:                                        NULL,
989:                                        MatSetUnfactored_MPISELL,
990:                                        NULL,
991:                                        NULL,
992:                                 /*59*/ NULL,
993:                                        MatDestroy_MPISELL,
994:                                        MatView_MPISELL,
995:                                        NULL,
996:                                        NULL,
997:                                 /*64*/ NULL,
998:                                        NULL,
999:                                        NULL,
1000:                                        NULL,
1001:                                        NULL,
1002:                                 /*69*/ NULL,
1003:                                        NULL,
1004:                                        NULL,
1005:                                        NULL,
1006:                                        NULL,
1007:                                        NULL,
1008:                                 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */
1009:                                        MatSetFromOptions_MPISELL,
1010:                                        NULL,
1011:                                        NULL,
1012:                                        NULL,
1013:                                 /*80*/ NULL,
1014:                                        NULL,
1015:                                        NULL,
1016:                                 /*83*/ NULL,
1017:                                        NULL,
1018:                                        NULL,
1019:                                        NULL,
1020:                                        NULL,
1021:                                        NULL,
1022:                                 /*89*/ NULL,
1023:                                        NULL,
1024:                                        NULL,
1025:                                        NULL,
1026:                                        NULL,
1027:                                 /*94*/ NULL,
1028:                                        NULL,
1029:                                        NULL,
1030:                                        NULL,
1031:                                        NULL,
1032:                                 /*99*/ NULL,
1033:                                        NULL,
1034:                                        NULL,
1035:                                        MatConjugate_MPISELL,
1036:                                        NULL,
1037:                                 /*104*/NULL,
1038:                                        MatRealPart_MPISELL,
1039:                                        MatImaginaryPart_MPISELL,
1040:                                        NULL,
1041:                                        NULL,
1042:                                 /*109*/NULL,
1043:                                        NULL,
1044:                                        NULL,
1045:                                        NULL,
1046:                                        MatMissingDiagonal_MPISELL,
1047:                                 /*114*/NULL,
1048:                                        NULL,
1049:                                        MatGetGhosts_MPISELL,
1050:                                        NULL,
1051:                                        NULL,
1052:                                 /*119*/NULL,
1053:                                        NULL,
1054:                                        NULL,
1055:                                        NULL,
1056:                                        NULL,
1057:                                 /*124*/NULL,
1058:                                        NULL,
1059:                                        MatInvertBlockDiagonal_MPISELL,
1060:                                        NULL,
1061:                                        NULL,
1062:                                 /*129*/NULL,
1063:                                        NULL,
1064:                                        NULL,
1065:                                        NULL,
1066:                                        NULL,
1067:                                 /*134*/NULL,
1068:                                        NULL,
1069:                                        NULL,
1070:                                        NULL,
1071:                                        NULL,
1072:                                 /*139*/NULL,
1073:                                        NULL,
1074:                                        NULL,
1075:                                        MatFDColoringSetUp_MPIXAIJ,
1076:                                        NULL,
1077:                                 /*144*/NULL,
1078:                                        NULL,
1079:                                        NULL,
1080:                                        NULL
1081: };

1083: /* ----------------------------------------------------------------------------------------*/

1085: PetscErrorCode MatStoreValues_MPISELL(Mat mat)
1086: {
1087:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

1089:   MatStoreValues(sell->A);
1090:   MatStoreValues(sell->B);
1091:   return 0;
1092: }

1094: PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
1095: {
1096:   Mat_MPISELL    *sell=(Mat_MPISELL*)mat->data;

1098:   MatRetrieveValues(sell->A);
1099:   MatRetrieveValues(sell->B);
1100:   return 0;
1101: }

1103: PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[])
1104: {
1105:   Mat_MPISELL    *b;

1107:   PetscLayoutSetUp(B->rmap);
1108:   PetscLayoutSetUp(B->cmap);
1109:   b = (Mat_MPISELL*)B->data;

1111:   if (!B->preallocated) {
1112:     /* Explicitly create 2 MATSEQSELL matrices. */
1113:     MatCreate(PETSC_COMM_SELF,&b->A);
1114:     MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
1115:     MatSetBlockSizesFromMats(b->A,B,B);
1116:     MatSetType(b->A,MATSEQSELL);
1117:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
1118:     MatCreate(PETSC_COMM_SELF,&b->B);
1119:     MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
1120:     MatSetBlockSizesFromMats(b->B,B,B);
1121:     MatSetType(b->B,MATSEQSELL);
1122:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
1123:   }

1125:   MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);
1126:   MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);
1127:   B->preallocated  = PETSC_TRUE;
1128:   B->was_assembled = PETSC_FALSE;
1129:   /*
1130:     critical for MatAssemblyEnd to work.
1131:     MatAssemblyBegin checks it to set up was_assembled
1132:     and MatAssemblyEnd checks was_assembled to determine whether to build garray
1133:   */
1134:   B->assembled     = PETSC_FALSE;
1135:   return 0;
1136: }

1138: PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat)
1139: {
1140:   Mat            mat;
1141:   Mat_MPISELL    *a,*oldmat=(Mat_MPISELL*)matin->data;

1143:   *newmat = NULL;
1144:   MatCreate(PetscObjectComm((PetscObject)matin),&mat);
1145:   MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);
1146:   MatSetBlockSizesFromMats(mat,matin,matin);
1147:   MatSetType(mat,((PetscObject)matin)->type_name);
1148:   a       = (Mat_MPISELL*)mat->data;

1150:   mat->factortype   = matin->factortype;
1151:   mat->assembled    = PETSC_TRUE;
1152:   mat->insertmode   = NOT_SET_VALUES;
1153:   mat->preallocated = PETSC_TRUE;

1155:   a->size         = oldmat->size;
1156:   a->rank         = oldmat->rank;
1157:   a->donotstash   = oldmat->donotstash;
1158:   a->roworiented  = oldmat->roworiented;
1159:   a->rowindices   = NULL;
1160:   a->rowvalues    = NULL;
1161:   a->getrowactive = PETSC_FALSE;

1163:   PetscLayoutReference(matin->rmap,&mat->rmap);
1164:   PetscLayoutReference(matin->cmap,&mat->cmap);

1166:   if (oldmat->colmap) {
1167: #if defined(PETSC_USE_CTABLE)
1168:     PetscTableCreateCopy(oldmat->colmap,&a->colmap);
1169: #else
1170:     PetscMalloc1(mat->cmap->N,&a->colmap);
1171:     PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));
1172:     PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N);
1173: #endif
1174:   } else a->colmap = NULL;
1175:   if (oldmat->garray) {
1176:     PetscInt len;
1177:     len  = oldmat->B->cmap->n;
1178:     PetscMalloc1(len+1,&a->garray);
1179:     PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));
1180:     if (len) PetscArraycpy(a->garray,oldmat->garray,len);
1181:   } else a->garray = NULL;

1183:   VecDuplicate(oldmat->lvec,&a->lvec);
1184:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);
1185:   VecScatterCopy(oldmat->Mvctx,&a->Mvctx);
1186:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);
1187:   MatDuplicate(oldmat->A,cpvalues,&a->A);
1188:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);
1189:   MatDuplicate(oldmat->B,cpvalues,&a->B);
1190:   PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);
1191:   PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);
1192:   *newmat = mat;
1193:   return 0;
1194: }

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

1201:    Collective

1203:    Input Parameters:
1204: +  B - the matrix
1205: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1206:            (same value is used for all local rows)
1207: .  d_nnz - array containing the number of nonzeros in the various rows of the
1208:            DIAGONAL portion of the local submatrix (possibly different for each row)
1209:            or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1210:            The size of this array is equal to the number of local rows, i.e 'm'.
1211:            For matrices that will be factored, you must leave room for (and set)
1212:            the diagonal entry even if it is zero.
1213: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1214:            submatrix (same value is used for all local rows).
1215: -  o_nnz - array containing the number of nonzeros in the various rows of the
1216:            OFF-DIAGONAL portion of the local submatrix (possibly different for
1217:            each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1218:            structure. The size of this array is equal to the number
1219:            of local rows, i.e 'm'.

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

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

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

1229:    The DIAGONAL portion of the local submatrix of a processor can be defined
1230:    as the submatrix which is obtained by extraction the part corresponding to
1231:    the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
1232:    first row that belongs to the processor, r2 is the last row belonging to
1233:    the this processor, and c1-c2 is range of indices of the local part of a
1234:    vector suitable for applying the matrix to.  This is an mxn matrix.  In the
1235:    common case of a square matrix, the row and column ranges are the same and
1236:    the DIAGONAL part is also square. The remaining portion of the local
1237:    submatrix (mxN) constitute the OFF-DIAGONAL portion.

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

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

1246:    Example usage:

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

1253: .vb
1254:             1  2  0  |  0  3  0  |  0  4
1255:     Proc0   0  5  6  |  7  0  0  |  8  0
1256:             9  0 10  | 11  0  0  | 12  0
1257:     -------------------------------------
1258:            13  0 14  | 15 16 17  |  0  0
1259:     Proc1   0 18  0  | 19 20 21  |  0  0
1260:             0  0  0  | 22 23  0  | 24  0
1261:     -------------------------------------
1262:     Proc2  25 26 27  |  0  0 28  | 29  0
1263:            30  0  0  | 31 32 33  |  0 34
1264: .ve

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

1268: .vb
1269:       A B C
1270:       D E F
1271:       G H I
1272: .ve

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

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

1281:    The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1282:    submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1283:    corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1284:    Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1285:    part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1286:    matrix, ans [DF] as another SeqSELL matrix.

1288:    When d_nz, o_nz parameters are specified, d_nz storage elements are
1289:    allocated for every row of the local diagonal submatrix, and o_nz
1290:    storage locations are allocated for every row of the OFF-DIAGONAL submat.
1291:    One way to choose d_nz and o_nz is to use the max nonzerors per local
1292:    rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1293:    In this case, the values of d_nz,o_nz are:
1294: .vb
1295:      proc0 : dnz = 2, o_nz = 2
1296:      proc1 : dnz = 3, o_nz = 2
1297:      proc2 : dnz = 1, o_nz = 4
1298: .ve
1299:    We are allocating m*(d_nz+o_nz) storage locations for every proc. This
1300:    translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1301:    for proc3. i.e we are using 12+15+10=37 storage locations to store
1302:    34 values.

1304:    When d_nnz, o_nnz parameters are specified, the storage is specified
1305:    for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1306:    In the above case the values for d_nnz,o_nnz are:
1307: .vb
1308:      proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
1309:      proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
1310:      proc2: d_nnz = [1,1]   and o_nnz = [4,4]
1311: .ve
1312:    Here the space allocated is according to nz (or maximum values in the nnz
1313:    if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37

1315:    Level: intermediate

1317: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(),
1318:           MATMPISELL, MatGetInfo(), PetscSplitOwnership()
1319: @*/
1320: PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1321: {
1324:   PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1325:   return 0;
1326: }

1328: /*MC
1329:    MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices,
1330:    based on the sliced Ellpack format

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

1335:    Level: beginner

1337: .seealso: MatCreateSell(), MATSEQSELL, MATSELL, MATSEQAIJ, MATAIJ, MATMPIAIJ
1338: M*/

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

1343:    Collective

1345:    Input Parameters:
1346: +  comm - MPI communicator
1347: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
1348:            This value should be the same as the local size used in creating the
1349:            y vector for the matrix-vector product y = Ax.
1350: .  n - This value should be the same as the local size used in creating the
1351:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
1352:        calculated if N is given) For square matrices n is almost always m.
1353: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
1354: .  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
1355: .  d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix
1356:                (same value is used for all local rows)
1357: .  d_rlen - array containing the number of nonzeros in the various rows of the
1358:             DIAGONAL portion of the local submatrix (possibly different for each row)
1359:             or NULL, if d_rlenmax is used to specify the nonzero structure.
1360:             The size of this array is equal to the number of local rows, i.e 'm'.
1361: .  o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local
1362:                submatrix (same value is used for all local rows).
1363: -  o_rlen - array containing the number of nonzeros in the various rows of the
1364:             OFF-DIAGONAL portion of the local submatrix (possibly different for
1365:             each row) or NULL, if o_rlenmax is used to specify the nonzero
1366:             structure. The size of this array is equal to the number
1367:             of local rows, i.e 'm'.

1369:    Output Parameter:
1370: .  A - the matrix

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

1376:    Notes:
1377:    If the *_rlen parameter is given then the *_rlenmax parameter is ignored

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

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

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

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

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

1400:    The DIAGONAL portion of the local submatrix on any given processor
1401:    is the submatrix corresponding to the rows and columns m,n
1402:    corresponding to the given processor. i.e diagonal matrix on
1403:    process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
1404:    etc. The remaining portion of the local submatrix [m x (N-n)]
1405:    constitute the OFF-DIAGONAL portion. The example below better
1406:    illustrates this concept.

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

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

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

1420:    Options Database Keys:
1421: -  -mat_sell_oneindex - Internally use indexing starting at 1
1422:         rather than 0.  Note that when calling MatSetValues(),
1423:         the user still MUST index entries starting at 0!

1425:    Example usage:

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

1432: .vb
1433:             1  2  0  |  0  3  0  |  0  4
1434:     Proc0   0  5  6  |  7  0  0  |  8  0
1435:             9  0 10  | 11  0  0  | 12  0
1436:     -------------------------------------
1437:            13  0 14  | 15 16 17  |  0  0
1438:     Proc1   0 18  0  | 19 20 21  |  0  0
1439:             0  0  0  | 22 23  0  | 24  0
1440:     -------------------------------------
1441:     Proc2  25 26 27  |  0  0 28  | 29  0
1442:            30  0  0  | 31 32 33  |  0 34
1443: .ve

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

1447: .vb
1448:       A B C
1449:       D E F
1450:       G H I
1451: .ve

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

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

1460:    The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1461:    submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1462:    corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1463:    Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1464:    part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1465:    matrix, ans [DF] as another SeqSELL matrix.

1467:    When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are
1468:    allocated for every row of the local diagonal submatrix, and o_rlenmax
1469:    storage locations are allocated for every row of the OFF-DIAGONAL submat.
1470:    One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local
1471:    rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1472:    In this case, the values of d_rlenmax,o_rlenmax are:
1473: .vb
1474:      proc0 : d_rlenmax = 2, o_rlenmax = 2
1475:      proc1 : d_rlenmax = 3, o_rlenmax = 2
1476:      proc2 : d_rlenmax = 1, o_rlenmax = 4
1477: .ve
1478:    We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This
1479:    translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1480:    for proc3. i.e we are using 12+15+10=37 storage locations to store
1481:    34 values.

1483:    When d_rlen, o_rlen parameters are specified, the storage is specified
1484:    for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1485:    In the above case the values for d_nnz,o_nnz are:
1486: .vb
1487:      proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
1488:      proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
1489:      proc2: d_nnz = [1,1]   and o_nnz = [4,4]
1490: .ve
1491:    Here the space allocated is still 37 though there are 34 nonzeros because
1492:    the allocation is always done according to rlenmax.

1494:    Level: intermediate

1496: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(),
1497:           MATMPISELL, MatCreateMPISELLWithArrays()
1498: @*/
1499: 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)
1500: {
1501:   PetscMPIInt    size;

1503:   MatCreate(comm,A);
1504:   MatSetSizes(*A,m,n,M,N);
1505:   MPI_Comm_size(comm,&size);
1506:   if (size > 1) {
1507:     MatSetType(*A,MATMPISELL);
1508:     MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);
1509:   } else {
1510:     MatSetType(*A,MATSEQSELL);
1511:     MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);
1512:   }
1513:   return 0;
1514: }

1516: PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[])
1517: {
1518:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1519:   PetscBool      flg;

1521:   PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);
1523:   if (Ad)     *Ad     = a->A;
1524:   if (Ao)     *Ao     = a->B;
1525:   if (colmap) *colmap = a->garray;
1526:   return 0;
1527: }

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

1532:     Not Collective

1534:    Input Parameters:
1535: +    A - the matrix
1536: .    scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
1537: -    row, col - index sets of rows and columns to extract (or NULL)

1539:    Output Parameter:
1540: .    A_loc - the local sequential matrix generated

1542:     Level: developer

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

1546: @*/
1547: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc)
1548: {
1549:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1550:   PetscInt       i,start,end,ncols,nzA,nzB,*cmap,imark,*idx;
1551:   IS             isrowa,iscola;
1552:   Mat            *aloc;
1553:   PetscBool      match;

1555:   PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match);
1557:   PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0);
1558:   if (!row) {
1559:     start = A->rmap->rstart; end = A->rmap->rend;
1560:     ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa);
1561:   } else {
1562:     isrowa = *row;
1563:   }
1564:   if (!col) {
1565:     start = A->cmap->rstart;
1566:     cmap  = a->garray;
1567:     nzA   = a->A->cmap->n;
1568:     nzB   = a->B->cmap->n;
1569:     PetscMalloc1(nzA+nzB, &idx);
1570:     ncols = 0;
1571:     for (i=0; i<nzB; i++) {
1572:       if (cmap[i] < start) idx[ncols++] = cmap[i];
1573:       else break;
1574:     }
1575:     imark = i;
1576:     for (i=0; i<nzA; i++) idx[ncols++] = start + i;
1577:     for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i];
1578:     ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola);
1579:   } else {
1580:     iscola = *col;
1581:   }
1582:   if (scall != MAT_INITIAL_MATRIX) {
1583:     PetscMalloc1(1,&aloc);
1584:     aloc[0] = *A_loc;
1585:   }
1586:   MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc);
1587:   *A_loc = aloc[0];
1588:   PetscFree(aloc);
1589:   if (!row) {
1590:     ISDestroy(&isrowa);
1591:   }
1592:   if (!col) {
1593:     ISDestroy(&iscola);
1594:   }
1595:   PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0);
1596:   return 0;
1597: }

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

1601: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1602: {
1603:   Mat_MPISELL    *a=(Mat_MPISELL*)A->data;
1604:   Mat            B;
1605:   Mat_MPIAIJ     *b;


1609:   if (reuse == MAT_REUSE_MATRIX) {
1610:     B = *newmat;
1611:   } else {
1612:     MatCreate(PetscObjectComm((PetscObject)A),&B);
1613:     MatSetType(B,MATMPIAIJ);
1614:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1615:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1616:     MatSeqAIJSetPreallocation(B,0,NULL);
1617:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1618:   }
1619:   b    = (Mat_MPIAIJ*) B->data;

1621:   if (reuse == MAT_REUSE_MATRIX) {
1622:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);
1623:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);
1624:   } else {
1625:     MatDestroy(&b->A);
1626:     MatDestroy(&b->B);
1627:     MatDisAssemble_MPISELL(A);
1628:     MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);
1629:     MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);
1630:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1631:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1632:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1633:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1634:   }

1636:   if (reuse == MAT_INPLACE_MATRIX) {
1637:     MatHeaderReplace(A,&B);
1638:   } else {
1639:     *newmat = B;
1640:   }
1641:   return 0;
1642: }

1644: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1645: {
1646:   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)A->data;
1647:   Mat            B;
1648:   Mat_MPISELL    *b;


1652:   if (reuse == MAT_REUSE_MATRIX) {
1653:     B = *newmat;
1654:   } else {
1655:     MatCreate(PetscObjectComm((PetscObject)A),&B);
1656:     MatSetType(B,MATMPISELL);
1657:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1658:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1659:     MatSeqAIJSetPreallocation(B,0,NULL);
1660:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1661:   }
1662:   b    = (Mat_MPISELL*) B->data;

1664:   if (reuse == MAT_REUSE_MATRIX) {
1665:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);
1666:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);
1667:   } else {
1668:     MatDestroy(&b->A);
1669:     MatDestroy(&b->B);
1670:     MatDisAssemble_MPIAIJ(A);
1671:     MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);
1672:     MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);
1673:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1674:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1675:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1676:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1677:   }

1679:   if (reuse == MAT_INPLACE_MATRIX) {
1680:     MatHeaderReplace(A,&B);
1681:   } else {
1682:     *newmat = B;
1683:   }
1684:   return 0;
1685: }

1687: PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1688: {
1689:   Mat_MPISELL    *mat=(Mat_MPISELL*)matin->data;
1690:   Vec            bb1=NULL;

1692:   if (flag == SOR_APPLY_UPPER) {
1693:     (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1694:     return 0;
1695:   }

1697:   if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) {
1698:     VecDuplicate(bb,&bb1);
1699:   }

1701:   if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1702:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1703:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1704:       its--;
1705:     }

1707:     while (its--) {
1708:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1709:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1711:       /* update rhs: bb1 = bb - B*x */
1712:       VecScale(mat->lvec,-1.0);
1713:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

1715:       /* local sweep */
1716:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);
1717:     }
1718:   } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
1719:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1720:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1721:       its--;
1722:     }
1723:     while (its--) {
1724:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1725:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1727:       /* update rhs: bb1 = bb - B*x */
1728:       VecScale(mat->lvec,-1.0);
1729:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

1731:       /* local sweep */
1732:       (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);
1733:     }
1734:   } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1735:     if (flag & SOR_ZERO_INITIAL_GUESS) {
1736:       (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1737:       its--;
1738:     }
1739:     while (its--) {
1740:       VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1741:       VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);

1743:       /* update rhs: bb1 = bb - B*x */
1744:       VecScale(mat->lvec,-1.0);
1745:       (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);

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

1752:   VecDestroy(&bb1);

1754:   matin->factorerrortype = mat->A->factorerrortype;
1755:   return 0;
1756: }

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

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

1764:   Level: beginner

1766: .seealso: MatCreateSELL()
1767: M*/
1768: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1769: {
1770:   Mat_MPISELL    *b;
1771:   PetscMPIInt    size;

1773:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1774:   PetscNewLog(B,&b);
1775:   B->data       = (void*)b;
1776:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1777:   B->assembled  = PETSC_FALSE;
1778:   B->insertmode = NOT_SET_VALUES;
1779:   b->size       = size;
1780:   MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);
1781:   /* build cache for off array entries formed */
1782:   MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);

1784:   b->donotstash  = PETSC_FALSE;
1785:   b->colmap      = NULL;
1786:   b->garray      = NULL;
1787:   b->roworiented = PETSC_TRUE;

1789:   /* stuff used for matrix vector multiply */
1790:   b->lvec  = NULL;
1791:   b->Mvctx = NULL;

1793:   /* stuff for MatGetRow() */
1794:   b->rowindices   = NULL;
1795:   b->rowvalues    = NULL;
1796:   b->getrowactive = PETSC_FALSE;

1798:   PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);
1799:   PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);
1800:   PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);
1801:   PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);
1802:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);
1803:   PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);
1804:   PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);
1805:   return 0;
1806: }