Actual source code: matis.c
petsc-3.8.4 2018-03-24
1: /*
2: Creates a matrix class for using the Neumann-Neumann type preconditioners.
3: This stores the matrices in globally unassembled form. Each processor
4: assembles only its local Neumann problem and the parallel matrix vector
5: product is handled "implicitly".
7: Currently this allows for only one subdomain per processor.
8: */
10: #include <../src/mat/impls/is/matis.h>
11: #include <../src/mat/impls/aij/mpi/mpiaij.h>
12: #include <petsc/private/sfimpl.h>
14: #define MATIS_MAX_ENTRIES_INSERTION 2048
16: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
17: {
18: MatISLocalFields lf = (MatISLocalFields)ptr;
19: PetscInt i;
20: PetscErrorCode ierr;
23: for (i=0;i<lf->nr;i++) {
24: ISDestroy(&lf->rf[i]);
25: }
26: for (i=0;i<lf->nc;i++) {
27: ISDestroy(&lf->cf[i]);
28: }
29: PetscFree2(lf->rf,lf->cf);
30: PetscFree(lf);
31: return(0);
32: }
34: static PetscErrorCode MatISContainerDestroyArray_Private(void *ptr)
35: {
39: PetscFree(ptr);
40: return(0);
41: }
43: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
44: {
45: Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data;
46: Mat_SeqAIJ *diag = (Mat_SeqAIJ*)(aij->A->data);
47: Mat_SeqAIJ *offd = (Mat_SeqAIJ*)(aij->B->data);
48: Mat lA;
49: ISLocalToGlobalMapping rl2g,cl2g;
50: IS is;
51: MPI_Comm comm;
52: void *ptrs[2];
53: const char *names[2] = {"_convert_csr_aux","_convert_csr_data"};
54: PetscScalar *dd,*od,*aa,*data;
55: PetscInt *di,*dj,*oi,*oj;
56: PetscInt *aux,*ii,*jj;
57: PetscInt lc,dr,dc,oc,str,stc,nnz,i,jd,jo,cum;
58: PetscErrorCode ierr;
61: PetscObjectGetComm((PetscObject)A,&comm);
62: if (!aij->garray) SETERRQ(comm,PETSC_ERR_SUP,"garray not present");
64: /* access relevant information from MPIAIJ */
65: MatGetOwnershipRange(A,&str,NULL);
66: MatGetOwnershipRangeColumn(A,&stc,NULL);
67: MatGetLocalSize(A,&dr,&dc);
68: di = diag->i;
69: dj = diag->j;
70: dd = diag->a;
71: oc = aij->B->cmap->n;
72: oi = offd->i;
73: oj = offd->j;
74: od = offd->a;
75: nnz = diag->i[dr] + offd->i[dr];
77: /* generate l2g maps for rows and cols */
78: ISCreateStride(comm,dr,str,1,&is);
79: ISLocalToGlobalMappingCreateIS(is,&rl2g);
80: ISDestroy(&is);
81: if (dr) {
82: PetscMalloc1(dc+oc,&aux);
83: for (i=0; i<dc; i++) aux[i] = i+stc;
84: for (i=0; i<oc; i++) aux[i+dc] = aij->garray[i];
85: ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
86: lc = dc+oc;
87: } else {
88: ISCreateGeneral(comm,0,NULL,PETSC_OWN_POINTER,&is);
89: lc = 0;
90: }
91: ISLocalToGlobalMappingCreateIS(is,&cl2g);
92: ISDestroy(&is);
94: /* create MATIS object */
95: MatCreate(comm,newmat);
96: MatSetSizes(*newmat,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
97: MatSetType(*newmat,MATIS);
98: MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
99: ISLocalToGlobalMappingDestroy(&rl2g);
100: ISLocalToGlobalMappingDestroy(&cl2g);
102: /* merge local matrices */
103: PetscMalloc1(nnz+dr+1,&aux);
104: PetscMalloc1(nnz,&data);
105: ii = aux;
106: jj = aux+dr+1;
107: aa = data;
108: *ii = *(di++) + *(oi++);
109: for (jd=0,jo=0,cum=0;*ii<nnz;cum++)
110: {
111: for (;jd<*di;jd++) { *jj++ = *dj++; *aa++ = *dd++; }
112: for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; }
113: *(++ii) = *(di++) + *(oi++);
114: }
115: for (;cum<dr;cum++) *(++ii) = nnz;
116: ii = aux;
117: jj = aux+dr+1;
118: aa = data;
119: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);
121: /* create containers to destroy the data */
122: ptrs[0] = aux;
123: ptrs[1] = data;
124: for (i=0; i<2; i++) {
125: PetscContainer c;
127: PetscContainerCreate(PETSC_COMM_SELF,&c);
128: PetscContainerSetPointer(c,ptrs[i]);
129: PetscContainerSetUserDestroy(c,MatISContainerDestroyArray_Private);
130: PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
131: PetscContainerDestroy(&c);
132: }
134: /* finalize matrix */
135: MatISSetLocalMat(*newmat,lA);
136: MatDestroy(&lA);
137: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
139: return(0);
140: }
142: /*@
143: MatISSetUpSF - Setup star forest objects used by MatIS.
145: Collective on MPI_Comm
147: Input Parameters:
148: + A - the matrix
150: Level: advanced
152: Notes: This function does not need to be called by the user.
154: .keywords: matrix
156: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatISGetLocalMat()
157: @*/
158: PetscErrorCode MatISSetUpSF(Mat A)
159: {
165: PetscTryMethod(A,"MatISSetUpSF_C",(Mat),(A));
166: return(0);
167: }
169: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
170: {
171: Mat **nest,*snest,**rnest,lA,B;
172: IS *iscol,*isrow,*islrow,*islcol;
173: ISLocalToGlobalMapping rl2g,cl2g;
174: MPI_Comm comm;
175: PetscInt *lr,*lc,*l2gidxs;
176: PetscInt i,j,nr,nc,rbs,cbs;
177: PetscBool convert,lreuse,*istrans;
178: PetscErrorCode ierr;
181: MatNestGetSubMats(A,&nr,&nc,&nest);
182: lreuse = PETSC_FALSE;
183: rnest = NULL;
184: if (reuse == MAT_REUSE_MATRIX) {
185: PetscBool ismatis,isnest;
187: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
188: if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
189: MatISGetLocalMat(*newmat,&lA);
190: PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
191: if (isnest) {
192: MatNestGetSubMats(lA,&i,&j,&rnest);
193: lreuse = (PetscBool)(i == nr && j == nc);
194: if (!lreuse) rnest = NULL;
195: }
196: }
197: PetscObjectGetComm((PetscObject)A,&comm);
198: PetscCalloc2(nr,&lr,nc,&lc);
199: PetscCalloc6(nr,&isrow,nc,&iscol,
200: nr,&islrow,nc,&islcol,
201: nr*nc,&snest,nr*nc,&istrans);
202: MatNestGetISs(A,isrow,iscol);
203: for (i=0;i<nr;i++) {
204: for (j=0;j<nc;j++) {
205: PetscBool ismatis;
206: PetscInt l1,l2,lb1,lb2,ij=i*nc+j;
208: /* Null matrix pointers are allowed in MATNEST */
209: if (!nest[i][j]) continue;
211: /* Nested matrices should be of type MATIS */
212: PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
213: if (istrans[ij]) {
214: Mat T,lT;
215: MatTransposeGetMat(nest[i][j],&T);
216: PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
217: if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) (transposed) is not of type MATIS",i,j);
218: MatISGetLocalMat(T,&lT);
219: MatCreateTranspose(lT,&snest[ij]);
220: } else {
221: PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
222: if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
223: MatISGetLocalMat(nest[i][j],&snest[ij]);
224: }
226: /* Check compatibility of local sizes */
227: MatGetSize(snest[ij],&l1,&l2);
228: MatGetBlockSizes(snest[ij],&lb1,&lb2);
229: if (!l1 || !l2) continue;
230: if (lr[i] && l1 != lr[i]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lr[i],l1);
231: if (lc[j] && l2 != lc[j]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lc[j],l2);
232: lr[i] = l1;
233: lc[j] = l2;
235: /* check compatibilty for local matrix reusage */
236: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
237: }
238: }
240: #if defined (PETSC_USE_DEBUG)
241: /* Check compatibility of l2g maps for rows */
242: for (i=0;i<nr;i++) {
243: rl2g = NULL;
244: for (j=0;j<nc;j++) {
245: PetscInt n1,n2;
247: if (!nest[i][j]) continue;
248: if (istrans[i*nc+j]) {
249: Mat T;
251: MatTransposeGetMat(nest[i][j],&T);
252: MatGetLocalToGlobalMapping(T,NULL,&cl2g);
253: } else {
254: MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
255: }
256: ISLocalToGlobalMappingGetSize(cl2g,&n1);
257: if (!n1) continue;
258: if (!rl2g) {
259: rl2g = cl2g;
260: } else {
261: const PetscInt *idxs1,*idxs2;
262: PetscBool same;
264: ISLocalToGlobalMappingGetSize(rl2g,&n2);
265: if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap size %D != %D",i,j,n1,n2);
266: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
267: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
268: PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
269: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
270: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
271: if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap",i,j);
272: }
273: }
274: }
275: /* Check compatibility of l2g maps for columns */
276: for (i=0;i<nc;i++) {
277: rl2g = NULL;
278: for (j=0;j<nr;j++) {
279: PetscInt n1,n2;
281: if (!nest[j][i]) continue;
282: if (istrans[j*nc+i]) {
283: Mat T;
285: MatTransposeGetMat(nest[j][i],&T);
286: MatGetLocalToGlobalMapping(T,&cl2g,NULL);
287: } else {
288: MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
289: }
290: ISLocalToGlobalMappingGetSize(cl2g,&n1);
291: if (!n1) continue;
292: if (!rl2g) {
293: rl2g = cl2g;
294: } else {
295: const PetscInt *idxs1,*idxs2;
296: PetscBool same;
298: ISLocalToGlobalMappingGetSize(rl2g,&n2);
299: if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap size %D != %D",j,i,n1,n2);
300: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
301: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
302: PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
303: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
304: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
305: if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap",j,i);
306: }
307: }
308: }
309: #endif
311: B = NULL;
312: if (reuse != MAT_REUSE_MATRIX) {
313: PetscInt stl;
315: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
316: for (i=0,stl=0;i<nr;i++) stl += lr[i];
317: PetscMalloc1(stl,&l2gidxs);
318: for (i=0,stl=0;i<nr;i++) {
319: Mat usedmat;
320: Mat_IS *matis;
321: const PetscInt *idxs;
323: /* local IS for local NEST */
324: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
326: /* l2gmap */
327: j = 0;
328: usedmat = nest[i][j];
329: while (!usedmat && j < nc-1) usedmat = nest[i][++j];
330: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");
332: if (istrans[i*nc+j]) {
333: Mat T;
334: MatTransposeGetMat(usedmat,&T);
335: usedmat = T;
336: }
337: MatISSetUpSF(usedmat);
338: matis = (Mat_IS*)(usedmat->data);
339: ISGetIndices(isrow[i],&idxs);
340: if (istrans[i*nc+j]) {
341: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
342: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
343: } else {
344: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
345: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
346: }
347: ISRestoreIndices(isrow[i],&idxs);
348: stl += lr[i];
349: }
350: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);
352: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
353: for (i=0,stl=0;i<nc;i++) stl += lc[i];
354: PetscMalloc1(stl,&l2gidxs);
355: for (i=0,stl=0;i<nc;i++) {
356: Mat usedmat;
357: Mat_IS *matis;
358: const PetscInt *idxs;
360: /* local IS for local NEST */
361: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
363: /* l2gmap */
364: j = 0;
365: usedmat = nest[j][i];
366: while (!usedmat && j < nr-1) usedmat = nest[++j][i];
367: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat");
368: if (istrans[j*nc+i]) {
369: Mat T;
370: MatTransposeGetMat(usedmat,&T);
371: usedmat = T;
372: }
373: MatISSetUpSF(usedmat);
374: matis = (Mat_IS*)(usedmat->data);
375: ISGetIndices(iscol[i],&idxs);
376: if (istrans[j*nc+i]) {
377: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
378: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
379: } else {
380: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
381: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
382: }
383: ISRestoreIndices(iscol[i],&idxs);
384: stl += lc[i];
385: }
386: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);
388: /* Create MATIS */
389: MatCreate(comm,&B);
390: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
391: MatGetBlockSizes(A,&rbs,&cbs);
392: MatSetBlockSizes(B,rbs,cbs);
393: MatSetType(B,MATIS);
394: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
395: ISLocalToGlobalMappingDestroy(&rl2g);
396: ISLocalToGlobalMappingDestroy(&cl2g);
397: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
398: for (i=0;i<nr*nc;i++) {
399: if (istrans[i]) {
400: MatDestroy(&snest[i]);
401: }
402: }
403: MatISSetLocalMat(B,lA);
404: MatDestroy(&lA);
405: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
406: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
407: if (reuse == MAT_INPLACE_MATRIX) {
408: MatHeaderReplace(A,&B);
409: } else {
410: *newmat = B;
411: }
412: } else {
413: if (lreuse) {
414: MatISGetLocalMat(*newmat,&lA);
415: for (i=0;i<nr;i++) {
416: for (j=0;j<nc;j++) {
417: if (snest[i*nc+j]) {
418: MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
419: if (istrans[i*nc+j]) {
420: MatDestroy(&snest[i*nc+j]);
421: }
422: }
423: }
424: }
425: } else {
426: PetscInt stl;
427: for (i=0,stl=0;i<nr;i++) {
428: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
429: stl += lr[i];
430: }
431: for (i=0,stl=0;i<nc;i++) {
432: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
433: stl += lc[i];
434: }
435: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
436: for (i=0;i<nr*nc;i++) {
437: if (istrans[i]) {
438: MatDestroy(&snest[i]);
439: }
440: }
441: MatISSetLocalMat(*newmat,lA);
442: MatDestroy(&lA);
443: }
444: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
445: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
446: }
448: /* Create local matrix in MATNEST format */
449: convert = PETSC_FALSE;
450: PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
451: if (convert) {
452: Mat M;
453: MatISLocalFields lf;
454: PetscContainer c;
456: MatISGetLocalMat(*newmat,&lA);
457: MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
458: MatISSetLocalMat(*newmat,M);
459: MatDestroy(&M);
461: /* attach local fields to the matrix */
462: PetscNew(&lf);
463: PetscCalloc2(nr,&lf->rf,nc,&lf->cf);
464: for (i=0;i<nr;i++) {
465: PetscInt n,st;
467: ISGetLocalSize(islrow[i],&n);
468: ISStrideGetInfo(islrow[i],&st,NULL);
469: ISCreateStride(comm,n,st,1,&lf->rf[i]);
470: }
471: for (i=0;i<nc;i++) {
472: PetscInt n,st;
474: ISGetLocalSize(islcol[i],&n);
475: ISStrideGetInfo(islcol[i],&st,NULL);
476: ISCreateStride(comm,n,st,1,&lf->cf[i]);
477: }
478: lf->nr = nr;
479: lf->nc = nc;
480: PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
481: PetscContainerSetPointer(c,lf);
482: PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
483: PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
484: PetscContainerDestroy(&c);
485: }
487: /* Free workspace */
488: for (i=0;i<nr;i++) {
489: ISDestroy(&islrow[i]);
490: }
491: for (i=0;i<nc;i++) {
492: ISDestroy(&islcol[i]);
493: }
494: PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
495: PetscFree2(lr,lc);
496: return(0);
497: }
499: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
500: {
501: Mat_IS *matis = (Mat_IS*)A->data;
502: Vec ll,rr;
503: const PetscScalar *Y,*X;
504: PetscScalar *x,*y;
505: PetscErrorCode ierr;
508: MatISSetUpSF(A);
509: if (l) {
510: ll = matis->y;
511: VecGetArrayRead(l,&Y);
512: VecGetArray(ll,&y);
513: PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y);
514: } else {
515: ll = NULL;
516: }
517: if (r) {
518: rr = matis->x;
519: VecGetArrayRead(r,&X);
520: VecGetArray(rr,&x);
521: PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x);
522: } else {
523: rr = NULL;
524: }
525: if (ll) {
526: PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y);
527: VecRestoreArrayRead(l,&Y);
528: VecRestoreArray(ll,&y);
529: }
530: if (rr) {
531: PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x);
532: VecRestoreArrayRead(r,&X);
533: VecRestoreArray(rr,&x);
534: }
535: MatDiagonalScale(matis->A,ll,rr);
536: return(0);
537: }
539: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
540: {
541: Mat_IS *matis = (Mat_IS*)A->data;
542: MatInfo info;
543: PetscReal isend[6],irecv[6];
544: PetscInt bs;
548: MatGetBlockSize(A,&bs);
549: if (matis->A->ops->getinfo) {
550: MatGetInfo(matis->A,MAT_LOCAL,&info);
551: isend[0] = info.nz_used;
552: isend[1] = info.nz_allocated;
553: isend[2] = info.nz_unneeded;
554: isend[3] = info.memory;
555: isend[4] = info.mallocs;
556: } else {
557: isend[0] = 0.;
558: isend[1] = 0.;
559: isend[2] = 0.;
560: isend[3] = 0.;
561: isend[4] = 0.;
562: }
563: isend[5] = matis->A->num_ass;
564: if (flag == MAT_LOCAL) {
565: ginfo->nz_used = isend[0];
566: ginfo->nz_allocated = isend[1];
567: ginfo->nz_unneeded = isend[2];
568: ginfo->memory = isend[3];
569: ginfo->mallocs = isend[4];
570: ginfo->assemblies = isend[5];
571: } else if (flag == MAT_GLOBAL_MAX) {
572: MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));
574: ginfo->nz_used = irecv[0];
575: ginfo->nz_allocated = irecv[1];
576: ginfo->nz_unneeded = irecv[2];
577: ginfo->memory = irecv[3];
578: ginfo->mallocs = irecv[4];
579: ginfo->assemblies = irecv[5];
580: } else if (flag == MAT_GLOBAL_SUM) {
581: MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));
583: ginfo->nz_used = irecv[0];
584: ginfo->nz_allocated = irecv[1];
585: ginfo->nz_unneeded = irecv[2];
586: ginfo->memory = irecv[3];
587: ginfo->mallocs = irecv[4];
588: ginfo->assemblies = A->num_ass;
589: }
590: ginfo->block_size = bs;
591: ginfo->fill_ratio_given = 0;
592: ginfo->fill_ratio_needed = 0;
593: ginfo->factor_mallocs = 0;
594: return(0);
595: }
597: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
598: {
599: Mat C,lC,lA;
600: PetscErrorCode ierr;
603: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
604: ISLocalToGlobalMapping rl2g,cl2g;
605: MatCreate(PetscObjectComm((PetscObject)A),&C);
606: MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
607: MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
608: MatSetType(C,MATIS);
609: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
610: MatSetLocalToGlobalMapping(C,cl2g,rl2g);
611: } else {
612: C = *B;
613: }
615: /* perform local transposition */
616: MatISGetLocalMat(A,&lA);
617: MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
618: MatISSetLocalMat(C,lC);
619: MatDestroy(&lC);
621: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
622: *B = C;
623: } else {
624: MatHeaderMerge(A,&C);
625: }
626: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
627: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
628: return(0);
629: }
631: PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
632: {
633: Mat_IS *is = (Mat_IS*)A->data;
637: if (D) { /* MatShift_IS pass D = NULL */
638: VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
639: VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
640: }
641: VecPointwiseDivide(is->y,is->y,is->counter);
642: MatDiagonalSet(is->A,is->y,insmode);
643: return(0);
644: }
646: PetscErrorCode MatShift_IS(Mat A,PetscScalar a)
647: {
648: Mat_IS *is = (Mat_IS*)A->data;
652: VecSet(is->y,a);
653: MatDiagonalSet_IS(A,NULL,ADD_VALUES);
654: return(0);
655: }
657: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
658: {
660: Mat_IS *is = (Mat_IS*)A->data;
661: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
664: #if defined(PETSC_USE_DEBUG)
665: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
666: #endif
667: ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
668: ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
669: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
670: return(0);
671: }
673: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
674: {
676: Mat_IS *is = (Mat_IS*)A->data;
677: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
680: #if defined(PETSC_USE_DEBUG)
681: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
682: #endif
683: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
684: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
685: MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
686: return(0);
687: }
689: static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
690: {
691: PetscInt *owners = map->range;
692: PetscInt n = map->n;
693: PetscSF sf;
694: PetscInt *lidxs,*work = NULL;
695: PetscSFNode *ridxs;
696: PetscMPIInt rank;
697: PetscInt r, p = 0, len = 0;
701: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
702: /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */
703: MPI_Comm_rank(map->comm,&rank);
704: PetscMalloc1(n,&lidxs);
705: for (r = 0; r < n; ++r) lidxs[r] = -1;
706: PetscMalloc1(N,&ridxs);
707: for (r = 0; r < N; ++r) {
708: const PetscInt idx = idxs[r];
709: if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
710: if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
711: PetscLayoutFindOwner(map,idx,&p);
712: }
713: ridxs[r].rank = p;
714: ridxs[r].index = idxs[r] - owners[p];
715: }
716: PetscSFCreate(map->comm,&sf);
717: PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
718: PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
719: PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
720: if (ogidxs) { /* communicate global idxs */
721: PetscInt cum = 0,start,*work2;
723: PetscMalloc1(n,&work);
724: PetscCalloc1(N,&work2);
725: for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
726: MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
727: start -= cum;
728: cum = 0;
729: for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
730: PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);
731: PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);
732: PetscFree(work2);
733: }
734: PetscSFDestroy(&sf);
735: /* Compress and put in indices */
736: for (r = 0; r < n; ++r)
737: if (lidxs[r] >= 0) {
738: if (work) work[len] = work[r];
739: lidxs[len++] = r;
740: }
741: if (on) *on = len;
742: if (oidxs) *oidxs = lidxs;
743: if (ogidxs) *ogidxs = work;
744: return(0);
745: }
747: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
748: {
749: Mat locmat,newlocmat;
750: Mat_IS *newmatis;
751: #if defined(PETSC_USE_DEBUG)
752: Vec rtest,ltest;
753: const PetscScalar *array;
754: #endif
755: const PetscInt *idxs;
756: PetscInt i,m,n;
757: PetscErrorCode ierr;
760: if (scall == MAT_REUSE_MATRIX) {
761: PetscBool ismatis;
763: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
764: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
765: newmatis = (Mat_IS*)(*newmat)->data;
766: if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
767: if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
768: }
769: /* irow and icol may not have duplicate entries */
770: #if defined(PETSC_USE_DEBUG)
771: MatCreateVecs(mat,<est,&rtest);
772: ISGetLocalSize(irow,&n);
773: ISGetIndices(irow,&idxs);
774: for (i=0;i<n;i++) {
775: VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
776: }
777: VecAssemblyBegin(rtest);
778: VecAssemblyEnd(rtest);
779: VecGetLocalSize(rtest,&n);
780: VecGetOwnershipRange(rtest,&m,NULL);
781: VecGetArrayRead(rtest,&array);
782: for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Irow may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
783: VecRestoreArrayRead(rtest,&array);
784: ISRestoreIndices(irow,&idxs);
785: ISGetLocalSize(icol,&n);
786: ISGetIndices(icol,&idxs);
787: for (i=0;i<n;i++) {
788: VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
789: }
790: VecAssemblyBegin(ltest);
791: VecAssemblyEnd(ltest);
792: VecGetLocalSize(ltest,&n);
793: VecGetOwnershipRange(ltest,&m,NULL);
794: VecGetArrayRead(ltest,&array);
795: for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Icol may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
796: VecRestoreArrayRead(ltest,&array);
797: ISRestoreIndices(icol,&idxs);
798: VecDestroy(&rtest);
799: VecDestroy(<est);
800: #endif
801: if (scall == MAT_INITIAL_MATRIX) {
802: Mat_IS *matis = (Mat_IS*)mat->data;
803: ISLocalToGlobalMapping rl2g;
804: IS is;
805: PetscInt *lidxs,*lgidxs,*newgidxs;
806: PetscInt ll,newloc;
807: MPI_Comm comm;
809: PetscObjectGetComm((PetscObject)mat,&comm);
810: ISGetLocalSize(irow,&m);
811: ISGetLocalSize(icol,&n);
812: MatCreate(comm,newmat);
813: MatSetType(*newmat,MATIS);
814: MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
815: /* communicate irow to their owners in the layout */
816: ISGetIndices(irow,&idxs);
817: PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
818: ISRestoreIndices(irow,&idxs);
819: MatISSetUpSF(mat);
820: PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));
821: for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
822: PetscFree(lidxs);
823: PetscFree(lgidxs);
824: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
825: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
826: for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
827: PetscMalloc1(newloc,&newgidxs);
828: PetscMalloc1(newloc,&lidxs);
829: for (i=0,newloc=0;i<matis->sf->nleaves;i++)
830: if (matis->sf_leafdata[i]) {
831: lidxs[newloc] = i;
832: newgidxs[newloc++] = matis->sf_leafdata[i]-1;
833: }
834: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
835: ISLocalToGlobalMappingCreateIS(is,&rl2g);
836: ISDestroy(&is);
837: /* local is to extract local submatrix */
838: newmatis = (Mat_IS*)(*newmat)->data;
839: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
840: if (mat->congruentlayouts == -1) { /* first time we compare rows and cols layouts */
841: PetscBool cong;
843: PetscLayoutCompare(mat->rmap,mat->cmap,&cong);
844: if (cong) mat->congruentlayouts = 1;
845: else mat->congruentlayouts = 0;
846: }
847: if (mat->congruentlayouts && irow == icol && matis->csf == matis->sf) {
848: MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
849: PetscObjectReference((PetscObject)newmatis->getsub_ris);
850: newmatis->getsub_cis = newmatis->getsub_ris;
851: } else {
852: ISLocalToGlobalMapping cl2g;
854: /* communicate icol to their owners in the layout */
855: ISGetIndices(icol,&idxs);
856: PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
857: ISRestoreIndices(icol,&idxs);
858: PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));
859: for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
860: PetscFree(lidxs);
861: PetscFree(lgidxs);
862: PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
863: PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
864: for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
865: PetscMalloc1(newloc,&newgidxs);
866: PetscMalloc1(newloc,&lidxs);
867: for (i=0,newloc=0;i<matis->csf->nleaves;i++)
868: if (matis->csf_leafdata[i]) {
869: lidxs[newloc] = i;
870: newgidxs[newloc++] = matis->csf_leafdata[i]-1;
871: }
872: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
873: ISLocalToGlobalMappingCreateIS(is,&cl2g);
874: ISDestroy(&is);
875: /* local is to extract local submatrix */
876: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
877: MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
878: ISLocalToGlobalMappingDestroy(&cl2g);
879: }
880: ISLocalToGlobalMappingDestroy(&rl2g);
881: } else {
882: MatISGetLocalMat(*newmat,&newlocmat);
883: }
884: MatISGetLocalMat(mat,&locmat);
885: newmatis = (Mat_IS*)(*newmat)->data;
886: MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
887: if (scall == MAT_INITIAL_MATRIX) {
888: MatISSetLocalMat(*newmat,newlocmat);
889: MatDestroy(&newlocmat);
890: }
891: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
892: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
893: return(0);
894: }
896: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
897: {
898: Mat_IS *a = (Mat_IS*)A->data,*b;
899: PetscBool ismatis;
903: PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
904: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
905: b = (Mat_IS*)B->data;
906: MatCopy(a->A,b->A,str);
907: PetscObjectStateIncrease((PetscObject)B);
908: return(0);
909: }
911: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d)
912: {
913: Vec v;
914: const PetscScalar *array;
915: PetscInt i,n;
916: PetscErrorCode ierr;
919: *missing = PETSC_FALSE;
920: MatCreateVecs(A,NULL,&v);
921: MatGetDiagonal(A,v);
922: VecGetLocalSize(v,&n);
923: VecGetArrayRead(v,&array);
924: for (i=0;i<n;i++) if (array[i] == 0.) break;
925: VecRestoreArrayRead(v,&array);
926: VecDestroy(&v);
927: if (i != n) *missing = PETSC_TRUE;
928: if (d) {
929: *d = -1;
930: if (*missing) {
931: PetscInt rstart;
932: MatGetOwnershipRange(A,&rstart,NULL);
933: *d = i+rstart;
934: }
935: }
936: return(0);
937: }
939: static PetscErrorCode MatISSetUpSF_IS(Mat B)
940: {
941: Mat_IS *matis = (Mat_IS*)(B->data);
942: const PetscInt *gidxs;
943: PetscInt nleaves;
947: if (matis->sf) return(0);
948: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
949: ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
950: ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
951: PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
952: ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
953: PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
954: if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
955: ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
956: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
957: ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
958: PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
959: ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
960: PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
961: } else {
962: matis->csf = matis->sf;
963: matis->csf_leafdata = matis->sf_leafdata;
964: matis->csf_rootdata = matis->sf_rootdata;
965: }
966: return(0);
967: }
969: /*@
970: MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.
972: Collective on MPI_Comm
974: Input Parameters:
975: + B - the matrix
976: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
977: (same value is used for all local rows)
978: . d_nnz - array containing the number of nonzeros in the various rows of the
979: DIAGONAL portion of the local submatrix (possibly different for each row)
980: or NULL, if d_nz is used to specify the nonzero structure.
981: The size of this array is equal to the number of local rows, i.e 'm'.
982: For matrices that will be factored, you must leave room for (and set)
983: the diagonal entry even if it is zero.
984: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
985: submatrix (same value is used for all local rows).
986: - o_nnz - array containing the number of nonzeros in the various rows of the
987: OFF-DIAGONAL portion of the local submatrix (possibly different for
988: each row) or NULL, if o_nz is used to specify the nonzero
989: structure. The size of this array is equal to the number
990: of local rows, i.e 'm'.
992: If the *_nnz parameter is given then the *_nz parameter is ignored
994: Level: intermediate
996: Notes: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
997: from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
998: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1000: .keywords: matrix
1002: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1003: @*/
1004: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1005: {
1011: PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1012: return(0);
1013: }
1015: static PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1016: {
1017: Mat_IS *matis = (Mat_IS*)(B->data);
1018: PetscInt bs,i,nlocalcols;
1022: if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
1023: MatISSetUpSF(B);
1025: if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
1026: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
1028: if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz;
1029: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
1031: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1032: MatGetSize(matis->A,NULL,&nlocalcols);
1033: MatGetBlockSize(matis->A,&bs);
1034: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1036: for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1037: MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1038: #if defined(PETSC_HAVE_HYPRE)
1039: MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1040: #endif
1042: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
1043: MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1045: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i]-i;
1046: MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1048: /* for other matrix types */
1049: MatSetUp(matis->A);
1050: return(0);
1051: }
1053: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1054: {
1055: Mat_IS *matis = (Mat_IS*)(A->data);
1056: PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1057: const PetscInt *global_indices_r,*global_indices_c;
1058: PetscInt i,j,bs,rows,cols;
1059: PetscInt lrows,lcols;
1060: PetscInt local_rows,local_cols;
1061: PetscMPIInt nsubdomains;
1062: PetscBool isdense,issbaij;
1063: PetscErrorCode ierr;
1066: MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);
1067: MatGetSize(A,&rows,&cols);
1068: MatGetBlockSize(A,&bs);
1069: MatGetSize(matis->A,&local_rows,&local_cols);
1070: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1071: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1072: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
1073: if (A->rmap->mapping != A->cmap->mapping) {
1074: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
1075: } else {
1076: global_indices_c = global_indices_r;
1077: }
1079: if (issbaij) {
1080: MatGetRowUpperTriangular(matis->A);
1081: }
1082: /*
1083: An SF reduce is needed to sum up properly on shared rows.
1084: Note that generally preallocation is not exact, since it overestimates nonzeros
1085: */
1086: MatISSetUpSF(A);
1087: MatGetLocalSize(A,&lrows,&lcols);
1088: MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1089: /* All processes need to compute entire row ownership */
1090: PetscMalloc1(rows,&row_ownership);
1091: MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1092: for (i=0;i<nsubdomains;i++) {
1093: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
1094: row_ownership[j] = i;
1095: }
1096: }
1097: MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);
1099: /*
1100: my_dnz and my_onz contains exact contribution to preallocation from each local mat
1101: then, they will be summed up properly. This way, preallocation is always sufficient
1102: */
1103: PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1104: /* preallocation as a MATAIJ */
1105: if (isdense) { /* special case for dense local matrices */
1106: for (i=0;i<local_rows;i++) {
1107: PetscInt owner = row_ownership[global_indices_r[i]];
1108: for (j=0;j<local_cols;j++) {
1109: PetscInt index_col = global_indices_c[j];
1110: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1111: my_dnz[i] += 1;
1112: } else { /* offdiag block */
1113: my_onz[i] += 1;
1114: }
1115: }
1116: }
1117: } else if (matis->A->ops->getrowij) {
1118: const PetscInt *ii,*jj,*jptr;
1119: PetscBool done;
1120: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1121: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1122: jptr = jj;
1123: for (i=0;i<local_rows;i++) {
1124: PetscInt index_row = global_indices_r[i];
1125: for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1126: PetscInt owner = row_ownership[index_row];
1127: PetscInt index_col = global_indices_c[*jptr];
1128: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1129: my_dnz[i] += 1;
1130: } else { /* offdiag block */
1131: my_onz[i] += 1;
1132: }
1133: /* same as before, interchanging rows and cols */
1134: if (issbaij && index_col != index_row) {
1135: owner = row_ownership[index_col];
1136: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1137: my_dnz[*jptr] += 1;
1138: } else {
1139: my_onz[*jptr] += 1;
1140: }
1141: }
1142: }
1143: }
1144: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1145: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1146: } else { /* loop over rows and use MatGetRow */
1147: for (i=0;i<local_rows;i++) {
1148: const PetscInt *cols;
1149: PetscInt ncols,index_row = global_indices_r[i];
1150: MatGetRow(matis->A,i,&ncols,&cols,NULL);
1151: for (j=0;j<ncols;j++) {
1152: PetscInt owner = row_ownership[index_row];
1153: PetscInt index_col = global_indices_c[cols[j]];
1154: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1155: my_dnz[i] += 1;
1156: } else { /* offdiag block */
1157: my_onz[i] += 1;
1158: }
1159: /* same as before, interchanging rows and cols */
1160: if (issbaij && index_col != index_row) {
1161: owner = row_ownership[index_col];
1162: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1163: my_dnz[cols[j]] += 1;
1164: } else {
1165: my_onz[cols[j]] += 1;
1166: }
1167: }
1168: }
1169: MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1170: }
1171: }
1172: if (global_indices_c != global_indices_r) {
1173: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1174: }
1175: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1176: PetscFree(row_ownership);
1178: /* Reduce my_dnz and my_onz */
1179: if (maxreduce) {
1180: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1181: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1182: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1183: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1184: } else {
1185: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1186: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1187: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1188: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1189: }
1190: PetscFree2(my_dnz,my_onz);
1192: /* Resize preallocation if overestimated */
1193: for (i=0;i<lrows;i++) {
1194: dnz[i] = PetscMin(dnz[i],lcols);
1195: onz[i] = PetscMin(onz[i],cols-lcols);
1196: }
1198: /* Set preallocation */
1199: MatSeqAIJSetPreallocation(B,0,dnz);
1200: MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1201: for (i=0;i<lrows/bs;i++) {
1202: dnz[i] = dnz[i*bs]/bs;
1203: onz[i] = onz[i*bs]/bs;
1204: }
1205: MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1206: MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1207: MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1208: MatPreallocateFinalize(dnz,onz);
1209: if (issbaij) {
1210: MatRestoreRowUpperTriangular(matis->A);
1211: }
1212: return(0);
1213: }
1215: static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
1216: {
1217: Mat_IS *matis = (Mat_IS*)(mat->data);
1218: Mat local_mat;
1219: /* info on mat */
1220: PetscInt bs,rows,cols,lrows,lcols;
1221: PetscInt local_rows,local_cols;
1222: PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij;
1223: #if defined (PETSC_USE_DEBUG)
1224: PetscBool lb[4],bb[4];
1225: #endif
1226: PetscMPIInt nsubdomains;
1227: /* values insertion */
1228: PetscScalar *array;
1229: /* work */
1233: /* get info from mat */
1234: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
1235: if (nsubdomains == 1) {
1236: Mat B;
1237: IS rows,cols;
1238: IS irows,icols;
1239: const PetscInt *ridxs,*cidxs;
1241: ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);
1242: ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);
1243: ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);
1244: ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);
1245: ISSetPermutation(rows);
1246: ISSetPermutation(cols);
1247: ISInvertPermutation(rows,mat->rmap->n,&irows);
1248: ISInvertPermutation(cols,mat->cmap->n,&icols);
1249: ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);
1250: ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);
1251: ISDestroy(&cols);
1252: ISDestroy(&rows);
1253: MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
1254: MatCreateSubMatrix(B,irows,icols,reuse,M);
1255: MatDestroy(&B);
1256: ISDestroy(&icols);
1257: ISDestroy(&irows);
1258: return(0);
1259: }
1260: MatGetSize(mat,&rows,&cols);
1261: MatGetBlockSize(mat,&bs);
1262: MatGetLocalSize(mat,&lrows,&lcols);
1263: MatGetSize(matis->A,&local_rows,&local_cols);
1264: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
1265: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
1266: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
1267: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
1268: if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
1269: #if defined (PETSC_USE_DEBUG)
1270: lb[0] = isseqdense;
1271: lb[1] = isseqaij;
1272: lb[2] = isseqbaij;
1273: lb[3] = isseqsbaij;
1274: MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
1275: if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
1276: #endif
1278: if (reuse == MAT_INITIAL_MATRIX) {
1279: MatCreate(PetscObjectComm((PetscObject)mat),M);
1280: MatSetSizes(*M,lrows,lcols,rows,cols);
1281: MatSetBlockSize(*M,bs);
1282: if (!isseqsbaij) {
1283: MatSetType(*M,MATAIJ);
1284: } else {
1285: MatSetType(*M,MATSBAIJ);
1286: }
1287: MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);
1288: } else {
1289: PetscInt mbs,mrows,mcols,mlrows,mlcols;
1290: /* some checks */
1291: MatGetBlockSize(*M,&mbs);
1292: MatGetSize(*M,&mrows,&mcols);
1293: MatGetLocalSize(*M,&mlrows,&mlcols);
1294: if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
1295: if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
1296: if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
1297: if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
1298: if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
1299: MatZeroEntries(*M);
1300: }
1302: if (isseqsbaij) {
1303: MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
1304: } else {
1305: PetscObjectReference((PetscObject)matis->A);
1306: local_mat = matis->A;
1307: }
1309: /* Set values */
1310: MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);
1311: if (isseqdense) { /* special case for dense local matrices */
1312: PetscInt i,*dummy;
1314: PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
1315: for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
1316: MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
1317: MatDenseGetArray(local_mat,&array);
1318: MatSetValuesLocal(*M,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
1319: MatDenseRestoreArray(local_mat,&array);
1320: PetscFree(dummy);
1321: } else if (isseqaij) {
1322: PetscInt i,nvtxs,*xadj,*adjncy;
1323: PetscBool done;
1325: MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1326: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1327: MatSeqAIJGetArray(local_mat,&array);
1328: for (i=0;i<nvtxs;i++) {
1329: MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
1330: }
1331: MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1332: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1333: MatSeqAIJRestoreArray(local_mat,&array);
1334: } else { /* very basic values insertion for all other matrix types */
1335: PetscInt i;
1337: for (i=0;i<local_rows;i++) {
1338: PetscInt j;
1339: const PetscInt *local_indices_cols;
1341: MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
1342: MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);
1343: MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
1344: }
1345: }
1346: MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
1347: MatDestroy(&local_mat);
1348: MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
1349: if (isseqdense) {
1350: MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
1351: }
1352: return(0);
1353: }
1355: /*@
1356: MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format
1358: Input Parameter:
1359: . mat - the matrix (should be of type MATIS)
1360: . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
1362: Output Parameter:
1363: . newmat - the matrix in AIJ format
1365: Level: developer
1367: Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested.
1369: .seealso: MATIS
1370: @*/
1371: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
1372: {
1379: if (reuse != MAT_INITIAL_MATRIX) {
1382: if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
1383: }
1384: PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
1385: return(0);
1386: }
1388: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
1389: {
1391: Mat_IS *matis = (Mat_IS*)(mat->data);
1392: PetscInt bs,m,n,M,N;
1393: Mat B,localmat;
1396: MatGetBlockSize(mat,&bs);
1397: MatGetSize(mat,&M,&N);
1398: MatGetLocalSize(mat,&m,&n);
1399: MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);
1400: MatDuplicate(matis->A,op,&localmat);
1401: MatISSetLocalMat(B,localmat);
1402: MatDestroy(&localmat);
1403: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1404: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1405: *newmat = B;
1406: return(0);
1407: }
1409: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
1410: {
1412: Mat_IS *matis = (Mat_IS*)A->data;
1413: PetscBool local_sym;
1416: MatIsHermitian(matis->A,tol,&local_sym);
1417: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1418: return(0);
1419: }
1421: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
1422: {
1424: Mat_IS *matis = (Mat_IS*)A->data;
1425: PetscBool local_sym;
1428: MatIsSymmetric(matis->A,tol,&local_sym);
1429: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1430: return(0);
1431: }
1433: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg)
1434: {
1436: Mat_IS *matis = (Mat_IS*)A->data;
1437: PetscBool local_sym;
1440: if (A->rmap->mapping != A->cmap->mapping) {
1441: *flg = PETSC_FALSE;
1442: return(0);
1443: }
1444: MatIsStructurallySymmetric(matis->A,&local_sym);
1445: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1446: return(0);
1447: }
1449: static PetscErrorCode MatDestroy_IS(Mat A)
1450: {
1452: Mat_IS *b = (Mat_IS*)A->data;
1455: MatDestroy(&b->A);
1456: VecScatterDestroy(&b->cctx);
1457: VecScatterDestroy(&b->rctx);
1458: VecDestroy(&b->x);
1459: VecDestroy(&b->y);
1460: VecDestroy(&b->counter);
1461: ISDestroy(&b->getsub_ris);
1462: ISDestroy(&b->getsub_cis);
1463: if (b->sf != b->csf) {
1464: PetscSFDestroy(&b->csf);
1465: PetscFree2(b->csf_rootdata,b->csf_leafdata);
1466: }
1467: PetscSFDestroy(&b->sf);
1468: PetscFree2(b->sf_rootdata,b->sf_leafdata);
1469: PetscFree(A->data);
1470: PetscObjectChangeTypeName((PetscObject)A,0);
1471: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
1472: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
1473: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
1474: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
1475: PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);
1476: return(0);
1477: }
1479: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
1480: {
1482: Mat_IS *is = (Mat_IS*)A->data;
1483: PetscScalar zero = 0.0;
1486: /* scatter the global vector x into the local work vector */
1487: VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
1488: VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
1490: /* multiply the local matrix */
1491: MatMult(is->A,is->x,is->y);
1493: /* scatter product back into global memory */
1494: VecSet(y,zero);
1495: VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
1496: VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
1497: return(0);
1498: }
1500: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
1501: {
1502: Vec temp_vec;
1506: if (v3 != v2) {
1507: MatMult(A,v1,v3);
1508: VecAXPY(v3,1.0,v2);
1509: } else {
1510: VecDuplicate(v2,&temp_vec);
1511: MatMult(A,v1,temp_vec);
1512: VecAXPY(temp_vec,1.0,v2);
1513: VecCopy(temp_vec,v3);
1514: VecDestroy(&temp_vec);
1515: }
1516: return(0);
1517: }
1519: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
1520: {
1521: Mat_IS *is = (Mat_IS*)A->data;
1525: /* scatter the global vector x into the local work vector */
1526: VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
1527: VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
1529: /* multiply the local matrix */
1530: MatMultTranspose(is->A,is->y,is->x);
1532: /* scatter product back into global vector */
1533: VecSet(x,0);
1534: VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
1535: VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
1536: return(0);
1537: }
1539: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
1540: {
1541: Vec temp_vec;
1545: if (v3 != v2) {
1546: MatMultTranspose(A,v1,v3);
1547: VecAXPY(v3,1.0,v2);
1548: } else {
1549: VecDuplicate(v2,&temp_vec);
1550: MatMultTranspose(A,v1,temp_vec);
1551: VecAXPY(temp_vec,1.0,v2);
1552: VecCopy(temp_vec,v3);
1553: VecDestroy(&temp_vec);
1554: }
1555: return(0);
1556: }
1558: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
1559: {
1560: Mat_IS *a = (Mat_IS*)A->data;
1562: PetscViewer sviewer;
1563: PetscBool isascii,view = PETSC_TRUE;
1566: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1567: if (isascii) {
1568: PetscViewerFormat format;
1570: PetscViewerGetFormat(viewer,&format);
1571: if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
1572: }
1573: if (!view) return(0);
1574: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
1575: MatView(a->A,sviewer);
1576: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
1577: PetscViewerFlush(viewer);
1578: return(0);
1579: }
1581: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1582: {
1584: PetscInt nr,rbs,nc,cbs;
1585: Mat_IS *is = (Mat_IS*)A->data;
1586: IS from,to;
1587: Vec cglobal,rglobal;
1592: /* Destroy any previous data */
1593: VecDestroy(&is->x);
1594: VecDestroy(&is->y);
1595: VecDestroy(&is->counter);
1596: VecScatterDestroy(&is->rctx);
1597: VecScatterDestroy(&is->cctx);
1598: MatDestroy(&is->A);
1599: if (is->csf != is->sf) {
1600: PetscSFDestroy(&is->csf);
1601: PetscFree2(is->csf_rootdata,is->csf_leafdata);
1602: }
1603: PetscSFDestroy(&is->sf);
1604: PetscFree2(is->sf_rootdata,is->sf_leafdata);
1606: /* Setup Layout and set local to global maps */
1607: PetscLayoutSetUp(A->rmap);
1608: PetscLayoutSetUp(A->cmap);
1609: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
1610: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
1612: /* Create the local matrix A */
1613: ISLocalToGlobalMappingGetSize(rmapping,&nr);
1614: ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
1615: ISLocalToGlobalMappingGetSize(cmapping,&nc);
1616: ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
1617: MatCreate(PETSC_COMM_SELF,&is->A);
1618: MatSetType(is->A,MATAIJ);
1619: MatSetSizes(is->A,nr,nc,nr,nc);
1620: MatSetBlockSizes(is->A,rbs,cbs);
1621: MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
1622: MatAppendOptionsPrefix(is->A,"is_");
1623: MatSetFromOptions(is->A);
1624: PetscLayoutSetUp(is->A->rmap);
1625: PetscLayoutSetUp(is->A->cmap);
1627: if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
1628: /* Create the local work vectors */
1629: MatCreateVecs(is->A,&is->x,&is->y);
1631: /* setup the global to local scatters */
1632: MatCreateVecs(A,&cglobal,&rglobal);
1633: ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);
1634: ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
1635: VecScatterCreate(rglobal,from,is->y,to,&is->rctx);
1636: if (rmapping != cmapping) {
1637: ISDestroy(&to);
1638: ISDestroy(&from);
1639: ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);
1640: ISLocalToGlobalMappingApplyIS(cmapping,to,&from);
1641: VecScatterCreate(cglobal,from,is->x,to,&is->cctx);
1642: } else {
1643: PetscObjectReference((PetscObject)is->rctx);
1644: is->cctx = is->rctx;
1645: }
1647: /* interface counter vector (local) */
1648: VecDuplicate(is->y,&is->counter);
1649: VecSet(is->y,1.);
1650: VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
1651: VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
1652: VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
1653: VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
1655: /* free workspace */
1656: VecDestroy(&rglobal);
1657: VecDestroy(&cglobal);
1658: ISDestroy(&to);
1659: ISDestroy(&from);
1660: }
1661: MatSetUp(A);
1662: return(0);
1663: }
1665: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1666: {
1667: Mat_IS *is = (Mat_IS*)mat->data;
1669: #if defined(PETSC_USE_DEBUG)
1670: PetscInt i,zm,zn;
1671: #endif
1672: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1675: #if defined(PETSC_USE_DEBUG)
1676: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
1677: /* count negative indices */
1678: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
1679: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
1680: #endif
1681: ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
1682: ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
1683: #if defined(PETSC_USE_DEBUG)
1684: /* count negative indices (should be the same as before) */
1685: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
1686: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
1687: /* disable check when usesetlocal is true */
1688: if (!is->usesetlocal && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
1689: if (!is->usesetlocal && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
1690: #endif
1691: if (is->usesetlocal) {
1692: MatSetValuesLocal(is->A,m,rows_l,n,cols_l,values,addv);
1693: } else {
1694: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
1695: }
1696: return(0);
1697: }
1699: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1700: {
1701: Mat_IS *is = (Mat_IS*)mat->data;
1703: #if defined(PETSC_USE_DEBUG)
1704: PetscInt i,zm,zn;
1705: #endif
1706: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1709: #if defined(PETSC_USE_DEBUG)
1710: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
1711: /* count negative indices */
1712: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
1713: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
1714: #endif
1715: ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
1716: ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
1717: #if defined(PETSC_USE_DEBUG)
1718: /* count negative indices (should be the same as before) */
1719: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
1720: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
1721: if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
1722: if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
1723: #endif
1724: MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
1725: return(0);
1726: }
1728: static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1729: {
1731: Mat_IS *is = (Mat_IS*)A->data;
1734: if (is->usesetlocal) {
1735: MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
1736: } else {
1737: MatSetValues(is->A,m,rows,n,cols,values,addv);
1738: }
1739: return(0);
1740: }
1742: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1743: {
1745: Mat_IS *is = (Mat_IS*)A->data;
1748: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
1749: return(0);
1750: }
1752: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns)
1753: {
1754: Mat_IS *is = (Mat_IS*)A->data;
1758: if (!n) {
1759: is->pure_neumann = PETSC_TRUE;
1760: } else {
1761: PetscInt i;
1762: is->pure_neumann = PETSC_FALSE;
1764: if (columns) {
1765: MatZeroRowsColumns(is->A,n,rows,diag,0,0);
1766: } else {
1767: MatZeroRows(is->A,n,rows,diag,0,0);
1768: }
1769: if (diag != 0.) {
1770: const PetscScalar *array;
1771: VecGetArrayRead(is->counter,&array);
1772: for (i=0; i<n; i++) {
1773: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
1774: }
1775: VecRestoreArrayRead(is->counter,&array);
1776: }
1777: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
1778: MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
1779: }
1780: return(0);
1781: }
1783: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
1784: {
1785: Mat_IS *matis = (Mat_IS*)A->data;
1786: PetscInt nr,nl,len,i;
1787: PetscInt *lrows;
1791: #if defined(PETSC_USE_DEBUG)
1792: if (columns || diag != 0. || (x && b)) {
1793: PetscBool cong;
1795: PetscLayoutCompare(A->rmap,A->cmap,&cong);
1796: cong = (PetscBool)(cong && matis->sf == matis->csf);
1797: if (!cong && columns) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
1798: if (!cong && diag != 0.) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
1799: if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
1800: }
1801: #endif
1802: /* get locally owned rows */
1803: PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);
1804: /* fix right hand side if needed */
1805: if (x && b) {
1806: const PetscScalar *xx;
1807: PetscScalar *bb;
1809: VecGetArrayRead(x, &xx);
1810: VecGetArray(b, &bb);
1811: for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
1812: VecRestoreArrayRead(x, &xx);
1813: VecRestoreArray(b, &bb);
1814: }
1815: /* get rows associated to the local matrices */
1816: MatISSetUpSF(A);
1817: MatGetSize(matis->A,&nl,NULL);
1818: PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));
1819: PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));
1820: for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
1821: PetscFree(lrows);
1822: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1823: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1824: PetscMalloc1(nl,&lrows);
1825: for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
1826: MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
1827: PetscFree(lrows);
1828: return(0);
1829: }
1831: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1832: {
1836: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
1837: return(0);
1838: }
1840: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1841: {
1845: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
1846: return(0);
1847: }
1849: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
1850: {
1851: Mat_IS *is = (Mat_IS*)A->data;
1855: MatAssemblyBegin(is->A,type);
1856: return(0);
1857: }
1859: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
1860: {
1861: Mat_IS *is = (Mat_IS*)A->data;
1865: MatAssemblyEnd(is->A,type);
1866: /* fix for local empty rows/cols */
1867: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
1868: Mat newlA;
1869: ISLocalToGlobalMapping l2g;
1870: IS tis;
1871: const PetscScalar *v;
1872: PetscInt i,n,cf,*loce,*locf;
1873: PetscBool sym;
1875: MatGetRowMaxAbs(is->A,is->y,NULL);
1876: MatIsSymmetric(is->A,PETSC_SMALL,&sym);
1877: if (!sym) SETERRQ(PetscObjectComm((PetscObject)is->A),PETSC_ERR_SUP,"Not yet implemented for unsymmetric case");
1878: VecGetLocalSize(is->y,&n);
1879: PetscMalloc1(n,&loce);
1880: PetscMalloc1(n,&locf);
1881: VecGetArrayRead(is->y,&v);
1882: for (i=0,cf=0;i<n;i++) {
1883: if (v[i] == 0.0) {
1884: loce[i] = -1;
1885: } else {
1886: loce[i] = cf;
1887: locf[cf++] = i;
1888: }
1889: }
1890: VecRestoreArrayRead(is->y,&v);
1891: /* extract valid submatrix */
1892: ISCreateGeneral(PetscObjectComm((PetscObject)is->A),cf,locf,PETSC_USE_POINTER,&tis);
1893: MatCreateSubMatrix(is->A,tis,tis,MAT_INITIAL_MATRIX,&newlA);
1894: ISDestroy(&tis);
1895: /* attach local l2g map for successive calls of MatSetValues */
1896: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is->A),1,n,loce,PETSC_OWN_POINTER,&l2g);
1897: MatSetLocalToGlobalMapping(newlA,l2g,l2g);
1898: ISLocalToGlobalMappingDestroy(&l2g);
1899: /* flag MatSetValues */
1900: is->usesetlocal = PETSC_TRUE;
1901: /* attach new global l2g map */
1902: ISLocalToGlobalMappingApply(A->rmap->mapping,cf,locf,locf);
1903: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,cf,locf,PETSC_OWN_POINTER,&l2g);
1904: MatSetLocalToGlobalMapping(A,l2g,l2g);
1905: MatISSetLocalMat(A,newlA);
1906: MatDestroy(&newlA);
1907: ISLocalToGlobalMappingDestroy(&l2g);
1908: }
1909: is->locempty = PETSC_FALSE;
1910: return(0);
1911: }
1913: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
1914: {
1915: Mat_IS *is = (Mat_IS*)mat->data;
1918: *local = is->A;
1919: return(0);
1920: }
1922: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
1923: {
1925: *local = NULL;
1926: return(0);
1927: }
1929: /*@
1930: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
1932: Input Parameter:
1933: . mat - the matrix
1935: Output Parameter:
1936: . local - the local matrix
1938: Level: advanced
1940: Notes:
1941: This can be called if you have precomputed the nonzero structure of the
1942: matrix and want to provide it to the inner matrix object to improve the performance
1943: of the MatSetValues() operation.
1945: Call MatISRestoreLocalMat() when finished with the local matrix.
1947: .seealso: MATIS
1948: @*/
1949: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
1950: {
1956: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
1957: return(0);
1958: }
1960: /*@
1961: MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()
1963: Input Parameter:
1964: . mat - the matrix
1966: Output Parameter:
1967: . local - the local matrix
1969: Level: advanced
1971: .seealso: MATIS
1972: @*/
1973: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
1974: {
1980: PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
1981: return(0);
1982: }
1984: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
1985: {
1986: Mat_IS *is = (Mat_IS*)mat->data;
1987: PetscInt nrows,ncols,orows,ocols;
1991: if (is->A) {
1992: MatGetSize(is->A,&orows,&ocols);
1993: MatGetSize(local,&nrows,&ncols);
1994: if (orows != nrows || ocols != ncols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local MATIS matrix should be of size %Dx%D (you passed a %Dx%D matrix)",orows,ocols,nrows,ncols);
1995: }
1996: PetscObjectReference((PetscObject)local);
1997: MatDestroy(&is->A);
1998: is->A = local;
1999: return(0);
2000: }
2002: /*@
2003: MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.
2005: Input Parameter:
2006: . mat - the matrix
2007: . local - the local matrix
2009: Output Parameter:
2011: Level: advanced
2013: Notes:
2014: This can be called if you have precomputed the local matrix and
2015: want to provide it to the matrix object MATIS.
2017: .seealso: MATIS
2018: @*/
2019: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
2020: {
2026: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
2027: return(0);
2028: }
2030: static PetscErrorCode MatZeroEntries_IS(Mat A)
2031: {
2032: Mat_IS *a = (Mat_IS*)A->data;
2036: MatZeroEntries(a->A);
2037: return(0);
2038: }
2040: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
2041: {
2042: Mat_IS *is = (Mat_IS*)A->data;
2046: MatScale(is->A,a);
2047: return(0);
2048: }
2050: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2051: {
2052: Mat_IS *is = (Mat_IS*)A->data;
2056: /* get diagonal of the local matrix */
2057: MatGetDiagonal(is->A,is->y);
2059: /* scatter diagonal back into global vector */
2060: VecSet(v,0);
2061: VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2062: VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2063: return(0);
2064: }
2066: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
2067: {
2068: Mat_IS *a = (Mat_IS*)A->data;
2072: MatSetOption(a->A,op,flg);
2073: return(0);
2074: }
2076: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
2077: {
2078: Mat_IS *y = (Mat_IS*)Y->data;
2079: Mat_IS *x;
2080: #if defined(PETSC_USE_DEBUG)
2081: PetscBool ismatis;
2082: #endif
2086: #if defined(PETSC_USE_DEBUG)
2087: PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
2088: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
2089: #endif
2090: x = (Mat_IS*)X->data;
2091: MatAXPY(y->A,a,x->A,str);
2092: return(0);
2093: }
2095: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
2096: {
2097: Mat lA;
2098: Mat_IS *matis;
2099: ISLocalToGlobalMapping rl2g,cl2g;
2100: IS is;
2101: const PetscInt *rg,*rl;
2102: PetscInt nrg;
2103: PetscInt N,M,nrl,i,*idxs;
2104: PetscErrorCode ierr;
2107: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
2108: ISGetLocalSize(row,&nrl);
2109: ISGetIndices(row,&rl);
2110: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
2111: #if defined(PETSC_USE_DEBUG)
2112: for (i=0;i<nrl;i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater than maximum possible %D",i,rl[i],nrg);
2113: #endif
2114: PetscMalloc1(nrg,&idxs);
2115: /* map from [0,nrl) to row */
2116: for (i=0;i<nrl;i++) idxs[i] = rl[i];
2117: #if defined(PETSC_USE_DEBUG)
2118: for (i=nrl;i<nrg;i++) idxs[i] = nrg;
2119: #else
2120: for (i=nrl;i<nrg;i++) idxs[i] = -1;
2121: #endif
2122: ISRestoreIndices(row,&rl);
2123: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
2124: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
2125: ISLocalToGlobalMappingCreateIS(is,&rl2g);
2126: ISDestroy(&is);
2127: /* compute new l2g map for columns */
2128: if (col != row || A->rmap->mapping != A->cmap->mapping) {
2129: const PetscInt *cg,*cl;
2130: PetscInt ncg;
2131: PetscInt ncl;
2133: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
2134: ISGetLocalSize(col,&ncl);
2135: ISGetIndices(col,&cl);
2136: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
2137: #if defined(PETSC_USE_DEBUG)
2138: for (i=0;i<ncl;i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater than maximum possible %D",i,cl[i],ncg);
2139: #endif
2140: PetscMalloc1(ncg,&idxs);
2141: /* map from [0,ncl) to col */
2142: for (i=0;i<ncl;i++) idxs[i] = cl[i];
2143: #if defined(PETSC_USE_DEBUG)
2144: for (i=ncl;i<ncg;i++) idxs[i] = ncg;
2145: #else
2146: for (i=ncl;i<ncg;i++) idxs[i] = -1;
2147: #endif
2148: ISRestoreIndices(col,&cl);
2149: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
2150: ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
2151: ISLocalToGlobalMappingCreateIS(is,&cl2g);
2152: ISDestroy(&is);
2153: } else {
2154: PetscObjectReference((PetscObject)rl2g);
2155: cl2g = rl2g;
2156: }
2157: /* create the MATIS submatrix */
2158: MatGetSize(A,&M,&N);
2159: MatCreate(PetscObjectComm((PetscObject)A),submat);
2160: MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
2161: MatSetType(*submat,MATIS);
2162: matis = (Mat_IS*)((*submat)->data);
2163: matis->islocalref = PETSC_TRUE;
2164: MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
2165: MatISGetLocalMat(A,&lA);
2166: MatISSetLocalMat(*submat,lA);
2167: MatSetUp(*submat);
2168: MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
2169: MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
2170: ISLocalToGlobalMappingDestroy(&rl2g);
2171: ISLocalToGlobalMappingDestroy(&cl2g);
2172: /* remove unsupported ops */
2173: PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
2174: (*submat)->ops->destroy = MatDestroy_IS;
2175: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
2176: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
2177: (*submat)->ops->assemblybegin = MatAssemblyBegin_IS;
2178: (*submat)->ops->assemblyend = MatAssemblyEnd_IS;
2179: return(0);
2180: }
2182: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
2183: {
2184: Mat_IS *a = (Mat_IS*)A->data;
2188: PetscOptionsHead(PetscOptionsObject,"MATIS options");
2189: PetscObjectOptionsBegin((PetscObject)A);
2190: PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns",NULL,a->locempty,&a->locempty,NULL);
2191: PetscOptionsEnd();
2192: return(0);
2193: }
2195: /*@
2196: MatCreateIS - Creates a "process" unassembled matrix, assembled on each
2197: process but not across processes.
2199: Input Parameters:
2200: + comm - MPI communicator that will share the matrix
2201: . bs - block size of the matrix
2202: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
2203: . rmap - local to global map for rows
2204: - cmap - local to global map for cols
2206: Output Parameter:
2207: . A - the resulting matrix
2209: Level: advanced
2211: Notes: See MATIS for more details.
2212: m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors
2213: used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
2214: If either rmap or cmap are NULL, then the matrix is assumed to be square.
2216: .seealso: MATIS, MatSetLocalToGlobalMapping()
2217: @*/
2218: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
2219: {
2223: if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
2224: MatCreate(comm,A);
2225: MatSetSizes(*A,m,n,M,N);
2226: if (bs > 0) {
2227: MatSetBlockSize(*A,bs);
2228: }
2229: MatSetType(*A,MATIS);
2230: if (rmap && cmap) {
2231: MatSetLocalToGlobalMapping(*A,rmap,cmap);
2232: } else if (!rmap) {
2233: MatSetLocalToGlobalMapping(*A,cmap,cmap);
2234: } else {
2235: MatSetLocalToGlobalMapping(*A,rmap,rmap);
2236: }
2237: return(0);
2238: }
2240: /*MC
2241: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP).
2242: This stores the matrices in globally unassembled form. Each processor
2243: assembles only its local Neumann problem and the parallel matrix vector
2244: product is handled "implicitly".
2246: Operations Provided:
2247: + MatMult()
2248: . MatMultAdd()
2249: . MatMultTranspose()
2250: . MatMultTransposeAdd()
2251: . MatZeroEntries()
2252: . MatSetOption()
2253: . MatZeroRows()
2254: . MatSetValues()
2255: . MatSetValuesBlocked()
2256: . MatSetValuesLocal()
2257: . MatSetValuesBlockedLocal()
2258: . MatScale()
2259: . MatGetDiagonal()
2260: . MatMissingDiagonal()
2261: . MatDuplicate()
2262: . MatCopy()
2263: . MatAXPY()
2264: . MatCreateSubMatrix()
2265: . MatGetLocalSubMatrix()
2266: . MatTranspose()
2267: - MatSetLocalToGlobalMapping()
2269: Options Database Keys:
2270: . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
2272: Notes: Options prefix for the inner matrix are given by -is_mat_xxx
2274: You must call MatSetLocalToGlobalMapping() before using this matrix type.
2276: You can do matrix preallocation on the local matrix after you obtain it with
2277: MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation()
2279: Level: advanced
2281: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP
2283: M*/
2285: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
2286: {
2288: Mat_IS *b;
2291: PetscNewLog(A,&b);
2292: A->data = (void*)b;
2294: /* matrix ops */
2295: PetscMemzero(A->ops,sizeof(struct _MatOps));
2296: A->ops->mult = MatMult_IS;
2297: A->ops->multadd = MatMultAdd_IS;
2298: A->ops->multtranspose = MatMultTranspose_IS;
2299: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
2300: A->ops->destroy = MatDestroy_IS;
2301: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
2302: A->ops->setvalues = MatSetValues_IS;
2303: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
2304: A->ops->setvalueslocal = MatSetValuesLocal_IS;
2305: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
2306: A->ops->zerorows = MatZeroRows_IS;
2307: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
2308: A->ops->assemblybegin = MatAssemblyBegin_IS;
2309: A->ops->assemblyend = MatAssemblyEnd_IS;
2310: A->ops->view = MatView_IS;
2311: A->ops->zeroentries = MatZeroEntries_IS;
2312: A->ops->scale = MatScale_IS;
2313: A->ops->getdiagonal = MatGetDiagonal_IS;
2314: A->ops->setoption = MatSetOption_IS;
2315: A->ops->ishermitian = MatIsHermitian_IS;
2316: A->ops->issymmetric = MatIsSymmetric_IS;
2317: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
2318: A->ops->duplicate = MatDuplicate_IS;
2319: A->ops->missingdiagonal = MatMissingDiagonal_IS;
2320: A->ops->copy = MatCopy_IS;
2321: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
2322: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
2323: A->ops->axpy = MatAXPY_IS;
2324: A->ops->diagonalset = MatDiagonalSet_IS;
2325: A->ops->shift = MatShift_IS;
2326: A->ops->transpose = MatTranspose_IS;
2327: A->ops->getinfo = MatGetInfo_IS;
2328: A->ops->diagonalscale = MatDiagonalScale_IS;
2329: A->ops->setfromoptions = MatSetFromOptions_IS;
2331: /* special MATIS functions */
2332: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
2333: PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
2334: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
2335: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
2336: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
2337: PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);
2338: PetscObjectChangeTypeName((PetscObject)A,MATIS);
2339: return(0);
2340: }