Actual source code: mpisell.c
petsc-3.9.4 2018-09-11
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: Subclasses include MATSELLCUSP, MATSELLCUSPARSE, MATSELLPERM, MATSELLCRL, and also automatically switches over to use inodes when
21: enough exist.
23: Level: beginner
25: .seealso: MatCreateSELL(), MatCreateSeqSELL(), MATSEQSELL, MATMPISELL
26: M*/
28: PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is)
29: {
31: Mat_MPISELL *sell=(Mat_MPISELL*)Y->data;
34: if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) {
35: MatDiagonalSet(sell->A,D,is);
36: } else {
37: MatDiagonalSet_Default(Y,D,is);
38: }
39: return(0);
40: }
42: /*
43: Local utility routine that creates a mapping from the global column
44: number to the local number in the off-diagonal part of the local
45: storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at
46: a slightly higher hash table cost; without it it is not scalable (each processor
47: has an order N integer array but is fast to acess.
48: */
49: PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat)
50: {
51: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
53: PetscInt n=sell->B->cmap->n,i;
56: if (!sell->garray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MPISELL Matrix was assembled but is missing garray");
57: #if defined(PETSC_USE_CTABLE)
58: PetscTableCreate(n,mat->cmap->N+1,&sell->colmap);
59: for (i=0; i<n; i++) {
60: PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES);
61: }
62: #else
63: PetscCalloc1(mat->cmap->N+1,&sell->colmap);
64: PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt));
65: for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1;
66: #endif
67: return(0);
68: }
70: #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \
71: { \
72: if (col <= lastcol1) low1 = 0; \
73: else high1 = nrow1; \
74: lastcol1 = col; \
75: while (high1-low1 > 5) { \
76: t = (low1+high1)/2; \
77: if (*(cp1+8*t) > col) high1 = t; \
78: else low1 = t; \
79: } \
80: for (_i=low1; _i<high1; _i++) { \
81: if (*(cp1+8*_i) > col) break; \
82: if (*(cp1+8*_i) == col) { \
83: if (addv == ADD_VALUES) *(vp1+8*_i) += value; \
84: else *(vp1+8*_i) = value; \
85: goto a_noinsert; \
86: } \
87: } \
88: if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \
89: if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \
90: 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); \
91: MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \
92: /* shift up all the later entries in this row */ \
93: for (ii=nrow1-1; ii>=_i; ii--) { \
94: *(cp1+8*(ii+1)) = *(cp1+8*ii); \
95: *(vp1+8*(ii+1)) = *(vp1+8*ii); \
96: } \
97: *(cp1+8*_i) = col; \
98: *(vp1+8*_i) = value; \
99: a->nz++; nrow1++; A->nonzerostate++; \
100: a_noinsert: ; \
101: a->rlen[row] = nrow1; \
102: }
104: #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \
105: { \
106: if (col <= lastcol2) low2 = 0; \
107: else high2 = nrow2; \
108: lastcol2 = col; \
109: while (high2-low2 > 5) { \
110: t = (low2+high2)/2; \
111: if (*(cp2+8*t) > col) high2 = t; \
112: else low2 = t; \
113: } \
114: for (_i=low2; _i<high2; _i++) { \
115: if (*(cp2+8*_i) > col) break; \
116: if (*(cp2+8*_i) == col) { \
117: if (addv == ADD_VALUES) *(vp2+8*_i) += value; \
118: else *(vp2+8*_i) = value; \
119: goto b_noinsert; \
120: } \
121: } \
122: if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
123: if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
124: 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); \
125: MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \
126: /* shift up all the later entries in this row */ \
127: for (ii=nrow2-1; ii>=_i; ii--) { \
128: *(cp2+8*(ii+1)) = *(cp2+8*ii); \
129: *(vp2+8*(ii+1)) = *(vp2+8*ii); \
130: } \
131: *(cp2+8*_i) = col; \
132: *(vp2+8*_i) = value; \
133: b->nz++; nrow2++; B->nonzerostate++; \
134: b_noinsert: ; \
135: b->rlen[row] = nrow2; \
136: }
138: PetscErrorCode MatSetValues_MPISELL(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
139: {
140: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
141: PetscScalar value;
143: PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend,shift1,shift2;
144: PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;
145: PetscBool roworiented=sell->roworiented;
147: /* Some Variables required in the macro */
148: Mat A=sell->A;
149: Mat_SeqSELL *a=(Mat_SeqSELL*)A->data;
150: PetscBool ignorezeroentries=a->ignorezeroentries,found;
151: Mat B=sell->B;
152: Mat_SeqSELL *b=(Mat_SeqSELL*)B->data;
153: PetscInt *cp1,*cp2,ii,_i,nrow1,nrow2,low1,high1,low2,high2,t,lastcol1,lastcol2;
154: MatScalar *vp1,*vp2;
157: for (i=0; i<m; i++) {
158: if (im[i] < 0) continue;
159: #if defined(PETSC_USE_DEBUG)
160: 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);
161: #endif
162: if (im[i] >= rstart && im[i] < rend) {
163: row = im[i] - rstart;
164: lastcol1 = -1;
165: shift1 = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */
166: cp1 = a->colidx+shift1;
167: vp1 = a->val+shift1;
168: nrow1 = a->rlen[row];
169: low1 = 0;
170: high1 = nrow1;
171: lastcol2 = -1;
172: shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
173: cp2 = b->colidx+shift2;
174: vp2 = b->val+shift2;
175: nrow2 = b->rlen[row];
176: low2 = 0;
177: high2 = nrow2;
179: for (j=0; j<n; j++) {
180: if (roworiented) value = v[i*n+j];
181: else value = v[i+j*m];
182: if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue;
183: if (in[j] >= cstart && in[j] < cend) {
184: col = in[j] - cstart;
185: MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */
186: } else if (in[j] < 0) continue;
187: #if defined(PETSC_USE_DEBUG)
188: 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);
189: #endif
190: else {
191: if (mat->was_assembled) {
192: if (!sell->colmap) {
193: MatCreateColmap_MPISELL_Private(mat);
194: }
195: #if defined(PETSC_USE_CTABLE)
196: PetscTableFind(sell->colmap,in[j]+1,&col);
197: col--;
198: #else
199: col = sell->colmap[in[j]] - 1;
200: #endif
201: if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) {
202: MatDisAssemble_MPISELL(mat);
203: col = in[j];
204: /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */
205: B = sell->B;
206: b = (Mat_SeqSELL*)B->data;
207: shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
208: cp2 = b->colidx+shift2;
209: vp2 = b->val+shift2;
210: nrow2 = b->rlen[row];
211: low2 = 0;
212: high2 = nrow2;
213: } 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]);
214: } else col = in[j];
215: MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */
216: }
217: }
218: } else {
219: 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]);
220: if (!sell->donotstash) {
221: mat->assembled = PETSC_FALSE;
222: if (roworiented) {
223: MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
224: } else {
225: MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
226: }
227: }
228: }
229: }
230: return(0);
231: }
233: PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
234: {
235: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
237: PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend;
238: PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;
241: for (i=0; i<m; i++) {
242: if (idxm[i] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",idxm[i]);*/
243: 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);
244: if (idxm[i] >= rstart && idxm[i] < rend) {
245: row = idxm[i] - rstart;
246: for (j=0; j<n; j++) {
247: if (idxn[j] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",idxn[j]); */
248: 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);
249: if (idxn[j] >= cstart && idxn[j] < cend) {
250: col = idxn[j] - cstart;
251: MatGetValues(sell->A,1,&row,1,&col,v+i*n+j);
252: } else {
253: if (!sell->colmap) {
254: MatCreateColmap_MPISELL_Private(mat);
255: }
256: #if defined(PETSC_USE_CTABLE)
257: PetscTableFind(sell->colmap,idxn[j]+1,&col);
258: col--;
259: #else
260: col = sell->colmap[idxn[j]] - 1;
261: #endif
262: if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0;
263: else {
264: MatGetValues(sell->B,1,&row,1,&col,v+i*n+j);
265: }
266: }
267: }
268: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported");
269: }
270: return(0);
271: }
273: extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec);
275: PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode)
276: {
277: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
279: PetscInt nstash,reallocs;
282: if (sell->donotstash || mat->nooffprocentries) return(0);
284: MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
285: MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);
286: PetscInfo2(sell->A,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);
287: return(0);
288: }
290: PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode)
291: {
292: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
294: PetscMPIInt n;
295: PetscInt i,flg;
296: PetscInt *row,*col;
297: PetscScalar *val;
298: PetscBool other_disassembled;
299: /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */
301: if (!sell->donotstash && !mat->nooffprocentries) {
302: while (1) {
303: MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);
304: if (!flg) break;
306: for (i=0; i<n; i++) { /* assemble one by one */
307: MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode);
308: }
309: }
310: MatStashScatterEnd_Private(&mat->stash);
311: }
312: MatAssemblyBegin(sell->A,mode);
313: MatAssemblyEnd(sell->A,mode);
315: /*
316: determine if any processor has disassembled, if so we must
317: also disassemble ourselfs, in order that we may reassemble.
318: */
319: /*
320: if nonzero structure of submatrix B cannot change then we know that
321: no processor disassembled thus we can skip this stuff
322: */
323: if (!((Mat_SeqSELL*)sell->B->data)->nonew) {
324: MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat));
325: if (mat->was_assembled && !other_disassembled) {
326: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet\n");
327: MatDisAssemble_MPISELL(mat);
328: }
329: }
330: if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) {
331: MatSetUpMultiply_MPISELL(mat);
332: }
333: /*
334: MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE);
335: */
336: MatAssemblyBegin(sell->B,mode);
337: MatAssemblyEnd(sell->B,mode);
338: PetscFree2(sell->rowvalues,sell->rowindices);
339: sell->rowvalues = 0;
340: VecDestroy(&sell->diag);
342: /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
343: if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) {
344: PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate;
345: MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));
346: }
347: return(0);
348: }
350: PetscErrorCode MatZeroEntries_MPISELL(Mat A)
351: {
352: Mat_MPISELL *l=(Mat_MPISELL*)A->data;
356: MatZeroEntries(l->A);
357: MatZeroEntries(l->B);
358: return(0);
359: }
361: PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy)
362: {
363: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
365: PetscInt nt;
368: VecGetLocalSize(xx,&nt);
369: 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);
370: VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
371: (*a->A->ops->mult)(a->A,xx,yy);
372: VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
373: (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);
374: return(0);
375: }
377: PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx)
378: {
379: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
383: MatMultDiagonalBlock(a->A,bb,xx);
384: return(0);
385: }
387: PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
388: {
389: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
393: VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
394: (*a->A->ops->multadd)(a->A,xx,yy,zz);
395: VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
396: (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
397: return(0);
398: }
400: PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy)
401: {
402: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
404: PetscBool merged;
407: VecScatterGetMerged(a->Mvctx,&merged);
408: /* do nondiagonal part */
409: (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
410: if (!merged) {
411: /* send it on its way */
412: VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
413: /* do local part */
414: (*a->A->ops->multtranspose)(a->A,xx,yy);
415: /* receive remote parts: note this assumes the values are not actually */
416: /* added in yy until the next line, */
417: VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
418: } else {
419: /* do local part */
420: (*a->A->ops->multtranspose)(a->A,xx,yy);
421: /* send it on its way */
422: VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
423: /* values actually were received in the Begin() but we need to call this nop */
424: VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
425: }
426: return(0);
427: }
429: PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f)
430: {
431: MPI_Comm comm;
432: Mat_MPISELL *Asell=(Mat_MPISELL*)Amat->data,*Bsell;
433: Mat Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs;
434: IS Me,Notme;
436: PetscInt M,N,first,last,*notme,i;
437: PetscMPIInt size;
440: /* Easy test: symmetric diagonal block */
441: Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A;
442: MatIsTranspose(Adia,Bdia,tol,f);
443: if (!*f) return(0);
444: PetscObjectGetComm((PetscObject)Amat,&comm);
445: MPI_Comm_size(comm,&size);
446: if (size == 1) return(0);
448: /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */
449: MatGetSize(Amat,&M,&N);
450: MatGetOwnershipRange(Amat,&first,&last);
451: PetscMalloc1(N-last+first,¬me);
452: for (i=0; i<first; i++) notme[i] = i;
453: for (i=last; i<M; i++) notme[i-last+first] = i;
454: ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme);
455: ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me);
456: MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs);
457: Aoff = Aoffs[0];
458: MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs);
459: Boff = Boffs[0];
460: MatIsTranspose(Aoff,Boff,tol,f);
461: MatDestroyMatrices(1,&Aoffs);
462: MatDestroyMatrices(1,&Boffs);
463: ISDestroy(&Me);
464: ISDestroy(&Notme);
465: PetscFree(notme);
466: return(0);
467: }
469: PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
470: {
471: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
475: /* do nondiagonal part */
476: (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
477: /* send it on its way */
478: VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
479: /* do local part */
480: (*a->A->ops->multtransposeadd)(a->A,xx,yy,zz);
481: /* receive remote parts */
482: VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
483: return(0);
484: }
486: /*
487: This only works correctly for square matrices where the subblock A->A is the
488: diagonal block
489: */
490: PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v)
491: {
493: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
496: if (A->rmap->N != A->cmap->N) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block");
497: 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");
498: MatGetDiagonal(a->A,v);
499: return(0);
500: }
502: PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa)
503: {
504: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
508: MatScale(a->A,aa);
509: MatScale(a->B,aa);
510: return(0);
511: }
513: PetscErrorCode MatDestroy_MPISELL(Mat mat)
514: {
515: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
519: #if defined(PETSC_USE_LOG)
520: PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D",mat->rmap->N,mat->cmap->N);
521: #endif
522: MatStashDestroy_Private(&mat->stash);
523: VecDestroy(&sell->diag);
524: MatDestroy(&sell->A);
525: MatDestroy(&sell->B);
526: #if defined(PETSC_USE_CTABLE)
527: PetscTableDestroy(&sell->colmap);
528: #else
529: PetscFree(sell->colmap);
530: #endif
531: PetscFree(sell->garray);
532: VecDestroy(&sell->lvec);
533: VecScatterDestroy(&sell->Mvctx);
534: PetscFree2(sell->rowvalues,sell->rowindices);
535: PetscFree(sell->ld);
536: PetscFree(mat->data);
538: PetscObjectChangeTypeName((PetscObject)mat,0);
539: PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);
540: PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);
541: PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL);
542: PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL);
543: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL);
544: PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);
545: return(0);
546: }
548: #include <petscdraw.h>
549: PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer)
550: {
551: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
552: PetscErrorCode ierr;
553: PetscMPIInt rank=sell->rank,size=sell->size;
554: PetscBool isdraw,iascii,isbinary;
555: PetscViewer sviewer;
556: PetscViewerFormat format;
559: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
560: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
561: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
562: if (iascii) {
563: PetscViewerGetFormat(viewer,&format);
564: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
565: MatInfo info;
566: PetscBool inodes;
568: MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank);
569: MatGetInfo(mat,MAT_LOCAL,&info);
570: MatInodeGetInodeSizes(sell->A,NULL,(PetscInt**)&inodes,NULL);
571: PetscViewerASCIIPushSynchronized(viewer);
572: if (!inodes) {
573: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D mem %D, not using I-node routines\n",
574: rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
575: } else {
576: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D mem %D, using I-node routines\n",
577: rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
578: }
579: MatGetInfo(sell->A,MAT_LOCAL,&info);
580: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);
581: MatGetInfo(sell->B,MAT_LOCAL,&info);
582: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);
583: PetscViewerFlush(viewer);
584: PetscViewerASCIIPopSynchronized(viewer);
585: PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n");
586: VecScatterView(sell->Mvctx,viewer);
587: return(0);
588: } else if (format == PETSC_VIEWER_ASCII_INFO) {
589: PetscInt inodecount,inodelimit,*inodes;
590: MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit);
591: if (inodes) {
592: PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %D nodes, limit used is %D\n",inodecount,inodelimit);
593: } else {
594: PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n");
595: }
596: return(0);
597: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
598: return(0);
599: }
600: } else if (isbinary) {
601: if (size == 1) {
602: PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name);
603: MatView(sell->A,viewer);
604: } else {
605: /* MatView_MPISELL_Binary(mat,viewer); */
606: }
607: return(0);
608: } else if (isdraw) {
609: PetscDraw draw;
610: PetscBool isnull;
611: PetscViewerDrawGetDraw(viewer,0,&draw);
612: PetscDrawIsNull(draw,&isnull);
613: if (isnull) return(0);
614: }
616: {
617: /* assemble the entire matrix onto first processor. */
618: Mat A;
619: Mat_SeqSELL *Aloc;
620: PetscInt M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j;
621: MatScalar *aval;
622: PetscBool isnonzero;
624: MatCreate(PetscObjectComm((PetscObject)mat),&A);
625: if (!rank) {
626: MatSetSizes(A,M,N,M,N);
627: } else {
628: MatSetSizes(A,0,0,M,N);
629: }
630: /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */
631: MatSetType(A,MATMPISELL);
632: MatMPISELLSetPreallocation(A,0,NULL,0,NULL);
633: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
634: PetscLogObjectParent((PetscObject)mat,(PetscObject)A);
636: /* copy over the A part */
637: Aloc = (Mat_SeqSELL*)sell->A->data;
638: acolidx = Aloc->colidx; aval = Aloc->val;
639: for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */
640: for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
641: isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
642: if (isnonzero) { /* check the mask bit */
643: row = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */
644: col = *acolidx + mat->rmap->rstart;
645: MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
646: }
647: aval++; acolidx++;
648: }
649: }
651: /* copy over the B part */
652: Aloc = (Mat_SeqSELL*)sell->B->data;
653: acolidx = Aloc->colidx; aval = Aloc->val;
654: for (i=0; i<Aloc->totalslices; i++) {
655: for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
656: isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
657: if (isnonzero) {
658: row = (i<<3)+(j&0x07) + mat->rmap->rstart;
659: col = sell->garray[*acolidx];
660: MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
661: }
662: aval++; acolidx++;
663: }
664: }
666: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
667: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
668: /*
669: Everyone has to call to draw the matrix since the graphics waits are
670: synchronized across all processors that share the PetscDraw object
671: */
672: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
673: if (!rank) {
674: PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name);
675: MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer);
676: }
677: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
678: PetscViewerFlush(viewer);
679: MatDestroy(&A);
680: }
681: return(0);
682: }
684: PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer)
685: {
687: PetscBool iascii,isdraw,issocket,isbinary;
690: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
691: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
692: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
693: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
694: if (iascii || isdraw || isbinary || issocket) {
695: MatView_MPISELL_ASCIIorDraworSocket(mat,viewer);
696: }
697: return(0);
698: }
700: PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
701: {
702: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
706: MatGetSize(sell->B,NULL,nghosts);
707: if (ghosts) *ghosts = sell->garray;
708: return(0);
709: }
711: PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info)
712: {
713: Mat_MPISELL *mat=(Mat_MPISELL*)matin->data;
714: Mat A=mat->A,B=mat->B;
716: PetscReal isend[5],irecv[5];
719: info->block_size = 1.0;
720: MatGetInfo(A,MAT_LOCAL,info);
722: isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded;
723: isend[3] = info->memory; isend[4] = info->mallocs;
725: MatGetInfo(B,MAT_LOCAL,info);
727: isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded;
728: isend[3] += info->memory; isend[4] += info->mallocs;
729: if (flag == MAT_LOCAL) {
730: info->nz_used = isend[0];
731: info->nz_allocated = isend[1];
732: info->nz_unneeded = isend[2];
733: info->memory = isend[3];
734: info->mallocs = isend[4];
735: } else if (flag == MAT_GLOBAL_MAX) {
736: MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin));
738: info->nz_used = irecv[0];
739: info->nz_allocated = irecv[1];
740: info->nz_unneeded = irecv[2];
741: info->memory = irecv[3];
742: info->mallocs = irecv[4];
743: } else if (flag == MAT_GLOBAL_SUM) {
744: MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)matin));
746: info->nz_used = irecv[0];
747: info->nz_allocated = irecv[1];
748: info->nz_unneeded = irecv[2];
749: info->memory = irecv[3];
750: info->mallocs = irecv[4];
751: }
752: info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */
753: info->fill_ratio_needed = 0;
754: info->factor_mallocs = 0;
755: return(0);
756: }
758: PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg)
759: {
760: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
764: switch (op) {
765: case MAT_NEW_NONZERO_LOCATIONS:
766: case MAT_NEW_NONZERO_ALLOCATION_ERR:
767: case MAT_UNUSED_NONZERO_LOCATION_ERR:
768: case MAT_KEEP_NONZERO_PATTERN:
769: case MAT_NEW_NONZERO_LOCATION_ERR:
770: case MAT_USE_INODES:
771: case MAT_IGNORE_ZERO_ENTRIES:
772: MatCheckPreallocated(A,1);
773: MatSetOption(a->A,op,flg);
774: MatSetOption(a->B,op,flg);
775: break;
776: case MAT_ROW_ORIENTED:
777: MatCheckPreallocated(A,1);
778: a->roworiented = flg;
780: MatSetOption(a->A,op,flg);
781: MatSetOption(a->B,op,flg);
782: break;
783: case MAT_NEW_DIAGONALS:
784: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
785: break;
786: case MAT_IGNORE_OFF_PROC_ENTRIES:
787: a->donotstash = flg;
788: break;
789: case MAT_SPD:
790: A->spd_set = PETSC_TRUE;
791: A->spd = flg;
792: if (flg) {
793: A->symmetric = PETSC_TRUE;
794: A->structurally_symmetric = PETSC_TRUE;
795: A->symmetric_set = PETSC_TRUE;
796: A->structurally_symmetric_set = PETSC_TRUE;
797: }
798: break;
799: case MAT_SYMMETRIC:
800: MatCheckPreallocated(A,1);
801: MatSetOption(a->A,op,flg);
802: break;
803: case MAT_STRUCTURALLY_SYMMETRIC:
804: MatCheckPreallocated(A,1);
805: MatSetOption(a->A,op,flg);
806: break;
807: case MAT_HERMITIAN:
808: MatCheckPreallocated(A,1);
809: MatSetOption(a->A,op,flg);
810: break;
811: case MAT_SYMMETRY_ETERNAL:
812: MatCheckPreallocated(A,1);
813: MatSetOption(a->A,op,flg);
814: break;
815: default:
816: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
817: }
818: return(0);
819: }
822: PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr)
823: {
824: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
825: Mat a=sell->A,b=sell->B;
827: PetscInt s1,s2,s3;
830: MatGetLocalSize(mat,&s2,&s3);
831: if (rr) {
832: VecGetLocalSize(rr,&s1);
833: if (s1!=s3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size");
834: /* Overlap communication with computation. */
835: VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
836: }
837: if (ll) {
838: VecGetLocalSize(ll,&s1);
839: if (s1!=s2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size");
840: (*b->ops->diagonalscale)(b,ll,0);
841: }
842: /* scale the diagonal block */
843: (*a->ops->diagonalscale)(a,ll,rr);
845: if (rr) {
846: /* Do a scatter end and then right scale the off-diagonal block */
847: VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
848: (*b->ops->diagonalscale)(b,0,sell->lvec);
849: }
850: return(0);
851: }
853: PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
854: {
855: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
859: MatSetUnfactored(a->A);
860: return(0);
861: }
863: PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool *flag)
864: {
865: Mat_MPISELL *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data;
866: Mat a,b,c,d;
867: PetscBool flg;
871: a = matA->A; b = matA->B;
872: c = matB->A; d = matB->B;
874: MatEqual(a,c,&flg);
875: if (flg) {
876: MatEqual(b,d,&flg);
877: }
878: MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
879: return(0);
880: }
882: PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str)
883: {
885: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
886: Mat_MPISELL *b=(Mat_MPISELL*)B->data;
889: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
890: if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
891: /* because of the column compression in the off-processor part of the matrix a->B,
892: the number of columns in a->B and b->B may be different, hence we cannot call
893: the MatCopy() directly on the two parts. If need be, we can provide a more
894: efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices
895: then copying the submatrices */
896: MatCopy_Basic(A,B,str);
897: } else {
898: MatCopy(a->A,b->A,str);
899: MatCopy(a->B,b->B,str);
900: }
901: return(0);
902: }
904: PetscErrorCode MatSetUp_MPISELL(Mat A)
905: {
909: MatMPISELLSetPreallocation(A,PETSC_DEFAULT,0,PETSC_DEFAULT,0);
910: return(0);
911: }
914: extern PetscErrorCode MatConjugate_SeqSELL(Mat);
916: PetscErrorCode MatConjugate_MPISELL(Mat mat)
917: {
918: #if defined(PETSC_USE_COMPLEX)
920: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
923: MatConjugate_SeqSELL(sell->A);
924: MatConjugate_SeqSELL(sell->B);
925: #else
927: #endif
928: return(0);
929: }
931: PetscErrorCode MatRealPart_MPISELL(Mat A)
932: {
933: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
937: MatRealPart(a->A);
938: MatRealPart(a->B);
939: return(0);
940: }
942: PetscErrorCode MatImaginaryPart_MPISELL(Mat A)
943: {
944: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
948: MatImaginaryPart(a->A);
949: MatImaginaryPart(a->B);
950: return(0);
951: }
953: PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values)
954: {
955: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
959: MatInvertBlockDiagonal(a->A,values);
960: A->factorerrortype = a->A->factorerrortype;
961: return(0);
962: }
964: static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx)
965: {
967: Mat_MPISELL *sell=(Mat_MPISELL*)x->data;
970: MatSetRandom(sell->A,rctx);
971: MatSetRandom(sell->B,rctx);
972: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
973: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
974: return(0);
975: }
977: PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A)
978: {
982: PetscOptionsHead(PetscOptionsObject,"MPISELL options");
983: PetscObjectOptionsBegin((PetscObject)A);
984: PetscOptionsEnd();
985: return(0);
986: }
988: PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a)
989: {
991: Mat_MPISELL *msell=(Mat_MPISELL*)Y->data;
992: Mat_SeqSELL *sell=(Mat_SeqSELL*)msell->A->data;
995: if (!Y->preallocated) {
996: MatMPISELLSetPreallocation(Y,1,NULL,0,NULL);
997: } else if (!sell->nz) {
998: PetscInt nonew = sell->nonew;
999: MatSeqSELLSetPreallocation(msell->A,1,NULL);
1000: sell->nonew = nonew;
1001: }
1002: MatShift_Basic(Y,a);
1003: return(0);
1004: }
1006: PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool *missing,PetscInt *d)
1007: {
1008: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1012: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices");
1013: MatMissingDiagonal(a->A,missing,d);
1014: if (d) {
1015: PetscInt rstart;
1016: MatGetOwnershipRange(A,&rstart,NULL);
1017: *d += rstart;
1019: }
1020: return(0);
1021: }
1023: PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a)
1024: {
1026: *a = ((Mat_MPISELL*)A->data)->A;
1027: return(0);
1028: }
1030: /* -------------------------------------------------------------------*/
1031: static struct _MatOps MatOps_Values = {MatSetValues_MPISELL,
1032: 0,
1033: 0,
1034: MatMult_MPISELL,
1035: /* 4*/ MatMultAdd_MPISELL,
1036: MatMultTranspose_MPISELL,
1037: MatMultTransposeAdd_MPISELL,
1038: 0,
1039: 0,
1040: 0,
1041: /*10*/ 0,
1042: 0,
1043: 0,
1044: MatSOR_MPISELL,
1045: 0,
1046: /*15*/ MatGetInfo_MPISELL,
1047: MatEqual_MPISELL,
1048: MatGetDiagonal_MPISELL,
1049: MatDiagonalScale_MPISELL,
1050: 0,
1051: /*20*/ MatAssemblyBegin_MPISELL,
1052: MatAssemblyEnd_MPISELL,
1053: MatSetOption_MPISELL,
1054: MatZeroEntries_MPISELL,
1055: /*24*/ 0,
1056: 0,
1057: 0,
1058: 0,
1059: 0,
1060: /*29*/ MatSetUp_MPISELL,
1061: 0,
1062: 0,
1063: MatGetDiagonalBlock_MPISELL,
1064: 0,
1065: /*34*/ MatDuplicate_MPISELL,
1066: 0,
1067: 0,
1068: 0,
1069: 0,
1070: /*39*/ 0,
1071: 0,
1072: 0,
1073: MatGetValues_MPISELL,
1074: MatCopy_MPISELL,
1075: /*44*/ 0,
1076: MatScale_MPISELL,
1077: MatShift_MPISELL,
1078: MatDiagonalSet_MPISELL,
1079: 0,
1080: /*49*/ MatSetRandom_MPISELL,
1081: 0,
1082: 0,
1083: 0,
1084: 0,
1085: /*54*/ MatFDColoringCreate_MPIXAIJ,
1086: 0,
1087: MatSetUnfactored_MPISELL,
1088: 0,
1089: 0,
1090: /*59*/ 0,
1091: MatDestroy_MPISELL,
1092: MatView_MPISELL,
1093: 0,
1094: 0,
1095: /*64*/ 0,
1096: 0,
1097: 0,
1098: 0,
1099: 0,
1100: /*69*/ 0,
1101: 0,
1102: 0,
1103: 0,
1104: 0,
1105: 0,
1106: /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */
1107: MatSetFromOptions_MPISELL,
1108: 0,
1109: 0,
1110: 0,
1111: /*80*/ 0,
1112: 0,
1113: 0,
1114: /*83*/ 0,
1115: 0,
1116: 0,
1117: 0,
1118: 0,
1119: 0,
1120: /*89*/ 0,
1121: 0,
1122: 0,
1123: 0,
1124: 0,
1125: /*94*/ 0,
1126: 0,
1127: 0,
1128: 0,
1129: 0,
1130: /*99*/ 0,
1131: 0,
1132: 0,
1133: MatConjugate_MPISELL,
1134: 0,
1135: /*104*/0,
1136: MatRealPart_MPISELL,
1137: MatImaginaryPart_MPISELL,
1138: 0,
1139: 0,
1140: /*109*/0,
1141: 0,
1142: 0,
1143: 0,
1144: MatMissingDiagonal_MPISELL,
1145: /*114*/0,
1146: 0,
1147: MatGetGhosts_MPISELL,
1148: 0,
1149: 0,
1150: /*119*/0,
1151: 0,
1152: 0,
1153: 0,
1154: 0,
1155: /*124*/0,
1156: 0,
1157: MatInvertBlockDiagonal_MPISELL,
1158: 0,
1159: 0,
1160: /*129*/0,
1161: 0,
1162: 0,
1163: 0,
1164: 0,
1165: /*134*/0,
1166: 0,
1167: 0,
1168: 0,
1169: 0,
1170: /*139*/0,
1171: 0,
1172: 0,
1173: MatFDColoringSetUp_MPIXAIJ,
1174: 0,
1175: /*144*/0
1176: };
1178: /* ----------------------------------------------------------------------------------------*/
1180: PetscErrorCode MatStoreValues_MPISELL(Mat mat)
1181: {
1182: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
1186: MatStoreValues(sell->A);
1187: MatStoreValues(sell->B);
1188: return(0);
1189: }
1191: PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
1192: {
1193: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
1197: MatRetrieveValues(sell->A);
1198: MatRetrieveValues(sell->B);
1199: return(0);
1200: }
1202: PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[])
1203: {
1204: Mat_MPISELL *b;
1208: PetscLayoutSetUp(B->rmap);
1209: PetscLayoutSetUp(B->cmap);
1210: b = (Mat_MPISELL*)B->data;
1212: if (!B->preallocated) {
1213: /* Explicitly create 2 MATSEQSELL matrices. */
1214: MatCreate(PETSC_COMM_SELF,&b->A);
1215: MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
1216: MatSetBlockSizesFromMats(b->A,B,B);
1217: MatSetType(b->A,MATSEQSELL);
1218: PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
1219: MatCreate(PETSC_COMM_SELF,&b->B);
1220: MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
1221: MatSetBlockSizesFromMats(b->B,B,B);
1222: MatSetType(b->B,MATSEQSELL);
1223: PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
1224: }
1226: MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);
1227: MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);
1228: B->preallocated = PETSC_TRUE;
1229: B->was_assembled = PETSC_FALSE;
1230: /*
1231: critical for MatAssemblyEnd to work.
1232: MatAssemblyBegin checks it to set up was_assembled
1233: and MatAssemblyEnd checks was_assembled to determine whether to build garray
1234: */
1235: B->assembled = PETSC_FALSE;
1236: return(0);
1237: }
1239: PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat)
1240: {
1241: Mat mat;
1242: Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data;
1246: *newmat = 0;
1247: MatCreate(PetscObjectComm((PetscObject)matin),&mat);
1248: MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);
1249: MatSetBlockSizesFromMats(mat,matin,matin);
1250: MatSetType(mat,((PetscObject)matin)->type_name);
1251: PetscMemcpy(mat->ops,matin->ops,sizeof(struct _MatOps));
1252: a = (Mat_MPISELL*)mat->data;
1254: mat->factortype = matin->factortype;
1255: mat->assembled = PETSC_TRUE;
1256: mat->insertmode = NOT_SET_VALUES;
1257: mat->preallocated = PETSC_TRUE;
1259: a->size = oldmat->size;
1260: a->rank = oldmat->rank;
1261: a->donotstash = oldmat->donotstash;
1262: a->roworiented = oldmat->roworiented;
1263: a->rowindices = 0;
1264: a->rowvalues = 0;
1265: a->getrowactive = PETSC_FALSE;
1267: PetscLayoutReference(matin->rmap,&mat->rmap);
1268: PetscLayoutReference(matin->cmap,&mat->cmap);
1270: if (oldmat->colmap) {
1271: #if defined(PETSC_USE_CTABLE)
1272: PetscTableCreateCopy(oldmat->colmap,&a->colmap);
1273: #else
1274: PetscMalloc1(mat->cmap->N,&a->colmap);
1275: PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));
1276: PetscMemcpy(a->colmap,oldmat->colmap,(mat->cmap->N)*sizeof(PetscInt));
1277: #endif
1278: } else a->colmap = 0;
1279: if (oldmat->garray) {
1280: PetscInt len;
1281: len = oldmat->B->cmap->n;
1282: PetscMalloc1(len+1,&a->garray);
1283: PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));
1284: if (len) { PetscMemcpy(a->garray,oldmat->garray,len*sizeof(PetscInt)); }
1285: } else a->garray = 0;
1287: VecDuplicate(oldmat->lvec,&a->lvec);
1288: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);
1289: VecScatterCopy(oldmat->Mvctx,&a->Mvctx);
1290: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);
1291: MatDuplicate(oldmat->A,cpvalues,&a->A);
1292: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);
1293: MatDuplicate(oldmat->B,cpvalues,&a->B);
1294: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);
1295: PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);
1296: *newmat = mat;
1297: return(0);
1298: }
1300: /*@C
1301: MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format.
1302: For good matrix assembly performance the user should preallocate the matrix storage by
1303: setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz).
1305: Collective on MPI_Comm
1307: Input Parameters:
1308: + B - the matrix
1309: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1310: (same value is used for all local rows)
1311: . d_nnz - array containing the number of nonzeros in the various rows of the
1312: DIAGONAL portion of the local submatrix (possibly different for each row)
1313: or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1314: The size of this array is equal to the number of local rows, i.e 'm'.
1315: For matrices that will be factored, you must leave room for (and set)
1316: the diagonal entry even if it is zero.
1317: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1318: submatrix (same value is used for all local rows).
1319: - o_nnz - array containing the number of nonzeros in the various rows of the
1320: OFF-DIAGONAL portion of the local submatrix (possibly different for
1321: each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1322: structure. The size of this array is equal to the number
1323: of local rows, i.e 'm'.
1325: If the *_nnz parameter is given then the *_nz parameter is ignored
1327: The stored row and column indices begin with zero.
1329: The parallel matrix is partitioned such that the first m0 rows belong to
1330: process 0, the next m1 rows belong to process 1, the next m2 rows belong
1331: to process 2 etc.. where m0,m1,m2... are the input parameter 'm'.
1333: The DIAGONAL portion of the local submatrix of a processor can be defined
1334: as the submatrix which is obtained by extraction the part corresponding to
1335: the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
1336: first row that belongs to the processor, r2 is the last row belonging to
1337: the this processor, and c1-c2 is range of indices of the local part of a
1338: vector suitable for applying the matrix to. This is an mxn matrix. In the
1339: common case of a square matrix, the row and column ranges are the same and
1340: the DIAGONAL part is also square. The remaining portion of the local
1341: submatrix (mxN) constitute the OFF-DIAGONAL portion.
1343: If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored.
1345: You can call MatGetInfo() to get information on how effective the preallocation was;
1346: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
1347: You can also run with the option -info and look for messages with the string
1348: malloc in them to see if additional memory allocation was needed.
1350: Example usage:
1352: Consider the following 8x8 matrix with 34 non-zero values, that is
1353: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1354: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1355: as follows:
1357: .vb
1358: 1 2 0 | 0 3 0 | 0 4
1359: Proc0 0 5 6 | 7 0 0 | 8 0
1360: 9 0 10 | 11 0 0 | 12 0
1361: -------------------------------------
1362: 13 0 14 | 15 16 17 | 0 0
1363: Proc1 0 18 0 | 19 20 21 | 0 0
1364: 0 0 0 | 22 23 0 | 24 0
1365: -------------------------------------
1366: Proc2 25 26 27 | 0 0 28 | 29 0
1367: 30 0 0 | 31 32 33 | 0 34
1368: .ve
1370: This can be represented as a collection of submatrices as:
1372: .vb
1373: A B C
1374: D E F
1375: G H I
1376: .ve
1378: Where the submatrices A,B,C are owned by proc0, D,E,F are
1379: owned by proc1, G,H,I are owned by proc2.
1381: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1382: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1383: The 'M','N' parameters are 8,8, and have the same values on all procs.
1385: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1386: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1387: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1388: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1389: part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1390: matrix, ans [DF] as another SeqSELL matrix.
1392: When d_nz, o_nz parameters are specified, d_nz storage elements are
1393: allocated for every row of the local diagonal submatrix, and o_nz
1394: storage locations are allocated for every row of the OFF-DIAGONAL submat.
1395: One way to choose d_nz and o_nz is to use the max nonzerors per local
1396: rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1397: In this case, the values of d_nz,o_nz are:
1398: .vb
1399: proc0 : dnz = 2, o_nz = 2
1400: proc1 : dnz = 3, o_nz = 2
1401: proc2 : dnz = 1, o_nz = 4
1402: .ve
1403: We are allocating m*(d_nz+o_nz) storage locations for every proc. This
1404: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1405: for proc3. i.e we are using 12+15+10=37 storage locations to store
1406: 34 values.
1408: When d_nnz, o_nnz parameters are specified, the storage is specified
1409: for every row, coresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1410: In the above case the values for d_nnz,o_nnz are:
1411: .vb
1412: proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
1413: proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
1414: proc2: d_nnz = [1,1] and o_nnz = [4,4]
1415: .ve
1416: Here the space allocated is according to nz (or maximum values in the nnz
1417: if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37
1419: Level: intermediate
1421: .keywords: matrix, sell, sparse, parallel
1423: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(),
1424: MATMPISELL, MatGetInfo(), PetscSplitOwnership()
1425: @*/
1426: PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1427: {
1433: PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1434: return(0);
1435: }
1437: /*@C
1438: MatCreateSELL - Creates a sparse parallel matrix in SELL format.
1440: Collective on MPI_Comm
1442: Input Parameters:
1443: + comm - MPI communicator
1444: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
1445: This value should be the same as the local size used in creating the
1446: y vector for the matrix-vector product y = Ax.
1447: . n - This value should be the same as the local size used in creating the
1448: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
1449: calculated if N is given) For square matrices n is almost always m.
1450: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
1451: . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
1452: . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix
1453: (same value is used for all local rows)
1454: . d_rlen - array containing the number of nonzeros in the various rows of the
1455: DIAGONAL portion of the local submatrix (possibly different for each row)
1456: or NULL, if d_rlenmax is used to specify the nonzero structure.
1457: The size of this array is equal to the number of local rows, i.e 'm'.
1458: . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local
1459: submatrix (same value is used for all local rows).
1460: - o_rlen - array containing the number of nonzeros in the various rows of the
1461: OFF-DIAGONAL portion of the local submatrix (possibly different for
1462: each row) or NULL, if o_rlenmax is used to specify the nonzero
1463: structure. The size of this array is equal to the number
1464: of local rows, i.e 'm'.
1466: Output Parameter:
1467: . A - the matrix
1469: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
1470: MatXXXXSetPreallocation() paradgm instead of this routine directly.
1471: [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation]
1473: Notes:
1474: If the *_rlen parameter is given then the *_rlenmax parameter is ignored
1476: m,n,M,N parameters specify the size of the matrix, and its partitioning across
1477: processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate
1478: storage requirements for this matrix.
1480: If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one
1481: processor than it must be used on all processors that share the object for
1482: that argument.
1484: The user MUST specify either the local or global matrix dimensions
1485: (possibly both).
1487: The parallel matrix is partitioned across processors such that the
1488: first m0 rows belong to process 0, the next m1 rows belong to
1489: process 1, the next m2 rows belong to process 2 etc.. where
1490: m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores
1491: values corresponding to [m x N] submatrix.
1493: The columns are logically partitioned with the n0 columns belonging
1494: to 0th partition, the next n1 columns belonging to the next
1495: partition etc.. where n0,n1,n2... are the input parameter 'n'.
1497: The DIAGONAL portion of the local submatrix on any given processor
1498: is the submatrix corresponding to the rows and columns m,n
1499: corresponding to the given processor. i.e diagonal matrix on
1500: process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
1501: etc. The remaining portion of the local submatrix [m x (N-n)]
1502: constitute the OFF-DIAGONAL portion. The example below better
1503: illustrates this concept.
1505: For a square global matrix we define each processor's diagonal portion
1506: to be its local rows and the corresponding columns (a square submatrix);
1507: each processor's off-diagonal portion encompasses the remainder of the
1508: local matrix (a rectangular submatrix).
1510: If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored.
1512: When calling this routine with a single process communicator, a matrix of
1513: type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this
1514: type of communicator, use the construction mechanism:
1515: MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...);
1517: Options Database Keys:
1518: - -mat_sell_oneindex - Internally use indexing starting at 1
1519: rather than 0. Note that when calling MatSetValues(),
1520: the user still MUST index entries starting at 0!
1523: Example usage:
1525: Consider the following 8x8 matrix with 34 non-zero values, that is
1526: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1527: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1528: as follows:
1530: .vb
1531: 1 2 0 | 0 3 0 | 0 4
1532: Proc0 0 5 6 | 7 0 0 | 8 0
1533: 9 0 10 | 11 0 0 | 12 0
1534: -------------------------------------
1535: 13 0 14 | 15 16 17 | 0 0
1536: Proc1 0 18 0 | 19 20 21 | 0 0
1537: 0 0 0 | 22 23 0 | 24 0
1538: -------------------------------------
1539: Proc2 25 26 27 | 0 0 28 | 29 0
1540: 30 0 0 | 31 32 33 | 0 34
1541: .ve
1543: This can be represented as a collection of submatrices as:
1545: .vb
1546: A B C
1547: D E F
1548: G H I
1549: .ve
1551: Where the submatrices A,B,C are owned by proc0, D,E,F are
1552: owned by proc1, G,H,I are owned by proc2.
1554: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1555: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1556: The 'M','N' parameters are 8,8, and have the same values on all procs.
1558: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1559: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1560: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1561: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1562: part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1563: matrix, ans [DF] as another SeqSELL matrix.
1565: When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are
1566: allocated for every row of the local diagonal submatrix, and o_rlenmax
1567: storage locations are allocated for every row of the OFF-DIAGONAL submat.
1568: One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local
1569: rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1570: In this case, the values of d_rlenmax,o_rlenmax are:
1571: .vb
1572: proc0 : d_rlenmax = 2, o_rlenmax = 2
1573: proc1 : d_rlenmax = 3, o_rlenmax = 2
1574: proc2 : d_rlenmax = 1, o_rlenmax = 4
1575: .ve
1576: We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This
1577: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1578: for proc3. i.e we are using 12+15+10=37 storage locations to store
1579: 34 values.
1581: When d_rlen, o_rlen parameters are specified, the storage is specified
1582: for every row, coresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1583: In the above case the values for d_nnz,o_nnz are:
1584: .vb
1585: proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
1586: proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
1587: proc2: d_nnz = [1,1] and o_nnz = [4,4]
1588: .ve
1589: Here the space allocated is still 37 though there are 34 nonzeros because
1590: the allocation is always done according to rlenmax.
1592: Level: intermediate
1594: .keywords: matrix, sell, sparse, parallel
1596: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(),
1597: MATMPISELL, MatCreateMPISELLWithArrays()
1598: @*/
1599: 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)
1600: {
1602: PetscMPIInt size;
1605: MatCreate(comm,A);
1606: MatSetSizes(*A,m,n,M,N);
1607: MPI_Comm_size(comm,&size);
1608: if (size > 1) {
1609: MatSetType(*A,MATMPISELL);
1610: MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);
1611: } else {
1612: MatSetType(*A,MATSEQSELL);
1613: MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);
1614: }
1615: return(0);
1616: }
1618: PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[])
1619: {
1620: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1621: PetscBool flg;
1625: PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);
1626: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input");
1627: if (Ad) *Ad = a->A;
1628: if (Ao) *Ao = a->B;
1629: if (colmap) *colmap = a->garray;
1630: return(0);
1631: }
1633: /*@C
1634: MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns
1636: Not Collective
1638: Input Parameters:
1639: + A - the matrix
1640: . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
1641: - row, col - index sets of rows and columns to extract (or NULL)
1643: Output Parameter:
1644: . A_loc - the local sequential matrix generated
1646: Level: developer
1648: .seealso: MatGetOwnershipRange(), MatMPISELLGetLocalMat()
1650: @*/
1651: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc)
1652: {
1653: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1655: PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx;
1656: IS isrowa,iscola;
1657: Mat *aloc;
1658: PetscBool match;
1661: PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match);
1662: if (!match) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input");
1663: PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0);
1664: if (!row) {
1665: start = A->rmap->rstart; end = A->rmap->rend;
1666: ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa);
1667: } else {
1668: isrowa = *row;
1669: }
1670: if (!col) {
1671: start = A->cmap->rstart;
1672: cmap = a->garray;
1673: nzA = a->A->cmap->n;
1674: nzB = a->B->cmap->n;
1675: PetscMalloc1(nzA+nzB, &idx);
1676: ncols = 0;
1677: for (i=0; i<nzB; i++) {
1678: if (cmap[i] < start) idx[ncols++] = cmap[i];
1679: else break;
1680: }
1681: imark = i;
1682: for (i=0; i<nzA; i++) idx[ncols++] = start + i;
1683: for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i];
1684: ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola);
1685: } else {
1686: iscola = *col;
1687: }
1688: if (scall != MAT_INITIAL_MATRIX) {
1689: PetscMalloc1(1,&aloc);
1690: aloc[0] = *A_loc;
1691: }
1692: MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc);
1693: *A_loc = aloc[0];
1694: PetscFree(aloc);
1695: if (!row) {
1696: ISDestroy(&isrowa);
1697: }
1698: if (!col) {
1699: ISDestroy(&iscola);
1700: }
1701: PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0);
1702: return(0);
1703: }
1705: #include <../src/mat/impls/aij/mpi/mpiaij.h>
1707: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1708: {
1710: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1711: Mat B;
1712: Mat_MPIAIJ *b;
1715: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled");
1717: if (reuse == MAT_REUSE_MATRIX) {
1718: B = *newmat;
1719: } else {
1720: MatCreate(PetscObjectComm((PetscObject)A),&B);
1721: MatSetType(B,MATMPIAIJ);
1722: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1723: MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1724: MatSeqAIJSetPreallocation(B,0,NULL);
1725: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1726: }
1727: b = (Mat_MPIAIJ*) B->data;
1729: if (reuse == MAT_REUSE_MATRIX) {
1730: MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);
1731: MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);
1732: } else {
1733: MatDestroy(&b->A);
1734: MatDestroy(&b->B);
1735: MatDisAssemble_MPISELL(A);
1736: MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);
1737: MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);
1738: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1739: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1740: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1741: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1742: }
1744: if (reuse == MAT_INPLACE_MATRIX) {
1745: MatHeaderReplace(A,&B);
1746: } else {
1747: *newmat = B;
1748: }
1749: return(0);
1750: }
1752: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1753: {
1755: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data;
1756: Mat B;
1757: Mat_MPISELL *b;
1760: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled");
1762: if (reuse == MAT_REUSE_MATRIX) {
1763: B = *newmat;
1764: } else {
1765: MatCreate(PetscObjectComm((PetscObject)A),&B);
1766: MatSetType(B,MATMPISELL);
1767: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1768: MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1769: MatSeqAIJSetPreallocation(B,0,NULL);
1770: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1771: }
1772: b = (Mat_MPISELL*) B->data;
1774: if (reuse == MAT_REUSE_MATRIX) {
1775: MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);
1776: MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);
1777: } else {
1778: MatDestroy(&b->A);
1779: MatDestroy(&b->B);
1780: MatDisAssemble_MPIAIJ(A);
1781: MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);
1782: MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);
1783: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1784: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1785: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1786: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1787: }
1789: if (reuse == MAT_INPLACE_MATRIX) {
1790: MatHeaderReplace(A,&B);
1791: } else {
1792: *newmat = B;
1793: }
1794: return(0);
1795: }
1797: PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1798: {
1799: Mat_MPISELL *mat=(Mat_MPISELL*)matin->data;
1801: Vec bb1=0;
1804: if (flag == SOR_APPLY_UPPER) {
1805: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1806: return(0);
1807: }
1809: if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) {
1810: VecDuplicate(bb,&bb1);
1811: }
1813: if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1814: if (flag & SOR_ZERO_INITIAL_GUESS) {
1815: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1816: its--;
1817: }
1819: while (its--) {
1820: VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1821: VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1823: /* update rhs: bb1 = bb - B*x */
1824: VecScale(mat->lvec,-1.0);
1825: (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);
1827: /* local sweep */
1828: (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);
1829: }
1830: } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
1831: if (flag & SOR_ZERO_INITIAL_GUESS) {
1832: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1833: its--;
1834: }
1835: while (its--) {
1836: VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1837: VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1839: /* update rhs: bb1 = bb - B*x */
1840: VecScale(mat->lvec,-1.0);
1841: (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);
1843: /* local sweep */
1844: (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);
1845: }
1846: } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1847: if (flag & SOR_ZERO_INITIAL_GUESS) {
1848: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1849: its--;
1850: }
1851: while (its--) {
1852: VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1853: VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1855: /* update rhs: bb1 = bb - B*x */
1856: VecScale(mat->lvec,-1.0);
1857: (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);
1859: /* local sweep */
1860: (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx);
1861: }
1862: } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported");
1864: VecDestroy(&bb1);
1866: matin->factorerrortype = mat->A->factorerrortype;
1867: return(0);
1868: }
1870: /*MC
1871: MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices.
1873: Options Database Keys:
1874: . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions()
1876: Level: beginner
1878: .seealso: MatCreateSELL()
1879: M*/
1880: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1881: {
1882: Mat_MPISELL *b;
1884: PetscMPIInt size;
1887: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1888: PetscNewLog(B,&b);
1889: B->data = (void*)b;
1890: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1891: B->assembled = PETSC_FALSE;
1892: B->insertmode = NOT_SET_VALUES;
1893: b->size = size;
1894: MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);
1895: /* build cache for off array entries formed */
1896: MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);
1898: b->donotstash = PETSC_FALSE;
1899: b->colmap = 0;
1900: b->garray = 0;
1901: b->roworiented = PETSC_TRUE;
1903: /* stuff used for matrix vector multiply */
1904: b->lvec = NULL;
1905: b->Mvctx = NULL;
1907: /* stuff for MatGetRow() */
1908: b->rowindices = 0;
1909: b->rowvalues = 0;
1910: b->getrowactive = PETSC_FALSE;
1912: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);
1913: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);
1914: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);
1915: PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);
1916: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);
1917: PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);
1918: PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);
1919: return(0);
1920: }