Actual source code: mpisell.c
petsc-3.11.4 2019-09-28
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,¬me);
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: }