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