Actual source code: mpisell.c
1: #include <../src/mat/impls/aij/mpi/mpiaij.h>
2: #include <../src/mat/impls/sell/mpi/mpisell.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petsc/private/isimpl.h>
5: #include <petscblaslapack.h>
6: #include <petscsf.h>
8: /*MC
9: MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices.
11: This matrix type is identical to MATSEQSELL when constructed with a single process communicator,
12: and MATMPISELL otherwise. As a result, for single process communicators,
13: MatSeqSELLSetPreallocation is supported, and similarly MatMPISELLSetPreallocation is supported
14: for communicators controlling multiple processes. It is recommended that you call both of
15: the above preallocation routines for simplicity.
17: Options Database Keys:
18: . -mat_type sell - sets the matrix type to "sell" during a call to MatSetFromOptions()
20: Level: beginner
22: .seealso: MatCreateSELL(), MatCreateSeqSELL(), MATSEQSELL, MATMPISELL
23: M*/
25: PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is)
26: {
28: Mat_MPISELL *sell=(Mat_MPISELL*)Y->data;
31: if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) {
32: MatDiagonalSet(sell->A,D,is);
33: } else {
34: MatDiagonalSet_Default(Y,D,is);
35: }
36: return(0);
37: }
39: /*
40: Local utility routine that creates a mapping from the global column
41: number to the local number in the off-diagonal part of the local
42: storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at
43: a slightly higher hash table cost; without it it is not scalable (each processor
44: has an order N integer array but is fast to acess.
45: */
46: PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat)
47: {
48: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
50: PetscInt n=sell->B->cmap->n,i;
53: if (!sell->garray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MPISELL Matrix was assembled but is missing garray");
54: #if defined(PETSC_USE_CTABLE)
55: PetscTableCreate(n,mat->cmap->N+1,&sell->colmap);
56: for (i=0; i<n; i++) {
57: PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES);
58: }
59: #else
60: PetscCalloc1(mat->cmap->N+1,&sell->colmap);
61: PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt));
62: for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1;
63: #endif
64: return(0);
65: }
67: #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \
68: { \
69: if (col <= lastcol1) low1 = 0; \
70: else high1 = nrow1; \
71: lastcol1 = col; \
72: while (high1-low1 > 5) { \
73: t = (low1+high1)/2; \
74: if (*(cp1+8*t) > col) high1 = t; \
75: else low1 = t; \
76: } \
77: for (_i=low1; _i<high1; _i++) { \
78: if (*(cp1+8*_i) > col) break; \
79: if (*(cp1+8*_i) == col) { \
80: if (addv == ADD_VALUES) *(vp1+8*_i) += value; \
81: else *(vp1+8*_i) = value; \
82: goto a_noinsert; \
83: } \
84: } \
85: if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \
86: if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \
87: 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); \
88: MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \
89: /* shift up all the later entries in this row */ \
90: for (ii=nrow1-1; ii>=_i; ii--) { \
91: *(cp1+8*(ii+1)) = *(cp1+8*ii); \
92: *(vp1+8*(ii+1)) = *(vp1+8*ii); \
93: } \
94: *(cp1+8*_i) = col; \
95: *(vp1+8*_i) = value; \
96: a->nz++; nrow1++; A->nonzerostate++; \
97: a_noinsert: ; \
98: a->rlen[row] = nrow1; \
99: }
101: #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \
102: { \
103: if (col <= lastcol2) low2 = 0; \
104: else high2 = nrow2; \
105: lastcol2 = col; \
106: while (high2-low2 > 5) { \
107: t = (low2+high2)/2; \
108: if (*(cp2+8*t) > col) high2 = t; \
109: else low2 = t; \
110: } \
111: for (_i=low2; _i<high2; _i++) { \
112: if (*(cp2+8*_i) > col) break; \
113: if (*(cp2+8*_i) == col) { \
114: if (addv == ADD_VALUES) *(vp2+8*_i) += value; \
115: else *(vp2+8*_i) = value; \
116: goto b_noinsert; \
117: } \
118: } \
119: if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
120: if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \
121: 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); \
122: MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \
123: /* shift up all the later entries in this row */ \
124: for (ii=nrow2-1; ii>=_i; ii--) { \
125: *(cp2+8*(ii+1)) = *(cp2+8*ii); \
126: *(vp2+8*(ii+1)) = *(vp2+8*ii); \
127: } \
128: *(cp2+8*_i) = col; \
129: *(vp2+8*_i) = value; \
130: b->nz++; nrow2++; B->nonzerostate++; \
131: b_noinsert: ; \
132: b->rlen[row] = nrow2; \
133: }
135: PetscErrorCode MatSetValues_MPISELL(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
136: {
137: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
138: PetscScalar value;
140: PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend,shift1,shift2;
141: PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;
142: PetscBool roworiented=sell->roworiented;
144: /* Some Variables required in the macro */
145: Mat A=sell->A;
146: Mat_SeqSELL *a=(Mat_SeqSELL*)A->data;
147: PetscBool ignorezeroentries=a->ignorezeroentries,found;
148: Mat B=sell->B;
149: Mat_SeqSELL *b=(Mat_SeqSELL*)B->data;
150: PetscInt *cp1,*cp2,ii,_i,nrow1,nrow2,low1,high1,low2,high2,t,lastcol1,lastcol2;
151: MatScalar *vp1,*vp2;
154: for (i=0; i<m; i++) {
155: if (im[i] < 0) continue;
156: if (PetscUnlikelyDebug(im[i] >= mat->rmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap->N-1);
157: if (im[i] >= rstart && im[i] < rend) {
158: row = im[i] - rstart;
159: lastcol1 = -1;
160: shift1 = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */
161: cp1 = a->colidx+shift1;
162: vp1 = a->val+shift1;
163: nrow1 = a->rlen[row];
164: low1 = 0;
165: high1 = nrow1;
166: lastcol2 = -1;
167: shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
168: cp2 = b->colidx+shift2;
169: vp2 = b->val+shift2;
170: nrow2 = b->rlen[row];
171: low2 = 0;
172: high2 = nrow2;
174: for (j=0; j<n; j++) {
175: if (roworiented) value = v[i*n+j];
176: else value = v[i+j*m];
177: if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue;
178: if (in[j] >= cstart && in[j] < cend) {
179: col = in[j] - cstart;
180: MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */
181: } else if (in[j] < 0) continue;
182: else if (PetscUnlikelyDebug(in[j] >= mat->cmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[j],mat->cmap->N-1);
183: else {
184: if (mat->was_assembled) {
185: if (!sell->colmap) {
186: MatCreateColmap_MPISELL_Private(mat);
187: }
188: #if defined(PETSC_USE_CTABLE)
189: PetscTableFind(sell->colmap,in[j]+1,&col);
190: col--;
191: #else
192: col = sell->colmap[in[j]] - 1;
193: #endif
194: if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) {
195: MatDisAssemble_MPISELL(mat);
196: col = in[j];
197: /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */
198: B = sell->B;
199: b = (Mat_SeqSELL*)B->data;
200: shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */
201: cp2 = b->colidx+shift2;
202: vp2 = b->val+shift2;
203: nrow2 = b->rlen[row];
204: low2 = 0;
205: high2 = nrow2;
206: } 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]);
207: } else col = in[j];
208: MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */
209: }
210: }
211: } else {
212: 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]);
213: if (!sell->donotstash) {
214: mat->assembled = PETSC_FALSE;
215: if (roworiented) {
216: MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
217: } else {
218: MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));
219: }
220: }
221: }
222: }
223: return(0);
224: }
226: PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
227: {
228: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
230: PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend;
231: PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col;
234: for (i=0; i<m; i++) {
235: if (idxm[i] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",idxm[i]);*/
236: 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);
237: if (idxm[i] >= rstart && idxm[i] < rend) {
238: row = idxm[i] - rstart;
239: for (j=0; j<n; j++) {
240: if (idxn[j] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",idxn[j]); */
241: 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);
242: if (idxn[j] >= cstart && idxn[j] < cend) {
243: col = idxn[j] - cstart;
244: MatGetValues(sell->A,1,&row,1,&col,v+i*n+j);
245: } else {
246: if (!sell->colmap) {
247: MatCreateColmap_MPISELL_Private(mat);
248: }
249: #if defined(PETSC_USE_CTABLE)
250: PetscTableFind(sell->colmap,idxn[j]+1,&col);
251: col--;
252: #else
253: col = sell->colmap[idxn[j]] - 1;
254: #endif
255: if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0;
256: else {
257: MatGetValues(sell->B,1,&row,1,&col,v+i*n+j);
258: }
259: }
260: }
261: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported");
262: }
263: return(0);
264: }
266: extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec);
268: PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode)
269: {
270: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
272: PetscInt nstash,reallocs;
275: if (sell->donotstash || mat->nooffprocentries) return(0);
277: MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
278: MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);
279: PetscInfo2(sell->A,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);
280: return(0);
281: }
283: PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode)
284: {
285: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
287: PetscMPIInt n;
288: PetscInt i,flg;
289: PetscInt *row,*col;
290: PetscScalar *val;
291: PetscBool other_disassembled;
292: /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */
294: if (!sell->donotstash && !mat->nooffprocentries) {
295: while (1) {
296: MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);
297: if (!flg) break;
299: for (i=0; i<n; i++) { /* assemble one by one */
300: MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode);
301: }
302: }
303: MatStashScatterEnd_Private(&mat->stash);
304: }
305: MatAssemblyBegin(sell->A,mode);
306: MatAssemblyEnd(sell->A,mode);
308: /*
309: determine if any processor has disassembled, if so we must
310: also disassemble ourselfs, in order that we may reassemble.
311: */
312: /*
313: if nonzero structure of submatrix B cannot change then we know that
314: no processor disassembled thus we can skip this stuff
315: */
316: if (!((Mat_SeqSELL*)sell->B->data)->nonew) {
317: MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat));
318: if (mat->was_assembled && !other_disassembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet\n");
319: }
320: if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) {
321: MatSetUpMultiply_MPISELL(mat);
322: }
323: /*
324: MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE);
325: */
326: MatAssemblyBegin(sell->B,mode);
327: MatAssemblyEnd(sell->B,mode);
328: PetscFree2(sell->rowvalues,sell->rowindices);
329: sell->rowvalues = NULL;
330: VecDestroy(&sell->diag);
332: /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
333: if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) {
334: PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate;
335: MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));
336: }
337: return(0);
338: }
340: PetscErrorCode MatZeroEntries_MPISELL(Mat A)
341: {
342: Mat_MPISELL *l=(Mat_MPISELL*)A->data;
346: MatZeroEntries(l->A);
347: MatZeroEntries(l->B);
348: return(0);
349: }
351: PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy)
352: {
353: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
355: PetscInt nt;
358: VecGetLocalSize(xx,&nt);
359: 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);
360: VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
361: (*a->A->ops->mult)(a->A,xx,yy);
362: VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
363: (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);
364: return(0);
365: }
367: PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx)
368: {
369: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
373: MatMultDiagonalBlock(a->A,bb,xx);
374: return(0);
375: }
377: PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
378: {
379: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
383: VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
384: (*a->A->ops->multadd)(a->A,xx,yy,zz);
385: VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
386: (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
387: return(0);
388: }
390: PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy)
391: {
392: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
396: /* do nondiagonal part */
397: (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
398: /* do local part */
399: (*a->A->ops->multtranspose)(a->A,xx,yy);
400: /* add partial results together */
401: VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
402: VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
403: return(0);
404: }
406: PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f)
407: {
408: MPI_Comm comm;
409: Mat_MPISELL *Asell=(Mat_MPISELL*)Amat->data,*Bsell;
410: Mat Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs;
411: IS Me,Notme;
413: PetscInt M,N,first,last,*notme,i;
414: PetscMPIInt size;
417: /* Easy test: symmetric diagonal block */
418: Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A;
419: MatIsTranspose(Adia,Bdia,tol,f);
420: if (!*f) return(0);
421: PetscObjectGetComm((PetscObject)Amat,&comm);
422: MPI_Comm_size(comm,&size);
423: if (size == 1) return(0);
425: /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */
426: MatGetSize(Amat,&M,&N);
427: MatGetOwnershipRange(Amat,&first,&last);
428: PetscMalloc1(N-last+first,¬me);
429: for (i=0; i<first; i++) notme[i] = i;
430: for (i=last; i<M; i++) notme[i-last+first] = i;
431: ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme);
432: ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me);
433: MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs);
434: Aoff = Aoffs[0];
435: MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs);
436: Boff = Boffs[0];
437: MatIsTranspose(Aoff,Boff,tol,f);
438: MatDestroyMatrices(1,&Aoffs);
439: MatDestroyMatrices(1,&Boffs);
440: ISDestroy(&Me);
441: ISDestroy(&Notme);
442: PetscFree(notme);
443: return(0);
444: }
446: PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz)
447: {
448: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
452: /* do nondiagonal part */
453: (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
454: /* do local part */
455: (*a->A->ops->multtransposeadd)(a->A,xx,yy,zz);
456: /* add partial results together */
457: VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
458: VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);
459: return(0);
460: }
462: /*
463: This only works correctly for square matrices where the subblock A->A is the
464: diagonal block
465: */
466: PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v)
467: {
469: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
472: if (A->rmap->N != A->cmap->N) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block");
473: 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");
474: MatGetDiagonal(a->A,v);
475: return(0);
476: }
478: PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa)
479: {
480: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
484: MatScale(a->A,aa);
485: MatScale(a->B,aa);
486: return(0);
487: }
489: PetscErrorCode MatDestroy_MPISELL(Mat mat)
490: {
491: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
495: #if defined(PETSC_USE_LOG)
496: PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D",mat->rmap->N,mat->cmap->N);
497: #endif
498: MatStashDestroy_Private(&mat->stash);
499: VecDestroy(&sell->diag);
500: MatDestroy(&sell->A);
501: MatDestroy(&sell->B);
502: #if defined(PETSC_USE_CTABLE)
503: PetscTableDestroy(&sell->colmap);
504: #else
505: PetscFree(sell->colmap);
506: #endif
507: PetscFree(sell->garray);
508: VecDestroy(&sell->lvec);
509: VecScatterDestroy(&sell->Mvctx);
510: PetscFree2(sell->rowvalues,sell->rowindices);
511: PetscFree(sell->ld);
512: PetscFree(mat->data);
514: PetscObjectChangeTypeName((PetscObject)mat,NULL);
515: PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);
516: PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);
517: PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL);
518: PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL);
519: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL);
520: PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);
521: return(0);
522: }
524: #include <petscdraw.h>
525: PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer)
526: {
527: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
528: PetscErrorCode ierr;
529: PetscMPIInt rank=sell->rank,size=sell->size;
530: PetscBool isdraw,iascii,isbinary;
531: PetscViewer sviewer;
532: PetscViewerFormat format;
535: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
536: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
537: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
538: if (iascii) {
539: PetscViewerGetFormat(viewer,&format);
540: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
541: MatInfo info;
542: PetscInt *inodes;
544: MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank);
545: MatGetInfo(mat,MAT_LOCAL,&info);
546: MatInodeGetInodeSizes(sell->A,NULL,&inodes,NULL);
547: PetscViewerASCIIPushSynchronized(viewer);
548: if (!inodes) {
549: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D mem %D, not using I-node routines\n",
550: rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
551: } else {
552: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D mem %D, using I-node routines\n",
553: rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);
554: }
555: MatGetInfo(sell->A,MAT_LOCAL,&info);
556: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);
557: MatGetInfo(sell->B,MAT_LOCAL,&info);
558: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);
559: PetscViewerFlush(viewer);
560: PetscViewerASCIIPopSynchronized(viewer);
561: PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n");
562: VecScatterView(sell->Mvctx,viewer);
563: return(0);
564: } else if (format == PETSC_VIEWER_ASCII_INFO) {
565: PetscInt inodecount,inodelimit,*inodes;
566: MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit);
567: if (inodes) {
568: PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %D nodes, limit used is %D\n",inodecount,inodelimit);
569: } else {
570: PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n");
571: }
572: return(0);
573: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
574: return(0);
575: }
576: } else if (isbinary) {
577: if (size == 1) {
578: PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name);
579: MatView(sell->A,viewer);
580: } else {
581: /* MatView_MPISELL_Binary(mat,viewer); */
582: }
583: return(0);
584: } else if (isdraw) {
585: PetscDraw draw;
586: PetscBool isnull;
587: PetscViewerDrawGetDraw(viewer,0,&draw);
588: PetscDrawIsNull(draw,&isnull);
589: if (isnull) return(0);
590: }
592: {
593: /* assemble the entire matrix onto first processor. */
594: Mat A;
595: Mat_SeqSELL *Aloc;
596: PetscInt M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j;
597: MatScalar *aval;
598: PetscBool isnonzero;
600: MatCreate(PetscObjectComm((PetscObject)mat),&A);
601: if (rank == 0) {
602: MatSetSizes(A,M,N,M,N);
603: } else {
604: MatSetSizes(A,0,0,M,N);
605: }
606: /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */
607: MatSetType(A,MATMPISELL);
608: MatMPISELLSetPreallocation(A,0,NULL,0,NULL);
609: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
610: PetscLogObjectParent((PetscObject)mat,(PetscObject)A);
612: /* copy over the A part */
613: Aloc = (Mat_SeqSELL*)sell->A->data;
614: acolidx = Aloc->colidx; aval = Aloc->val;
615: for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */
616: for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
617: isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
618: if (isnonzero) { /* check the mask bit */
619: row = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */
620: col = *acolidx + mat->rmap->rstart;
621: MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
622: }
623: aval++; acolidx++;
624: }
625: }
627: /* copy over the B part */
628: Aloc = (Mat_SeqSELL*)sell->B->data;
629: acolidx = Aloc->colidx; aval = Aloc->val;
630: for (i=0; i<Aloc->totalslices; i++) {
631: for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) {
632: isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]);
633: if (isnonzero) {
634: row = (i<<3)+(j&0x07) + mat->rmap->rstart;
635: col = sell->garray[*acolidx];
636: MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES);
637: }
638: aval++; acolidx++;
639: }
640: }
642: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
643: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
644: /*
645: Everyone has to call to draw the matrix since the graphics waits are
646: synchronized across all processors that share the PetscDraw object
647: */
648: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
649: if (rank == 0) {
650: PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name);
651: MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer);
652: }
653: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
654: PetscViewerFlush(viewer);
655: MatDestroy(&A);
656: }
657: return(0);
658: }
660: PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer)
661: {
663: PetscBool iascii,isdraw,issocket,isbinary;
666: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
667: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
668: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
669: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
670: if (iascii || isdraw || isbinary || issocket) {
671: MatView_MPISELL_ASCIIorDraworSocket(mat,viewer);
672: }
673: return(0);
674: }
676: PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
677: {
678: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
682: MatGetSize(sell->B,NULL,nghosts);
683: if (ghosts) *ghosts = sell->garray;
684: return(0);
685: }
687: PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info)
688: {
689: Mat_MPISELL *mat=(Mat_MPISELL*)matin->data;
690: Mat A=mat->A,B=mat->B;
692: PetscLogDouble isend[5],irecv[5];
695: info->block_size = 1.0;
696: MatGetInfo(A,MAT_LOCAL,info);
698: isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded;
699: isend[3] = info->memory; isend[4] = info->mallocs;
701: MatGetInfo(B,MAT_LOCAL,info);
703: isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded;
704: isend[3] += info->memory; isend[4] += info->mallocs;
705: if (flag == MAT_LOCAL) {
706: info->nz_used = isend[0];
707: info->nz_allocated = isend[1];
708: info->nz_unneeded = isend[2];
709: info->memory = isend[3];
710: info->mallocs = isend[4];
711: } else if (flag == MAT_GLOBAL_MAX) {
712: MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin));
714: info->nz_used = irecv[0];
715: info->nz_allocated = irecv[1];
716: info->nz_unneeded = irecv[2];
717: info->memory = irecv[3];
718: info->mallocs = irecv[4];
719: } else if (flag == MAT_GLOBAL_SUM) {
720: MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin));
722: info->nz_used = irecv[0];
723: info->nz_allocated = irecv[1];
724: info->nz_unneeded = irecv[2];
725: info->memory = irecv[3];
726: info->mallocs = irecv[4];
727: }
728: info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */
729: info->fill_ratio_needed = 0;
730: info->factor_mallocs = 0;
731: return(0);
732: }
734: PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg)
735: {
736: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
740: switch (op) {
741: case MAT_NEW_NONZERO_LOCATIONS:
742: case MAT_NEW_NONZERO_ALLOCATION_ERR:
743: case MAT_UNUSED_NONZERO_LOCATION_ERR:
744: case MAT_KEEP_NONZERO_PATTERN:
745: case MAT_NEW_NONZERO_LOCATION_ERR:
746: case MAT_USE_INODES:
747: case MAT_IGNORE_ZERO_ENTRIES:
748: MatCheckPreallocated(A,1);
749: MatSetOption(a->A,op,flg);
750: MatSetOption(a->B,op,flg);
751: break;
752: case MAT_ROW_ORIENTED:
753: MatCheckPreallocated(A,1);
754: a->roworiented = flg;
756: MatSetOption(a->A,op,flg);
757: MatSetOption(a->B,op,flg);
758: break;
759: case MAT_FORCE_DIAGONAL_ENTRIES:
760: case MAT_SORTED_FULL:
761: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
762: break;
763: case MAT_IGNORE_OFF_PROC_ENTRIES:
764: a->donotstash = flg;
765: break;
766: case MAT_SPD:
767: A->spd_set = PETSC_TRUE;
768: A->spd = flg;
769: if (flg) {
770: A->symmetric = PETSC_TRUE;
771: A->structurally_symmetric = PETSC_TRUE;
772: A->symmetric_set = PETSC_TRUE;
773: A->structurally_symmetric_set = PETSC_TRUE;
774: }
775: break;
776: case MAT_SYMMETRIC:
777: MatCheckPreallocated(A,1);
778: MatSetOption(a->A,op,flg);
779: break;
780: case MAT_STRUCTURALLY_SYMMETRIC:
781: MatCheckPreallocated(A,1);
782: MatSetOption(a->A,op,flg);
783: break;
784: case MAT_HERMITIAN:
785: MatCheckPreallocated(A,1);
786: MatSetOption(a->A,op,flg);
787: break;
788: case MAT_SYMMETRY_ETERNAL:
789: MatCheckPreallocated(A,1);
790: MatSetOption(a->A,op,flg);
791: break;
792: default:
793: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
794: }
795: return(0);
796: }
798: PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr)
799: {
800: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
801: Mat a=sell->A,b=sell->B;
803: PetscInt s1,s2,s3;
806: MatGetLocalSize(mat,&s2,&s3);
807: if (rr) {
808: VecGetLocalSize(rr,&s1);
809: if (s1!=s3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size");
810: /* Overlap communication with computation. */
811: VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
812: }
813: if (ll) {
814: VecGetLocalSize(ll,&s1);
815: if (s1!=s2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size");
816: (*b->ops->diagonalscale)(b,ll,NULL);
817: }
818: /* scale the diagonal block */
819: (*a->ops->diagonalscale)(a,ll,rr);
821: if (rr) {
822: /* Do a scatter end and then right scale the off-diagonal block */
823: VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);
824: (*b->ops->diagonalscale)(b,NULL,sell->lvec);
825: }
826: return(0);
827: }
829: PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
830: {
831: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
835: MatSetUnfactored(a->A);
836: return(0);
837: }
839: PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool *flag)
840: {
841: Mat_MPISELL *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data;
842: Mat a,b,c,d;
843: PetscBool flg;
847: a = matA->A; b = matA->B;
848: c = matB->A; d = matB->B;
850: MatEqual(a,c,&flg);
851: if (flg) {
852: MatEqual(b,d,&flg);
853: }
854: MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
855: return(0);
856: }
858: PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str)
859: {
861: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
862: Mat_MPISELL *b=(Mat_MPISELL*)B->data;
865: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
866: if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
867: /* because of the column compression in the off-processor part of the matrix a->B,
868: the number of columns in a->B and b->B may be different, hence we cannot call
869: the MatCopy() directly on the two parts. If need be, we can provide a more
870: efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices
871: then copying the submatrices */
872: MatCopy_Basic(A,B,str);
873: } else {
874: MatCopy(a->A,b->A,str);
875: MatCopy(a->B,b->B,str);
876: }
877: return(0);
878: }
880: PetscErrorCode MatSetUp_MPISELL(Mat A)
881: {
885: MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
886: return(0);
887: }
889: extern PetscErrorCode MatConjugate_SeqSELL(Mat);
891: PetscErrorCode MatConjugate_MPISELL(Mat mat)
892: {
893: #if defined(PETSC_USE_COMPLEX)
895: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
898: MatConjugate_SeqSELL(sell->A);
899: MatConjugate_SeqSELL(sell->B);
900: #else
902: #endif
903: return(0);
904: }
906: PetscErrorCode MatRealPart_MPISELL(Mat A)
907: {
908: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
912: MatRealPart(a->A);
913: MatRealPart(a->B);
914: return(0);
915: }
917: PetscErrorCode MatImaginaryPart_MPISELL(Mat A)
918: {
919: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
923: MatImaginaryPart(a->A);
924: MatImaginaryPart(a->B);
925: return(0);
926: }
928: PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values)
929: {
930: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
934: MatInvertBlockDiagonal(a->A,values);
935: A->factorerrortype = a->A->factorerrortype;
936: return(0);
937: }
939: static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx)
940: {
942: Mat_MPISELL *sell=(Mat_MPISELL*)x->data;
945: MatSetRandom(sell->A,rctx);
946: MatSetRandom(sell->B,rctx);
947: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
948: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
949: return(0);
950: }
952: PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A)
953: {
957: PetscOptionsHead(PetscOptionsObject,"MPISELL options");
958: PetscOptionsTail();
959: return(0);
960: }
962: PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a)
963: {
965: Mat_MPISELL *msell=(Mat_MPISELL*)Y->data;
966: Mat_SeqSELL *sell=(Mat_SeqSELL*)msell->A->data;
969: if (!Y->preallocated) {
970: MatMPISELLSetPreallocation(Y,1,NULL,0,NULL);
971: } else if (!sell->nz) {
972: PetscInt nonew = sell->nonew;
973: MatSeqSELLSetPreallocation(msell->A,1,NULL);
974: sell->nonew = nonew;
975: }
976: MatShift_Basic(Y,a);
977: return(0);
978: }
980: PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool *missing,PetscInt *d)
981: {
982: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
986: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices");
987: MatMissingDiagonal(a->A,missing,d);
988: if (d) {
989: PetscInt rstart;
990: MatGetOwnershipRange(A,&rstart,NULL);
991: *d += rstart;
993: }
994: return(0);
995: }
997: PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a)
998: {
1000: *a = ((Mat_MPISELL*)A->data)->A;
1001: return(0);
1002: }
1004: /* -------------------------------------------------------------------*/
1005: static struct _MatOps MatOps_Values = {MatSetValues_MPISELL,
1006: NULL,
1007: NULL,
1008: MatMult_MPISELL,
1009: /* 4*/ MatMultAdd_MPISELL,
1010: MatMultTranspose_MPISELL,
1011: MatMultTransposeAdd_MPISELL,
1012: NULL,
1013: NULL,
1014: NULL,
1015: /*10*/ NULL,
1016: NULL,
1017: NULL,
1018: MatSOR_MPISELL,
1019: NULL,
1020: /*15*/ MatGetInfo_MPISELL,
1021: MatEqual_MPISELL,
1022: MatGetDiagonal_MPISELL,
1023: MatDiagonalScale_MPISELL,
1024: NULL,
1025: /*20*/ MatAssemblyBegin_MPISELL,
1026: MatAssemblyEnd_MPISELL,
1027: MatSetOption_MPISELL,
1028: MatZeroEntries_MPISELL,
1029: /*24*/ NULL,
1030: NULL,
1031: NULL,
1032: NULL,
1033: NULL,
1034: /*29*/ MatSetUp_MPISELL,
1035: NULL,
1036: NULL,
1037: MatGetDiagonalBlock_MPISELL,
1038: NULL,
1039: /*34*/ MatDuplicate_MPISELL,
1040: NULL,
1041: NULL,
1042: NULL,
1043: NULL,
1044: /*39*/ NULL,
1045: NULL,
1046: NULL,
1047: MatGetValues_MPISELL,
1048: MatCopy_MPISELL,
1049: /*44*/ NULL,
1050: MatScale_MPISELL,
1051: MatShift_MPISELL,
1052: MatDiagonalSet_MPISELL,
1053: NULL,
1054: /*49*/ MatSetRandom_MPISELL,
1055: NULL,
1056: NULL,
1057: NULL,
1058: NULL,
1059: /*54*/ MatFDColoringCreate_MPIXAIJ,
1060: NULL,
1061: MatSetUnfactored_MPISELL,
1062: NULL,
1063: NULL,
1064: /*59*/ NULL,
1065: MatDestroy_MPISELL,
1066: MatView_MPISELL,
1067: NULL,
1068: NULL,
1069: /*64*/ NULL,
1070: NULL,
1071: NULL,
1072: NULL,
1073: NULL,
1074: /*69*/ NULL,
1075: NULL,
1076: NULL,
1077: NULL,
1078: NULL,
1079: NULL,
1080: /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */
1081: MatSetFromOptions_MPISELL,
1082: NULL,
1083: NULL,
1084: NULL,
1085: /*80*/ NULL,
1086: NULL,
1087: NULL,
1088: /*83*/ NULL,
1089: NULL,
1090: NULL,
1091: NULL,
1092: NULL,
1093: NULL,
1094: /*89*/ NULL,
1095: NULL,
1096: NULL,
1097: NULL,
1098: NULL,
1099: /*94*/ NULL,
1100: NULL,
1101: NULL,
1102: NULL,
1103: NULL,
1104: /*99*/ NULL,
1105: NULL,
1106: NULL,
1107: MatConjugate_MPISELL,
1108: NULL,
1109: /*104*/NULL,
1110: MatRealPart_MPISELL,
1111: MatImaginaryPart_MPISELL,
1112: NULL,
1113: NULL,
1114: /*109*/NULL,
1115: NULL,
1116: NULL,
1117: NULL,
1118: MatMissingDiagonal_MPISELL,
1119: /*114*/NULL,
1120: NULL,
1121: MatGetGhosts_MPISELL,
1122: NULL,
1123: NULL,
1124: /*119*/NULL,
1125: NULL,
1126: NULL,
1127: NULL,
1128: NULL,
1129: /*124*/NULL,
1130: NULL,
1131: MatInvertBlockDiagonal_MPISELL,
1132: NULL,
1133: NULL,
1134: /*129*/NULL,
1135: NULL,
1136: NULL,
1137: NULL,
1138: NULL,
1139: /*134*/NULL,
1140: NULL,
1141: NULL,
1142: NULL,
1143: NULL,
1144: /*139*/NULL,
1145: NULL,
1146: NULL,
1147: MatFDColoringSetUp_MPIXAIJ,
1148: NULL,
1149: /*144*/NULL
1150: };
1152: /* ----------------------------------------------------------------------------------------*/
1154: PetscErrorCode MatStoreValues_MPISELL(Mat mat)
1155: {
1156: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
1160: MatStoreValues(sell->A);
1161: MatStoreValues(sell->B);
1162: return(0);
1163: }
1165: PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
1166: {
1167: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
1171: MatRetrieveValues(sell->A);
1172: MatRetrieveValues(sell->B);
1173: return(0);
1174: }
1176: PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[])
1177: {
1178: Mat_MPISELL *b;
1182: PetscLayoutSetUp(B->rmap);
1183: PetscLayoutSetUp(B->cmap);
1184: b = (Mat_MPISELL*)B->data;
1186: if (!B->preallocated) {
1187: /* Explicitly create 2 MATSEQSELL matrices. */
1188: MatCreate(PETSC_COMM_SELF,&b->A);
1189: MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
1190: MatSetBlockSizesFromMats(b->A,B,B);
1191: MatSetType(b->A,MATSEQSELL);
1192: PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
1193: MatCreate(PETSC_COMM_SELF,&b->B);
1194: MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
1195: MatSetBlockSizesFromMats(b->B,B,B);
1196: MatSetType(b->B,MATSEQSELL);
1197: PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
1198: }
1200: MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);
1201: MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);
1202: B->preallocated = PETSC_TRUE;
1203: B->was_assembled = PETSC_FALSE;
1204: /*
1205: critical for MatAssemblyEnd to work.
1206: MatAssemblyBegin checks it to set up was_assembled
1207: and MatAssemblyEnd checks was_assembled to determine whether to build garray
1208: */
1209: B->assembled = PETSC_FALSE;
1210: return(0);
1211: }
1213: PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat)
1214: {
1215: Mat mat;
1216: Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data;
1220: *newmat = NULL;
1221: MatCreate(PetscObjectComm((PetscObject)matin),&mat);
1222: MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);
1223: MatSetBlockSizesFromMats(mat,matin,matin);
1224: MatSetType(mat,((PetscObject)matin)->type_name);
1225: a = (Mat_MPISELL*)mat->data;
1227: mat->factortype = matin->factortype;
1228: mat->assembled = PETSC_TRUE;
1229: mat->insertmode = NOT_SET_VALUES;
1230: mat->preallocated = PETSC_TRUE;
1232: a->size = oldmat->size;
1233: a->rank = oldmat->rank;
1234: a->donotstash = oldmat->donotstash;
1235: a->roworiented = oldmat->roworiented;
1236: a->rowindices = NULL;
1237: a->rowvalues = NULL;
1238: a->getrowactive = PETSC_FALSE;
1240: PetscLayoutReference(matin->rmap,&mat->rmap);
1241: PetscLayoutReference(matin->cmap,&mat->cmap);
1243: if (oldmat->colmap) {
1244: #if defined(PETSC_USE_CTABLE)
1245: PetscTableCreateCopy(oldmat->colmap,&a->colmap);
1246: #else
1247: PetscMalloc1(mat->cmap->N,&a->colmap);
1248: PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));
1249: PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N);
1250: #endif
1251: } else a->colmap = NULL;
1252: if (oldmat->garray) {
1253: PetscInt len;
1254: len = oldmat->B->cmap->n;
1255: PetscMalloc1(len+1,&a->garray);
1256: PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));
1257: if (len) { PetscArraycpy(a->garray,oldmat->garray,len); }
1258: } else a->garray = NULL;
1260: VecDuplicate(oldmat->lvec,&a->lvec);
1261: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);
1262: VecScatterCopy(oldmat->Mvctx,&a->Mvctx);
1263: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);
1264: MatDuplicate(oldmat->A,cpvalues,&a->A);
1265: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);
1266: MatDuplicate(oldmat->B,cpvalues,&a->B);
1267: PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);
1268: PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);
1269: *newmat = mat;
1270: return(0);
1271: }
1273: /*@C
1274: MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format.
1275: For good matrix assembly performance the user should preallocate the matrix storage by
1276: setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz).
1278: Collective
1280: Input Parameters:
1281: + B - the matrix
1282: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1283: (same value is used for all local rows)
1284: . d_nnz - array containing the number of nonzeros in the various rows of the
1285: DIAGONAL portion of the local submatrix (possibly different for each row)
1286: or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1287: The size of this array is equal to the number of local rows, i.e 'm'.
1288: For matrices that will be factored, you must leave room for (and set)
1289: the diagonal entry even if it is zero.
1290: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1291: submatrix (same value is used for all local rows).
1292: - o_nnz - array containing the number of nonzeros in the various rows of the
1293: OFF-DIAGONAL portion of the local submatrix (possibly different for
1294: each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1295: structure. The size of this array is equal to the number
1296: of local rows, i.e 'm'.
1298: If the *_nnz parameter is given then the *_nz parameter is ignored
1300: The stored row and column indices begin with zero.
1302: The parallel matrix is partitioned such that the first m0 rows belong to
1303: process 0, the next m1 rows belong to process 1, the next m2 rows belong
1304: to process 2 etc.. where m0,m1,m2... are the input parameter 'm'.
1306: The DIAGONAL portion of the local submatrix of a processor can be defined
1307: as the submatrix which is obtained by extraction the part corresponding to
1308: the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
1309: first row that belongs to the processor, r2 is the last row belonging to
1310: the this processor, and c1-c2 is range of indices of the local part of a
1311: vector suitable for applying the matrix to. This is an mxn matrix. In the
1312: common case of a square matrix, the row and column ranges are the same and
1313: the DIAGONAL part is also square. The remaining portion of the local
1314: submatrix (mxN) constitute the OFF-DIAGONAL portion.
1316: If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored.
1318: You can call MatGetInfo() to get information on how effective the preallocation was;
1319: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
1320: You can also run with the option -info and look for messages with the string
1321: malloc in them to see if additional memory allocation was needed.
1323: Example usage:
1325: Consider the following 8x8 matrix with 34 non-zero values, that is
1326: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1327: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1328: as follows:
1330: .vb
1331: 1 2 0 | 0 3 0 | 0 4
1332: Proc0 0 5 6 | 7 0 0 | 8 0
1333: 9 0 10 | 11 0 0 | 12 0
1334: -------------------------------------
1335: 13 0 14 | 15 16 17 | 0 0
1336: Proc1 0 18 0 | 19 20 21 | 0 0
1337: 0 0 0 | 22 23 0 | 24 0
1338: -------------------------------------
1339: Proc2 25 26 27 | 0 0 28 | 29 0
1340: 30 0 0 | 31 32 33 | 0 34
1341: .ve
1343: This can be represented as a collection of submatrices as:
1345: .vb
1346: A B C
1347: D E F
1348: G H I
1349: .ve
1351: Where the submatrices A,B,C are owned by proc0, D,E,F are
1352: owned by proc1, G,H,I are owned by proc2.
1354: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1355: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1356: The 'M','N' parameters are 8,8, and have the same values on all procs.
1358: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1359: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1360: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1361: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1362: part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1363: matrix, ans [DF] as another SeqSELL matrix.
1365: When d_nz, o_nz parameters are specified, d_nz storage elements are
1366: allocated for every row of the local diagonal submatrix, and o_nz
1367: storage locations are allocated for every row of the OFF-DIAGONAL submat.
1368: One way to choose d_nz and o_nz is to use the max nonzerors per local
1369: rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1370: In this case, the values of d_nz,o_nz are:
1371: .vb
1372: proc0 : dnz = 2, o_nz = 2
1373: proc1 : dnz = 3, o_nz = 2
1374: proc2 : dnz = 1, o_nz = 4
1375: .ve
1376: We are allocating m*(d_nz+o_nz) storage locations for every proc. This
1377: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1378: for proc3. i.e we are using 12+15+10=37 storage locations to store
1379: 34 values.
1381: When d_nnz, o_nnz parameters are specified, the storage is specified
1382: for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1383: In the above case the values for d_nnz,o_nnz are:
1384: .vb
1385: proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
1386: proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
1387: proc2: d_nnz = [1,1] and o_nnz = [4,4]
1388: .ve
1389: Here the space allocated is according to nz (or maximum values in the nnz
1390: if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37
1392: Level: intermediate
1394: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(),
1395: MATMPISELL, MatGetInfo(), PetscSplitOwnership()
1396: @*/
1397: PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1398: {
1404: PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1405: return(0);
1406: }
1408: /*@C
1409: MatCreateSELL - Creates a sparse parallel matrix in SELL format.
1411: Collective
1413: Input Parameters:
1414: + comm - MPI communicator
1415: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
1416: This value should be the same as the local size used in creating the
1417: y vector for the matrix-vector product y = Ax.
1418: . n - This value should be the same as the local size used in creating the
1419: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
1420: calculated if N is given) For square matrices n is almost always m.
1421: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
1422: . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
1423: . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix
1424: (same value is used for all local rows)
1425: . d_rlen - array containing the number of nonzeros in the various rows of the
1426: DIAGONAL portion of the local submatrix (possibly different for each row)
1427: or NULL, if d_rlenmax is used to specify the nonzero structure.
1428: The size of this array is equal to the number of local rows, i.e 'm'.
1429: . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local
1430: submatrix (same value is used for all local rows).
1431: - o_rlen - array containing the number of nonzeros in the various rows of the
1432: OFF-DIAGONAL portion of the local submatrix (possibly different for
1433: each row) or NULL, if o_rlenmax is used to specify the nonzero
1434: structure. The size of this array is equal to the number
1435: of local rows, i.e 'm'.
1437: Output Parameter:
1438: . A - the matrix
1440: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
1441: MatXXXXSetPreallocation() paradigm instead of this routine directly.
1442: [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation]
1444: Notes:
1445: If the *_rlen parameter is given then the *_rlenmax parameter is ignored
1447: m,n,M,N parameters specify the size of the matrix, and its partitioning across
1448: processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate
1449: storage requirements for this matrix.
1451: If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one
1452: processor than it must be used on all processors that share the object for
1453: that argument.
1455: The user MUST specify either the local or global matrix dimensions
1456: (possibly both).
1458: The parallel matrix is partitioned across processors such that the
1459: first m0 rows belong to process 0, the next m1 rows belong to
1460: process 1, the next m2 rows belong to process 2 etc.. where
1461: m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores
1462: values corresponding to [m x N] submatrix.
1464: The columns are logically partitioned with the n0 columns belonging
1465: to 0th partition, the next n1 columns belonging to the next
1466: partition etc.. where n0,n1,n2... are the input parameter 'n'.
1468: The DIAGONAL portion of the local submatrix on any given processor
1469: is the submatrix corresponding to the rows and columns m,n
1470: corresponding to the given processor. i.e diagonal matrix on
1471: process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
1472: etc. The remaining portion of the local submatrix [m x (N-n)]
1473: constitute the OFF-DIAGONAL portion. The example below better
1474: illustrates this concept.
1476: For a square global matrix we define each processor's diagonal portion
1477: to be its local rows and the corresponding columns (a square submatrix);
1478: each processor's off-diagonal portion encompasses the remainder of the
1479: local matrix (a rectangular submatrix).
1481: If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored.
1483: When calling this routine with a single process communicator, a matrix of
1484: type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this
1485: type of communicator, use the construction mechanism:
1486: MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...);
1488: Options Database Keys:
1489: - -mat_sell_oneindex - Internally use indexing starting at 1
1490: rather than 0. Note that when calling MatSetValues(),
1491: the user still MUST index entries starting at 0!
1493: Example usage:
1495: Consider the following 8x8 matrix with 34 non-zero values, that is
1496: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1497: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1498: as follows:
1500: .vb
1501: 1 2 0 | 0 3 0 | 0 4
1502: Proc0 0 5 6 | 7 0 0 | 8 0
1503: 9 0 10 | 11 0 0 | 12 0
1504: -------------------------------------
1505: 13 0 14 | 15 16 17 | 0 0
1506: Proc1 0 18 0 | 19 20 21 | 0 0
1507: 0 0 0 | 22 23 0 | 24 0
1508: -------------------------------------
1509: Proc2 25 26 27 | 0 0 28 | 29 0
1510: 30 0 0 | 31 32 33 | 0 34
1511: .ve
1513: This can be represented as a collection of submatrices as:
1515: .vb
1516: A B C
1517: D E F
1518: G H I
1519: .ve
1521: Where the submatrices A,B,C are owned by proc0, D,E,F are
1522: owned by proc1, G,H,I are owned by proc2.
1524: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1525: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1526: The 'M','N' parameters are 8,8, and have the same values on all procs.
1528: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1529: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1530: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1531: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1532: part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL
1533: matrix, ans [DF] as another SeqSELL matrix.
1535: When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are
1536: allocated for every row of the local diagonal submatrix, and o_rlenmax
1537: storage locations are allocated for every row of the OFF-DIAGONAL submat.
1538: One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local
1539: rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1540: In this case, the values of d_rlenmax,o_rlenmax are:
1541: .vb
1542: proc0 : d_rlenmax = 2, o_rlenmax = 2
1543: proc1 : d_rlenmax = 3, o_rlenmax = 2
1544: proc2 : d_rlenmax = 1, o_rlenmax = 4
1545: .ve
1546: We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This
1547: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1548: for proc3. i.e we are using 12+15+10=37 storage locations to store
1549: 34 values.
1551: When d_rlen, o_rlen parameters are specified, the storage is specified
1552: for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1553: In the above case the values for d_nnz,o_nnz are:
1554: .vb
1555: proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2]
1556: proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1]
1557: proc2: d_nnz = [1,1] and o_nnz = [4,4]
1558: .ve
1559: Here the space allocated is still 37 though there are 34 nonzeros because
1560: the allocation is always done according to rlenmax.
1562: Level: intermediate
1564: .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(),
1565: MATMPISELL, MatCreateMPISELLWithArrays()
1566: @*/
1567: 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)
1568: {
1570: PetscMPIInt size;
1573: MatCreate(comm,A);
1574: MatSetSizes(*A,m,n,M,N);
1575: MPI_Comm_size(comm,&size);
1576: if (size > 1) {
1577: MatSetType(*A,MATMPISELL);
1578: MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);
1579: } else {
1580: MatSetType(*A,MATSEQSELL);
1581: MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);
1582: }
1583: return(0);
1584: }
1586: PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[])
1587: {
1588: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1589: PetscBool flg;
1593: PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);
1594: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input");
1595: if (Ad) *Ad = a->A;
1596: if (Ao) *Ao = a->B;
1597: if (colmap) *colmap = a->garray;
1598: return(0);
1599: }
1601: /*@C
1602: MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns
1604: Not Collective
1606: Input Parameters:
1607: + A - the matrix
1608: . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
1609: - row, col - index sets of rows and columns to extract (or NULL)
1611: Output Parameter:
1612: . A_loc - the local sequential matrix generated
1614: Level: developer
1616: .seealso: MatGetOwnershipRange(), MatMPISELLGetLocalMat()
1618: @*/
1619: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc)
1620: {
1621: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1623: PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx;
1624: IS isrowa,iscola;
1625: Mat *aloc;
1626: PetscBool match;
1629: PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match);
1630: if (!match) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input");
1631: PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0);
1632: if (!row) {
1633: start = A->rmap->rstart; end = A->rmap->rend;
1634: ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa);
1635: } else {
1636: isrowa = *row;
1637: }
1638: if (!col) {
1639: start = A->cmap->rstart;
1640: cmap = a->garray;
1641: nzA = a->A->cmap->n;
1642: nzB = a->B->cmap->n;
1643: PetscMalloc1(nzA+nzB, &idx);
1644: ncols = 0;
1645: for (i=0; i<nzB; i++) {
1646: if (cmap[i] < start) idx[ncols++] = cmap[i];
1647: else break;
1648: }
1649: imark = i;
1650: for (i=0; i<nzA; i++) idx[ncols++] = start + i;
1651: for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i];
1652: ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola);
1653: } else {
1654: iscola = *col;
1655: }
1656: if (scall != MAT_INITIAL_MATRIX) {
1657: PetscMalloc1(1,&aloc);
1658: aloc[0] = *A_loc;
1659: }
1660: MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc);
1661: *A_loc = aloc[0];
1662: PetscFree(aloc);
1663: if (!row) {
1664: ISDestroy(&isrowa);
1665: }
1666: if (!col) {
1667: ISDestroy(&iscola);
1668: }
1669: PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0);
1670: return(0);
1671: }
1673: #include <../src/mat/impls/aij/mpi/mpiaij.h>
1675: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1676: {
1678: Mat_MPISELL *a=(Mat_MPISELL*)A->data;
1679: Mat B;
1680: Mat_MPIAIJ *b;
1683: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled");
1685: if (reuse == MAT_REUSE_MATRIX) {
1686: B = *newmat;
1687: } else {
1688: MatCreate(PetscObjectComm((PetscObject)A),&B);
1689: MatSetType(B,MATMPIAIJ);
1690: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1691: MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1692: MatSeqAIJSetPreallocation(B,0,NULL);
1693: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1694: }
1695: b = (Mat_MPIAIJ*) B->data;
1697: if (reuse == MAT_REUSE_MATRIX) {
1698: MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);
1699: MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);
1700: } else {
1701: MatDestroy(&b->A);
1702: MatDestroy(&b->B);
1703: MatDisAssemble_MPISELL(A);
1704: MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);
1705: MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);
1706: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1707: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1708: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1709: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1710: }
1712: if (reuse == MAT_INPLACE_MATRIX) {
1713: MatHeaderReplace(A,&B);
1714: } else {
1715: *newmat = B;
1716: }
1717: return(0);
1718: }
1720: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1721: {
1723: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data;
1724: Mat B;
1725: Mat_MPISELL *b;
1728: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled");
1730: if (reuse == MAT_REUSE_MATRIX) {
1731: B = *newmat;
1732: } else {
1733: MatCreate(PetscObjectComm((PetscObject)A),&B);
1734: MatSetType(B,MATMPISELL);
1735: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1736: MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
1737: MatSeqAIJSetPreallocation(B,0,NULL);
1738: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1739: }
1740: b = (Mat_MPISELL*) B->data;
1742: if (reuse == MAT_REUSE_MATRIX) {
1743: MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);
1744: MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);
1745: } else {
1746: MatDestroy(&b->A);
1747: MatDestroy(&b->B);
1748: MatDisAssemble_MPIAIJ(A);
1749: MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);
1750: MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);
1751: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1752: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1753: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1754: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1755: }
1757: if (reuse == MAT_INPLACE_MATRIX) {
1758: MatHeaderReplace(A,&B);
1759: } else {
1760: *newmat = B;
1761: }
1762: return(0);
1763: }
1765: PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1766: {
1767: Mat_MPISELL *mat=(Mat_MPISELL*)matin->data;
1769: Vec bb1=NULL;
1772: if (flag == SOR_APPLY_UPPER) {
1773: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1774: return(0);
1775: }
1777: if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) {
1778: VecDuplicate(bb,&bb1);
1779: }
1781: if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1782: if (flag & SOR_ZERO_INITIAL_GUESS) {
1783: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1784: its--;
1785: }
1787: while (its--) {
1788: VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1789: VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1791: /* update rhs: bb1 = bb - B*x */
1792: VecScale(mat->lvec,-1.0);
1793: (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);
1795: /* local sweep */
1796: (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);
1797: }
1798: } else if (flag & SOR_LOCAL_FORWARD_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: }
1803: while (its--) {
1804: VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1805: VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);
1807: /* update rhs: bb1 = bb - B*x */
1808: VecScale(mat->lvec,-1.0);
1809: (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);
1811: /* local sweep */
1812: (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);
1813: }
1814: } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1815: if (flag & SOR_ZERO_INITIAL_GUESS) {
1816: (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);
1817: its--;
1818: }
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_BACKWARD_SWEEP,fshift,lits,1,xx);
1829: }
1830: } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported");
1832: VecDestroy(&bb1);
1834: matin->factorerrortype = mat->A->factorerrortype;
1835: return(0);
1836: }
1838: /*MC
1839: MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices.
1841: Options Database Keys:
1842: . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions()
1844: Level: beginner
1846: .seealso: MatCreateSELL()
1847: M*/
1848: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1849: {
1850: Mat_MPISELL *b;
1852: PetscMPIInt size;
1855: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1856: PetscNewLog(B,&b);
1857: B->data = (void*)b;
1858: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1859: B->assembled = PETSC_FALSE;
1860: B->insertmode = NOT_SET_VALUES;
1861: b->size = size;
1862: MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);
1863: /* build cache for off array entries formed */
1864: MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);
1866: b->donotstash = PETSC_FALSE;
1867: b->colmap = NULL;
1868: b->garray = NULL;
1869: b->roworiented = PETSC_TRUE;
1871: /* stuff used for matrix vector multiply */
1872: b->lvec = NULL;
1873: b->Mvctx = NULL;
1875: /* stuff for MatGetRow() */
1876: b->rowindices = NULL;
1877: b->rowvalues = NULL;
1878: b->getrowactive = PETSC_FALSE;
1880: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);
1881: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);
1882: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);
1883: PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);
1884: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);
1885: PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);
1886: PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);
1887: return(0);
1888: }