Actual source code: matis.c
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 <petsc/private/sfimpl.h>
12: #include <petsc/private/vecimpl.h>
14: #define MATIS_MAX_ENTRIES_INSERTION 2048
15: static PetscErrorCode MatSetValuesLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode);
16: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode);
17: static PetscErrorCode MatISSetUpScatters_Private(Mat);
19: static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr)
20: {
21: MatISPtAP ptap = (MatISPtAP)ptr;
25: MatDestroySubMatrices(ptap->ris1 ? 2 : 1,&ptap->lP);
26: ISDestroy(&ptap->cis0);
27: ISDestroy(&ptap->cis1);
28: ISDestroy(&ptap->ris0);
29: ISDestroy(&ptap->ris1);
30: PetscFree(ptap);
31: return(0);
32: }
34: static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
35: {
36: MatISPtAP ptap;
37: Mat_IS *matis = (Mat_IS*)A->data;
38: Mat lA,lC;
39: MatReuse reuse;
40: IS ris[2],cis[2];
41: PetscContainer c;
42: PetscInt n;
46: PetscObjectQuery((PetscObject)C,"_MatIS_PtAP",(PetscObject*)&c);
47: if (!c) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing PtAP information");
48: PetscContainerGetPointer(c,(void**)&ptap);
49: ris[0] = ptap->ris0;
50: ris[1] = ptap->ris1;
51: cis[0] = ptap->cis0;
52: cis[1] = ptap->cis1;
53: n = ptap->ris1 ? 2 : 1;
54: reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
55: MatCreateSubMatrices(P,n,ris,cis,reuse,&ptap->lP);
57: MatISGetLocalMat(A,&lA);
58: MatISGetLocalMat(C,&lC);
59: if (ptap->ris1) { /* unsymmetric A mapping */
60: Mat lPt;
62: MatTranspose(ptap->lP[1],MAT_INITIAL_MATRIX,&lPt);
63: MatMatMatMult(lPt,lA,ptap->lP[0],reuse,ptap->fill,&lC);
64: if (matis->storel2l) {
65: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",(PetscObject)lPt);
66: }
67: MatDestroy(&lPt);
68: } else {
69: MatPtAP(lA,ptap->lP[0],reuse,ptap->fill,&lC);
70: if (matis->storel2l) {
71: PetscObjectCompose((PetscObject)C,"_MatIS_PtAP_l2l",(PetscObject)ptap->lP[0]);
72: }
73: }
74: if (reuse == MAT_INITIAL_MATRIX) {
75: MatISSetLocalMat(C,lC);
76: MatDestroy(&lC);
77: }
78: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
80: return(0);
81: }
83: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT,IS *cis)
84: {
85: Mat Po,Pd;
86: IS zd,zo;
87: const PetscInt *garray;
88: PetscInt *aux,i,bs;
89: PetscInt dc,stc,oc,ctd,cto;
90: PetscBool ismpiaij,ismpibaij,isseqaij,isseqbaij;
91: MPI_Comm comm;
97: PetscObjectGetComm((PetscObject)PT,&comm);
98: bs = 1;
99: PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIAIJ,&ismpiaij);
100: PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIBAIJ,&ismpibaij);
101: PetscObjectBaseTypeCompare((PetscObject)PT,MATSEQAIJ,&isseqaij);
102: PetscObjectTypeCompare((PetscObject)PT,MATSEQBAIJ,&isseqbaij);
103: if (isseqaij || isseqbaij) {
104: Pd = PT;
105: Po = NULL;
106: garray = NULL;
107: } else if (ismpiaij) {
108: MatMPIAIJGetSeqAIJ(PT,&Pd,&Po,&garray);
109: } else if (ismpibaij) {
110: MatMPIBAIJGetSeqBAIJ(PT,&Pd,&Po,&garray);
111: MatGetBlockSize(PT,&bs);
112: } else SETERRQ1(comm,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(PT))->type_name);
114: /* identify any null columns in Pd or Po */
115: /* We use a tolerance comparison since it may happen that, with geometric multigrid,
116: some of the columns are not really zero, but very close to */
117: zo = zd = NULL;
118: if (Po) {
119: MatFindNonzeroRowsOrCols_Basic(Po,PETSC_TRUE,PETSC_SMALL,&zo);
120: }
121: MatFindNonzeroRowsOrCols_Basic(Pd,PETSC_TRUE,PETSC_SMALL,&zd);
123: MatGetLocalSize(PT,NULL,&dc);
124: MatGetOwnershipRangeColumn(PT,&stc,NULL);
125: if (Po) { MatGetLocalSize(Po,NULL,&oc); }
126: else oc = 0;
127: PetscMalloc1((dc+oc)/bs,&aux);
128: if (zd) {
129: const PetscInt *idxs;
130: PetscInt nz;
132: /* this will throw an error if bs is not valid */
133: ISSetBlockSize(zd,bs);
134: ISGetLocalSize(zd,&nz);
135: ISGetIndices(zd,&idxs);
136: ctd = nz/bs;
137: for (i=0; i<ctd; i++) aux[i] = (idxs[bs*i]+stc)/bs;
138: ISRestoreIndices(zd,&idxs);
139: } else {
140: ctd = dc/bs;
141: for (i=0; i<ctd; i++) aux[i] = i+stc/bs;
142: }
143: if (zo) {
144: const PetscInt *idxs;
145: PetscInt nz;
147: /* this will throw an error if bs is not valid */
148: ISSetBlockSize(zo,bs);
149: ISGetLocalSize(zo,&nz);
150: ISGetIndices(zo,&idxs);
151: cto = nz/bs;
152: for (i=0; i<cto; i++) aux[i+ctd] = garray[idxs[bs*i]/bs];
153: ISRestoreIndices(zo,&idxs);
154: } else {
155: cto = oc/bs;
156: for (i=0; i<cto; i++) aux[i+ctd] = garray[i];
157: }
158: ISCreateBlock(comm,bs,ctd+cto,aux,PETSC_OWN_POINTER,cis);
159: ISDestroy(&zd);
160: ISDestroy(&zo);
161: return(0);
162: }
164: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A,Mat P,PetscReal fill,Mat C)
165: {
166: Mat PT,lA;
167: MatISPtAP ptap;
168: ISLocalToGlobalMapping Crl2g,Ccl2g,rl2g,cl2g;
169: PetscContainer c;
170: MatType lmtype;
171: const PetscInt *garray;
172: PetscInt ibs,N,dc;
173: MPI_Comm comm;
174: PetscErrorCode ierr;
177: PetscObjectGetComm((PetscObject)A,&comm);
178: MatSetType(C,MATIS);
179: MatISGetLocalMat(A,&lA);
180: MatGetType(lA,&lmtype);
181: MatISSetLocalMatType(C,lmtype);
182: MatGetSize(P,NULL,&N);
183: MatGetLocalSize(P,NULL,&dc);
184: MatSetSizes(C,dc,dc,N,N);
185: /* Not sure about this
186: MatGetBlockSizes(P,NULL,&ibs);
187: MatSetBlockSize(*C,ibs);
188: */
190: PetscNew(&ptap);
191: PetscContainerCreate(PETSC_COMM_SELF,&c);
192: PetscContainerSetPointer(c,ptap);
193: PetscContainerSetUserDestroy(c,MatISContainerDestroyPtAP_Private);
194: PetscObjectCompose((PetscObject)C,"_MatIS_PtAP",(PetscObject)c);
195: PetscContainerDestroy(&c);
196: ptap->fill = fill;
198: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
200: ISLocalToGlobalMappingGetBlockSize(cl2g,&ibs);
201: ISLocalToGlobalMappingGetSize(cl2g,&N);
202: ISLocalToGlobalMappingGetBlockIndices(cl2g,&garray);
203: ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris0);
204: ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&garray);
206: MatCreateSubMatrix(P,ptap->ris0,NULL,MAT_INITIAL_MATRIX,&PT);
207: MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis0);
208: ISLocalToGlobalMappingCreateIS(ptap->cis0,&Ccl2g);
209: MatDestroy(&PT);
211: Crl2g = NULL;
212: if (rl2g != cl2g) { /* unsymmetric A mapping */
213: PetscBool same,lsame = PETSC_FALSE;
214: PetscInt N1,ibs1;
216: ISLocalToGlobalMappingGetSize(rl2g,&N1);
217: ISLocalToGlobalMappingGetBlockSize(rl2g,&ibs1);
218: ISLocalToGlobalMappingGetBlockIndices(rl2g,&garray);
219: ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris1);
220: ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&garray);
221: if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
222: const PetscInt *i1,*i2;
224: ISBlockGetIndices(ptap->ris0,&i1);
225: ISBlockGetIndices(ptap->ris1,&i2);
226: PetscArraycmp(i1,i2,N,&lsame);
227: }
228: MPIU_Allreduce(&lsame,&same,1,MPIU_BOOL,MPI_LAND,comm);
229: if (same) {
230: ISDestroy(&ptap->ris1);
231: } else {
232: MatCreateSubMatrix(P,ptap->ris1,NULL,MAT_INITIAL_MATRIX,&PT);
233: MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis1);
234: ISLocalToGlobalMappingCreateIS(ptap->cis1,&Crl2g);
235: MatDestroy(&PT);
236: }
237: }
238: /* Not sure about this
239: if (!Crl2g) {
240: MatGetBlockSize(C,&ibs);
241: ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs);
242: }
243: */
244: MatSetLocalToGlobalMapping(C,Crl2g ? Crl2g : Ccl2g,Ccl2g);
245: ISLocalToGlobalMappingDestroy(&Crl2g);
246: ISLocalToGlobalMappingDestroy(&Ccl2g);
248: C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
249: return(0);
250: }
252: /* ----------------------------------------- */
253: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
254: {
256: Mat_Product *product = C->product;
257: Mat A=product->A,P=product->B;
258: PetscReal fill=product->fill;
261: MatPtAPSymbolic_IS_XAIJ(A,P,fill,C);
262: C->ops->productnumeric = MatProductNumeric_PtAP;
263: return(0);
264: }
266: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
267: {
269: C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
270: return(0);
271: }
273: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
274: {
276: Mat_Product *product = C->product;
279: if (product->type == MATPRODUCT_PtAP) {
280: MatProductSetFromOptions_IS_XAIJ_PtAP(C);
281: }
282: return(0);
283: }
285: /* ----------------------------------------- */
286: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
287: {
288: MatISLocalFields lf = (MatISLocalFields)ptr;
289: PetscInt i;
290: PetscErrorCode ierr;
293: for (i=0;i<lf->nr;i++) {
294: ISDestroy(&lf->rf[i]);
295: }
296: for (i=0;i<lf->nc;i++) {
297: ISDestroy(&lf->cf[i]);
298: }
299: PetscFree2(lf->rf,lf->cf);
300: PetscFree(lf);
301: return(0);
302: }
304: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
305: {
306: Mat B,lB;
310: if (reuse != MAT_REUSE_MATRIX) {
311: ISLocalToGlobalMapping rl2g,cl2g;
312: PetscInt bs;
313: IS is;
315: MatGetBlockSize(A,&bs);
316: ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is);
317: if (bs > 1) {
318: IS is2;
319: PetscInt i,*aux;
321: ISGetLocalSize(is,&i);
322: ISGetIndices(is,(const PetscInt**)&aux);
323: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2);
324: ISRestoreIndices(is,(const PetscInt**)&aux);
325: ISDestroy(&is);
326: is = is2;
327: }
328: ISSetIdentity(is);
329: ISLocalToGlobalMappingCreateIS(is,&rl2g);
330: ISDestroy(&is);
331: ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->n/bs,0,1,&is);
332: if (bs > 1) {
333: IS is2;
334: PetscInt i,*aux;
336: ISGetLocalSize(is,&i);
337: ISGetIndices(is,(const PetscInt**)&aux);
338: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2);
339: ISRestoreIndices(is,(const PetscInt**)&aux);
340: ISDestroy(&is);
341: is = is2;
342: }
343: ISSetIdentity(is);
344: ISLocalToGlobalMappingCreateIS(is,&cl2g);
345: ISDestroy(&is);
346: MatCreateIS(PetscObjectComm((PetscObject)A),bs,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,rl2g,cl2g,&B);
347: ISLocalToGlobalMappingDestroy(&rl2g);
348: ISLocalToGlobalMappingDestroy(&cl2g);
349: MatDuplicate(A,MAT_COPY_VALUES,&lB);
350: if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
351: } else {
352: B = *newmat;
353: PetscObjectReference((PetscObject)A);
354: lB = A;
355: }
356: MatISSetLocalMat(B,lB);
357: MatDestroy(&lB);
358: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
359: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
360: if (reuse == MAT_INPLACE_MATRIX) {
361: MatHeaderReplace(A,&B);
362: }
363: return(0);
364: }
366: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
367: {
368: Mat_IS *matis = (Mat_IS*)(A->data);
369: PetscScalar *aa;
370: const PetscInt *ii,*jj;
371: PetscInt i,n,m;
372: PetscInt *ecount,**eneighs;
373: PetscBool flg;
377: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
378: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
379: ISLocalToGlobalMappingGetNodeInfo(A->rmap->mapping,&n,&ecount,&eneighs);
380: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D != %D",m,n);
381: MatSeqAIJGetArray(matis->A,&aa);
382: for (i=0;i<n;i++) {
383: if (ecount[i] > 1) {
384: PetscInt j;
386: for (j=ii[i];j<ii[i+1];j++) {
387: PetscInt i2 = jj[j],p,p2;
388: PetscReal scal = 0.0;
390: for (p=0;p<ecount[i];p++) {
391: for (p2=0;p2<ecount[i2];p2++) {
392: if (eneighs[i][p] == eneighs[i2][p2]) { scal += 1.0; break; }
393: }
394: }
395: if (scal) aa[j] /= scal;
396: }
397: }
398: }
399: ISLocalToGlobalMappingRestoreNodeInfo(A->rmap->mapping,&n,&ecount,&eneighs);
400: MatSeqAIJRestoreArray(matis->A,&aa);
401: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
402: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
403: return(0);
404: }
406: typedef enum {MAT_IS_DISASSEMBLE_L2G_NATURAL,MAT_IS_DISASSEMBLE_L2G_MAT, MAT_IS_DISASSEMBLE_L2G_ND} MatISDisassemblel2gType;
408: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
409: {
410: Mat Ad,Ao;
411: IS is,ndmap,ndsub;
412: MPI_Comm comm;
413: const PetscInt *garray,*ndmapi;
414: PetscInt bs,i,cnt,nl,*ncount,*ndmapc;
415: PetscBool ismpiaij,ismpibaij;
416: const char *const MatISDisassemblel2gTypes[] = {"NATURAL","MAT","ND","MatISDisassemblel2gType","MAT_IS_DISASSEMBLE_L2G_",NULL};
417: MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
418: MatPartitioning part;
419: PetscSF sf;
420: PetscErrorCode ierr;
423: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatIS l2g disassembling options","Mat");
424: PetscOptionsEnum("-mat_is_disassemble_l2g_type","Type of local-to-global mapping to be used for disassembling","MatISDisassemblel2gType",MatISDisassemblel2gTypes,(PetscEnum)mode,(PetscEnum*)&mode,NULL);
425: PetscOptionsEnd();
426: if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
427: MatGetLocalToGlobalMapping(A,l2g,NULL);
428: return(0);
429: }
430: PetscObjectGetComm((PetscObject)A,&comm);
431: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij);
432: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij);
433: MatGetBlockSize(A,&bs);
434: switch (mode) {
435: case MAT_IS_DISASSEMBLE_L2G_ND:
436: MatPartitioningCreate(comm,&part);
437: MatPartitioningSetAdjacency(part,A);
438: PetscObjectSetOptionsPrefix((PetscObject)part,((PetscObject)A)->prefix);
439: MatPartitioningSetFromOptions(part);
440: MatPartitioningApplyND(part,&ndmap);
441: MatPartitioningDestroy(&part);
442: ISBuildTwoSided(ndmap,NULL,&ndsub);
443: MatMPIAIJSetUseScalableIncreaseOverlap(A,PETSC_TRUE);
444: MatIncreaseOverlap(A,1,&ndsub,1);
445: ISLocalToGlobalMappingCreateIS(ndsub,l2g);
447: /* it may happen that a separator node is not properly shared */
448: ISLocalToGlobalMappingGetNodeInfo(*l2g,&nl,&ncount,NULL);
449: PetscSFCreate(comm,&sf);
450: ISLocalToGlobalMappingGetIndices(*l2g,&garray);
451: PetscSFSetGraphLayout(sf,A->rmap,nl,NULL,PETSC_OWN_POINTER,garray);
452: ISLocalToGlobalMappingRestoreIndices(*l2g,&garray);
453: PetscCalloc1(A->rmap->n,&ndmapc);
454: PetscSFReduceBegin(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE);
455: PetscSFReduceEnd(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE);
456: ISLocalToGlobalMappingRestoreNodeInfo(*l2g,NULL,&ncount,NULL);
457: ISGetIndices(ndmap,&ndmapi);
458: for (i = 0, cnt = 0; i < A->rmap->n; i++)
459: if (ndmapi[i] < 0 && ndmapc[i] < 2)
460: cnt++;
462: MPIU_Allreduce(&cnt,&i,1,MPIU_INT,MPI_MAX,comm);
463: if (i) { /* we detected isolated separator nodes */
464: Mat A2,A3;
465: IS *workis,is2;
466: PetscScalar *vals;
467: PetscInt gcnt = i,*dnz,*onz,j,*lndmapi;
468: ISLocalToGlobalMapping ll2g;
469: PetscBool flg;
470: const PetscInt *ii,*jj;
472: /* communicate global id of separators */
473: MatPreallocateInitialize(comm,A->rmap->n,A->cmap->n,dnz,onz);
474: for (i = 0, cnt = 0; i < A->rmap->n; i++)
475: dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;
477: PetscMalloc1(nl,&lndmapi);
478: PetscSFBcastBegin(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE);
480: /* compute adjacency of isolated separators node */
481: PetscMalloc1(gcnt,&workis);
482: for (i = 0, cnt = 0; i < A->rmap->n; i++) {
483: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
484: ISCreateStride(comm,1,i+A->rmap->rstart,1,&workis[cnt++]);
485: }
486: }
487: for (i = cnt; i < gcnt; i++) {
488: ISCreateStride(comm,0,0,1,&workis[i]);
489: }
490: for (i = 0; i < gcnt; i++) {
491: PetscObjectSetName((PetscObject)workis[i],"ISOLATED");
492: ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
493: }
495: /* no communications since all the ISes correspond to locally owned rows */
496: MatIncreaseOverlap(A,gcnt,workis,1);
498: /* end communicate global id of separators */
499: PetscSFBcastEnd(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE);
501: /* communicate new layers : create a matrix and transpose it */
502: PetscArrayzero(dnz,A->rmap->n);
503: PetscArrayzero(onz,A->rmap->n);
504: for (i = 0, j = 0; i < A->rmap->n; i++) {
505: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
506: const PetscInt* idxs;
507: PetscInt s;
509: ISGetLocalSize(workis[j],&s);
510: ISGetIndices(workis[j],&idxs);
511: MatPreallocateSet(i+A->rmap->rstart,s,idxs,dnz,onz);
512: j++;
513: }
514: }
515: if (j != cnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %D != %D",j,cnt);
517: for (i = 0; i < gcnt; i++) {
518: PetscObjectSetName((PetscObject)workis[i],"EXTENDED");
519: ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
520: }
522: for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j,dnz[i]+onz[i]);
523: PetscMalloc1(j,&vals);
524: for (i = 0; i < j; i++) vals[i] = 1.0;
526: MatCreate(comm,&A2);
527: MatSetType(A2,MATMPIAIJ);
528: MatSetSizes(A2,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
529: MatMPIAIJSetPreallocation(A2,0,dnz,0,onz);
530: MatSetOption(A2,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
531: for (i = 0, j = 0; i < A2->rmap->n; i++) {
532: PetscInt row = i+A2->rmap->rstart,s = dnz[i] + onz[i];
533: const PetscInt* idxs;
535: if (s) {
536: ISGetIndices(workis[j],&idxs);
537: MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES);
538: ISRestoreIndices(workis[j],&idxs);
539: j++;
540: }
541: }
542: if (j != cnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %D != %D",j,cnt);
543: PetscFree(vals);
544: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
545: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
546: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
548: /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
549: for (i = 0, j = 0; i < nl; i++)
550: if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
551: ISCreateGeneral(comm,j,lndmapi,PETSC_USE_POINTER,&is);
552: MatMPIAIJGetLocalMatCondensed(A2,MAT_INITIAL_MATRIX,&is,NULL,&A3);
553: ISDestroy(&is);
554: MatDestroy(&A2);
556: /* extend local to global map to include connected isolated separators */
557: PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is);
558: if (!is) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing column map");
559: ISLocalToGlobalMappingCreateIS(is,&ll2g);
560: MatGetRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
561: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
562: ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is);
563: MatRestoreRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
564: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
565: ISLocalToGlobalMappingApplyIS(ll2g,is,&is2);
566: ISDestroy(&is);
567: ISLocalToGlobalMappingDestroy(&ll2g);
569: /* add new nodes to the local-to-global map */
570: ISLocalToGlobalMappingDestroy(l2g);
571: ISExpand(ndsub,is2,&is);
572: ISDestroy(&is2);
573: ISLocalToGlobalMappingCreateIS(is,l2g);
574: ISDestroy(&is);
576: MatDestroy(&A3);
577: PetscFree(lndmapi);
578: MatPreallocateFinalize(dnz,onz);
579: for (i = 0; i < gcnt; i++) {
580: ISDestroy(&workis[i]);
581: }
582: PetscFree(workis);
583: }
584: ISRestoreIndices(ndmap,&ndmapi);
585: PetscSFDestroy(&sf);
586: PetscFree(ndmapc);
587: ISDestroy(&ndmap);
588: ISDestroy(&ndsub);
589: ISLocalToGlobalMappingSetBlockSize(*l2g,bs);
590: ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view");
591: break;
592: case MAT_IS_DISASSEMBLE_L2G_NATURAL:
593: if (ismpiaij) {
594: MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
595: } else if (ismpibaij) {
596: MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
597: } else SETERRQ1(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
598: if (!garray) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present");
599: if (A->rmap->n) {
600: PetscInt dc,oc,stc,*aux;
602: MatGetLocalSize(A,NULL,&dc);
603: MatGetLocalSize(Ao,NULL,&oc);
604: MatGetOwnershipRangeColumn(A,&stc,NULL);
605: PetscMalloc1((dc+oc)/bs,&aux);
606: for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs;
607: for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
608: ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
609: } else {
610: ISCreateBlock(comm,1,0,NULL,PETSC_OWN_POINTER,&is);
611: }
612: ISLocalToGlobalMappingCreateIS(is,l2g);
613: ISDestroy(&is);
614: break;
615: default:
616: SETERRQ1(comm,PETSC_ERR_ARG_WRONG,"Unsupported l2g disassembling type %D",mode);
617: }
618: return(0);
619: }
621: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
622: {
623: Mat lA,Ad,Ao,B = NULL;
624: ISLocalToGlobalMapping rl2g,cl2g;
625: IS is;
626: MPI_Comm comm;
627: void *ptrs[2];
628: const char *names[2] = {"_convert_csr_aux","_convert_csr_data"};
629: const PetscInt *garray;
630: PetscScalar *dd,*od,*aa,*data;
631: const PetscInt *di,*dj,*oi,*oj;
632: const PetscInt *odi,*odj,*ooi,*ooj;
633: PetscInt *aux,*ii,*jj;
634: PetscInt bs,lc,dr,dc,oc,str,stc,nnz,i,jd,jo,cum;
635: PetscBool flg,ismpiaij,ismpibaij,was_inplace = PETSC_FALSE;
636: PetscMPIInt size;
637: PetscErrorCode ierr;
640: PetscObjectGetComm((PetscObject)A,&comm);
641: MPI_Comm_size(comm,&size);
642: if (size == 1) {
643: MatConvert_SeqXAIJ_IS(A,type,reuse,newmat);
644: return(0);
645: }
646: if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
647: MatMPIXAIJComputeLocalToGlobalMapping_Private(A,&rl2g);
648: MatCreate(comm,&B);
649: MatSetType(B,MATIS);
650: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
651: MatSetLocalToGlobalMapping(B,rl2g,rl2g);
652: MatGetBlockSize(A,&bs);
653: MatSetBlockSize(B,bs);
654: ISLocalToGlobalMappingDestroy(&rl2g);
655: if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
656: reuse = MAT_REUSE_MATRIX;
657: }
658: if (reuse == MAT_REUSE_MATRIX) {
659: Mat *newlA, lA;
660: IS rows, cols;
661: const PetscInt *ridx, *cidx;
662: PetscInt rbs, cbs, nr, nc;
664: if (!B) B = *newmat;
665: MatGetLocalToGlobalMapping(B,&rl2g,&cl2g);
666: ISLocalToGlobalMappingGetBlockIndices(rl2g,&ridx);
667: ISLocalToGlobalMappingGetBlockIndices(cl2g,&cidx);
668: ISLocalToGlobalMappingGetSize(rl2g,&nr);
669: ISLocalToGlobalMappingGetSize(cl2g,&nc);
670: ISLocalToGlobalMappingGetBlockSize(rl2g,&rbs);
671: ISLocalToGlobalMappingGetBlockSize(cl2g,&cbs);
672: ISCreateBlock(comm,rbs,nr/rbs,ridx,PETSC_USE_POINTER,&rows);
673: if (rl2g != cl2g) {
674: ISCreateBlock(comm,cbs,nc/cbs,cidx,PETSC_USE_POINTER,&cols);
675: } else {
676: PetscObjectReference((PetscObject)rows);
677: cols = rows;
678: }
679: MatISGetLocalMat(B,&lA);
680: MatCreateSubMatrices(A,1,&rows,&cols,MAT_INITIAL_MATRIX,&newlA);
681: MatConvert(newlA[0],MATSEQAIJ,MAT_INPLACE_MATRIX,&newlA[0]);
682: ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&ridx);
683: ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&cidx);
684: ISDestroy(&rows);
685: ISDestroy(&cols);
686: if (!lA->preallocated) { /* first time */
687: MatDuplicate(newlA[0],MAT_COPY_VALUES,&lA);
688: MatISSetLocalMat(B,lA);
689: PetscObjectDereference((PetscObject)lA);
690: }
691: MatCopy(newlA[0],lA,SAME_NONZERO_PATTERN);
692: MatDestroySubMatrices(1,&newlA);
693: MatISScaleDisassembling_Private(B);
694: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
695: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
696: if (was_inplace) { MatHeaderReplace(A,&B); }
697: else *newmat = B;
698: return(0);
699: }
700: /* rectangular case, just compress out the column space */
701: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij);
702: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij);
703: if (ismpiaij) {
704: bs = 1;
705: MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
706: } else if (ismpibaij) {
707: MatGetBlockSize(A,&bs);
708: MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
709: MatConvert(Ad,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ad);
710: MatConvert(Ao,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ao);
711: } else SETERRQ1(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
712: MatSeqAIJGetArray(Ad,&dd);
713: MatSeqAIJGetArray(Ao,&od);
714: if (!garray) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present");
716: /* access relevant information from MPIAIJ */
717: MatGetOwnershipRange(A,&str,NULL);
718: MatGetOwnershipRangeColumn(A,&stc,NULL);
719: MatGetLocalSize(A,&dr,&dc);
720: MatGetLocalSize(Ao,NULL,&oc);
721: MatGetRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&di,&dj,&flg);
722: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
723: MatGetRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&oi,&oj,&flg);
724: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
725: nnz = di[dr] + oi[dr];
726: /* store original pointers to be restored later */
727: odi = di; odj = dj; ooi = oi; ooj = oj;
729: /* generate l2g maps for rows and cols */
730: ISCreateStride(comm,dr/bs,str/bs,1,&is);
731: if (bs > 1) {
732: IS is2;
734: ISGetLocalSize(is,&i);
735: ISGetIndices(is,(const PetscInt**)&aux);
736: ISCreateBlock(comm,bs,i,aux,PETSC_COPY_VALUES,&is2);
737: ISRestoreIndices(is,(const PetscInt**)&aux);
738: ISDestroy(&is);
739: is = is2;
740: }
741: ISLocalToGlobalMappingCreateIS(is,&rl2g);
742: ISDestroy(&is);
743: if (dr) {
744: PetscMalloc1((dc+oc)/bs,&aux);
745: for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs;
746: for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
747: ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
748: lc = dc+oc;
749: } else {
750: ISCreateBlock(comm,bs,0,NULL,PETSC_OWN_POINTER,&is);
751: lc = 0;
752: }
753: ISLocalToGlobalMappingCreateIS(is,&cl2g);
754: ISDestroy(&is);
756: /* create MATIS object */
757: MatCreate(comm,&B);
758: MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
759: MatSetType(B,MATIS);
760: MatSetBlockSize(B,bs);
761: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
762: ISLocalToGlobalMappingDestroy(&rl2g);
763: ISLocalToGlobalMappingDestroy(&cl2g);
765: /* merge local matrices */
766: PetscMalloc1(nnz+dr+1,&aux);
767: PetscMalloc1(nnz,&data);
768: ii = aux;
769: jj = aux+dr+1;
770: aa = data;
771: *ii = *(di++) + *(oi++);
772: for (jd=0,jo=0,cum=0;*ii<nnz;cum++)
773: {
774: for (;jd<*di;jd++) { *jj++ = *dj++; *aa++ = *dd++; }
775: for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; }
776: *(++ii) = *(di++) + *(oi++);
777: }
778: for (;cum<dr;cum++) *(++ii) = nnz;
780: MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg);
781: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
782: MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg);
783: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
784: MatSeqAIJRestoreArray(Ad,&dd);
785: MatSeqAIJRestoreArray(Ao,&od);
787: ii = aux;
788: jj = aux+dr+1;
789: aa = data;
790: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);
792: /* create containers to destroy the data */
793: ptrs[0] = aux;
794: ptrs[1] = data;
795: for (i=0; i<2; i++) {
796: PetscContainer c;
798: PetscContainerCreate(PETSC_COMM_SELF,&c);
799: PetscContainerSetPointer(c,ptrs[i]);
800: PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);
801: PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
802: PetscContainerDestroy(&c);
803: }
804: if (ismpibaij) { /* destroy converted local matrices */
805: MatDestroy(&Ad);
806: MatDestroy(&Ao);
807: }
809: /* finalize matrix */
810: MatISSetLocalMat(B,lA);
811: MatDestroy(&lA);
812: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
813: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
814: if (reuse == MAT_INPLACE_MATRIX) {
815: MatHeaderReplace(A,&B);
816: } else *newmat = B;
817: return(0);
818: }
820: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
821: {
822: Mat **nest,*snest,**rnest,lA,B;
823: IS *iscol,*isrow,*islrow,*islcol;
824: ISLocalToGlobalMapping rl2g,cl2g;
825: MPI_Comm comm;
826: PetscInt *lr,*lc,*l2gidxs;
827: PetscInt i,j,nr,nc,rbs,cbs;
828: PetscBool convert,lreuse,*istrans;
829: PetscErrorCode ierr;
832: MatNestGetSubMats(A,&nr,&nc,&nest);
833: lreuse = PETSC_FALSE;
834: rnest = NULL;
835: if (reuse == MAT_REUSE_MATRIX) {
836: PetscBool ismatis,isnest;
838: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
839: if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
840: MatISGetLocalMat(*newmat,&lA);
841: PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
842: if (isnest) {
843: MatNestGetSubMats(lA,&i,&j,&rnest);
844: lreuse = (PetscBool)(i == nr && j == nc);
845: if (!lreuse) rnest = NULL;
846: }
847: }
848: PetscObjectGetComm((PetscObject)A,&comm);
849: PetscCalloc2(nr,&lr,nc,&lc);
850: PetscCalloc6(nr,&isrow,nc,&iscol,nr,&islrow,nc,&islcol,nr*nc,&snest,nr*nc,&istrans);
851: MatNestGetISs(A,isrow,iscol);
852: for (i=0;i<nr;i++) {
853: for (j=0;j<nc;j++) {
854: PetscBool ismatis;
855: PetscInt l1,l2,lb1,lb2,ij=i*nc+j;
857: /* Null matrix pointers are allowed in MATNEST */
858: if (!nest[i][j]) continue;
860: /* Nested matrices should be of type MATIS */
861: PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
862: if (istrans[ij]) {
863: Mat T,lT;
864: MatTransposeGetMat(nest[i][j],&T);
865: PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
866: 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);
867: MatISGetLocalMat(T,&lT);
868: MatCreateTranspose(lT,&snest[ij]);
869: } else {
870: PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
871: if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
872: MatISGetLocalMat(nest[i][j],&snest[ij]);
873: }
875: /* Check compatibility of local sizes */
876: MatGetSize(snest[ij],&l1,&l2);
877: MatGetBlockSizes(snest[ij],&lb1,&lb2);
878: if (!l1 || !l2) continue;
879: 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);
880: 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);
881: lr[i] = l1;
882: lc[j] = l2;
884: /* check compatibilty for local matrix reusage */
885: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
886: }
887: }
889: if (PetscDefined (USE_DEBUG)) {
890: /* Check compatibility of l2g maps for rows */
891: for (i=0;i<nr;i++) {
892: rl2g = NULL;
893: for (j=0;j<nc;j++) {
894: PetscInt n1,n2;
896: if (!nest[i][j]) continue;
897: if (istrans[i*nc+j]) {
898: Mat T;
900: MatTransposeGetMat(nest[i][j],&T);
901: MatGetLocalToGlobalMapping(T,NULL,&cl2g);
902: } else {
903: MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
904: }
905: ISLocalToGlobalMappingGetSize(cl2g,&n1);
906: if (!n1) continue;
907: if (!rl2g) {
908: rl2g = cl2g;
909: } else {
910: const PetscInt *idxs1,*idxs2;
911: PetscBool same;
913: ISLocalToGlobalMappingGetSize(rl2g,&n2);
914: 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);
915: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
916: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
917: PetscArraycmp(idxs1,idxs2,n1,&same);
918: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
919: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
920: 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);
921: }
922: }
923: }
924: /* Check compatibility of l2g maps for columns */
925: for (i=0;i<nc;i++) {
926: rl2g = NULL;
927: for (j=0;j<nr;j++) {
928: PetscInt n1,n2;
930: if (!nest[j][i]) continue;
931: if (istrans[j*nc+i]) {
932: Mat T;
934: MatTransposeGetMat(nest[j][i],&T);
935: MatGetLocalToGlobalMapping(T,&cl2g,NULL);
936: } else {
937: MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
938: }
939: ISLocalToGlobalMappingGetSize(cl2g,&n1);
940: if (!n1) continue;
941: if (!rl2g) {
942: rl2g = cl2g;
943: } else {
944: const PetscInt *idxs1,*idxs2;
945: PetscBool same;
947: ISLocalToGlobalMappingGetSize(rl2g,&n2);
948: 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);
949: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
950: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
951: PetscArraycmp(idxs1,idxs2,n1,&same);
952: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
953: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
954: 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);
955: }
956: }
957: }
958: }
960: B = NULL;
961: if (reuse != MAT_REUSE_MATRIX) {
962: PetscInt stl;
964: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
965: for (i=0,stl=0;i<nr;i++) stl += lr[i];
966: PetscMalloc1(stl,&l2gidxs);
967: for (i=0,stl=0;i<nr;i++) {
968: Mat usedmat;
969: Mat_IS *matis;
970: const PetscInt *idxs;
972: /* local IS for local NEST */
973: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
975: /* l2gmap */
976: j = 0;
977: usedmat = nest[i][j];
978: while (!usedmat && j < nc-1) usedmat = nest[i][++j];
979: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");
981: if (istrans[i*nc+j]) {
982: Mat T;
983: MatTransposeGetMat(usedmat,&T);
984: usedmat = T;
985: }
986: matis = (Mat_IS*)(usedmat->data);
987: ISGetIndices(isrow[i],&idxs);
988: if (istrans[i*nc+j]) {
989: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
990: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
991: } else {
992: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
993: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
994: }
995: ISRestoreIndices(isrow[i],&idxs);
996: stl += lr[i];
997: }
998: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);
1000: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
1001: for (i=0,stl=0;i<nc;i++) stl += lc[i];
1002: PetscMalloc1(stl,&l2gidxs);
1003: for (i=0,stl=0;i<nc;i++) {
1004: Mat usedmat;
1005: Mat_IS *matis;
1006: const PetscInt *idxs;
1008: /* local IS for local NEST */
1009: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
1011: /* l2gmap */
1012: j = 0;
1013: usedmat = nest[j][i];
1014: while (!usedmat && j < nr-1) usedmat = nest[++j][i];
1015: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat");
1016: if (istrans[j*nc+i]) {
1017: Mat T;
1018: MatTransposeGetMat(usedmat,&T);
1019: usedmat = T;
1020: }
1021: matis = (Mat_IS*)(usedmat->data);
1022: ISGetIndices(iscol[i],&idxs);
1023: if (istrans[j*nc+i]) {
1024: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1025: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1026: } else {
1027: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1028: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1029: }
1030: ISRestoreIndices(iscol[i],&idxs);
1031: stl += lc[i];
1032: }
1033: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);
1035: /* Create MATIS */
1036: MatCreate(comm,&B);
1037: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1038: MatGetBlockSizes(A,&rbs,&cbs);
1039: MatSetBlockSizes(B,rbs,cbs);
1040: MatSetType(B,MATIS);
1041: MatISSetLocalMatType(B,MATNEST);
1042: { /* hack : avoid setup of scatters */
1043: Mat_IS *matis = (Mat_IS*)(B->data);
1044: matis->islocalref = PETSC_TRUE;
1045: }
1046: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
1047: ISLocalToGlobalMappingDestroy(&rl2g);
1048: ISLocalToGlobalMappingDestroy(&cl2g);
1049: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1050: MatNestSetVecType(lA,VECNEST);
1051: for (i=0;i<nr*nc;i++) {
1052: if (istrans[i]) {
1053: MatDestroy(&snest[i]);
1054: }
1055: }
1056: MatISSetLocalMat(B,lA);
1057: MatDestroy(&lA);
1058: { /* hack : setup of scatters done here */
1059: Mat_IS *matis = (Mat_IS*)(B->data);
1061: matis->islocalref = PETSC_FALSE;
1062: MatISSetUpScatters_Private(B);
1063: }
1064: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1065: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1066: if (reuse == MAT_INPLACE_MATRIX) {
1067: MatHeaderReplace(A,&B);
1068: } else {
1069: *newmat = B;
1070: }
1071: } else {
1072: if (lreuse) {
1073: MatISGetLocalMat(*newmat,&lA);
1074: for (i=0;i<nr;i++) {
1075: for (j=0;j<nc;j++) {
1076: if (snest[i*nc+j]) {
1077: MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
1078: if (istrans[i*nc+j]) {
1079: MatDestroy(&snest[i*nc+j]);
1080: }
1081: }
1082: }
1083: }
1084: } else {
1085: PetscInt stl;
1086: for (i=0,stl=0;i<nr;i++) {
1087: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
1088: stl += lr[i];
1089: }
1090: for (i=0,stl=0;i<nc;i++) {
1091: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
1092: stl += lc[i];
1093: }
1094: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1095: for (i=0;i<nr*nc;i++) {
1096: if (istrans[i]) {
1097: MatDestroy(&snest[i]);
1098: }
1099: }
1100: MatISSetLocalMat(*newmat,lA);
1101: MatDestroy(&lA);
1102: }
1103: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1104: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1105: }
1107: /* Create local matrix in MATNEST format */
1108: convert = PETSC_FALSE;
1109: PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
1110: if (convert) {
1111: Mat M;
1112: MatISLocalFields lf;
1113: PetscContainer c;
1115: MatISGetLocalMat(*newmat,&lA);
1116: MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
1117: MatISSetLocalMat(*newmat,M);
1118: MatDestroy(&M);
1120: /* attach local fields to the matrix */
1121: PetscNew(&lf);
1122: PetscMalloc2(nr,&lf->rf,nc,&lf->cf);
1123: for (i=0;i<nr;i++) {
1124: PetscInt n,st;
1126: ISGetLocalSize(islrow[i],&n);
1127: ISStrideGetInfo(islrow[i],&st,NULL);
1128: ISCreateStride(comm,n,st,1,&lf->rf[i]);
1129: }
1130: for (i=0;i<nc;i++) {
1131: PetscInt n,st;
1133: ISGetLocalSize(islcol[i],&n);
1134: ISStrideGetInfo(islcol[i],&st,NULL);
1135: ISCreateStride(comm,n,st,1,&lf->cf[i]);
1136: }
1137: lf->nr = nr;
1138: lf->nc = nc;
1139: PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
1140: PetscContainerSetPointer(c,lf);
1141: PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
1142: PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
1143: PetscContainerDestroy(&c);
1144: }
1146: /* Free workspace */
1147: for (i=0;i<nr;i++) {
1148: ISDestroy(&islrow[i]);
1149: }
1150: for (i=0;i<nc;i++) {
1151: ISDestroy(&islcol[i]);
1152: }
1153: PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
1154: PetscFree2(lr,lc);
1155: return(0);
1156: }
1158: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1159: {
1160: Mat_IS *matis = (Mat_IS*)A->data;
1161: Vec ll,rr;
1162: const PetscScalar *Y,*X;
1163: PetscScalar *x,*y;
1164: PetscErrorCode ierr;
1167: if (l) {
1168: ll = matis->y;
1169: VecGetArrayRead(l,&Y);
1170: VecGetArray(ll,&y);
1171: PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE);
1172: } else {
1173: ll = NULL;
1174: }
1175: if (r) {
1176: rr = matis->x;
1177: VecGetArrayRead(r,&X);
1178: VecGetArray(rr,&x);
1179: PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE);
1180: } else {
1181: rr = NULL;
1182: }
1183: if (ll) {
1184: PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE);
1185: VecRestoreArrayRead(l,&Y);
1186: VecRestoreArray(ll,&y);
1187: }
1188: if (rr) {
1189: PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE);
1190: VecRestoreArrayRead(r,&X);
1191: VecRestoreArray(rr,&x);
1192: }
1193: MatDiagonalScale(matis->A,ll,rr);
1194: return(0);
1195: }
1197: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
1198: {
1199: Mat_IS *matis = (Mat_IS*)A->data;
1200: MatInfo info;
1201: PetscLogDouble isend[6],irecv[6];
1202: PetscInt bs;
1206: MatGetBlockSize(A,&bs);
1207: if (matis->A->ops->getinfo) {
1208: MatGetInfo(matis->A,MAT_LOCAL,&info);
1209: isend[0] = info.nz_used;
1210: isend[1] = info.nz_allocated;
1211: isend[2] = info.nz_unneeded;
1212: isend[3] = info.memory;
1213: isend[4] = info.mallocs;
1214: } else {
1215: isend[0] = 0.;
1216: isend[1] = 0.;
1217: isend[2] = 0.;
1218: isend[3] = 0.;
1219: isend[4] = 0.;
1220: }
1221: isend[5] = matis->A->num_ass;
1222: if (flag == MAT_LOCAL) {
1223: ginfo->nz_used = isend[0];
1224: ginfo->nz_allocated = isend[1];
1225: ginfo->nz_unneeded = isend[2];
1226: ginfo->memory = isend[3];
1227: ginfo->mallocs = isend[4];
1228: ginfo->assemblies = isend[5];
1229: } else if (flag == MAT_GLOBAL_MAX) {
1230: MPIU_Allreduce(isend,irecv,6,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)A));
1232: ginfo->nz_used = irecv[0];
1233: ginfo->nz_allocated = irecv[1];
1234: ginfo->nz_unneeded = irecv[2];
1235: ginfo->memory = irecv[3];
1236: ginfo->mallocs = irecv[4];
1237: ginfo->assemblies = irecv[5];
1238: } else if (flag == MAT_GLOBAL_SUM) {
1239: MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)A));
1241: ginfo->nz_used = irecv[0];
1242: ginfo->nz_allocated = irecv[1];
1243: ginfo->nz_unneeded = irecv[2];
1244: ginfo->memory = irecv[3];
1245: ginfo->mallocs = irecv[4];
1246: ginfo->assemblies = A->num_ass;
1247: }
1248: ginfo->block_size = bs;
1249: ginfo->fill_ratio_given = 0;
1250: ginfo->fill_ratio_needed = 0;
1251: ginfo->factor_mallocs = 0;
1252: return(0);
1253: }
1255: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
1256: {
1257: Mat C,lC,lA;
1258: PetscErrorCode ierr;
1261: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1262: ISLocalToGlobalMapping rl2g,cl2g;
1263: MatCreate(PetscObjectComm((PetscObject)A),&C);
1264: MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
1265: MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
1266: MatSetType(C,MATIS);
1267: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
1268: MatSetLocalToGlobalMapping(C,cl2g,rl2g);
1269: } else {
1270: C = *B;
1271: }
1273: /* perform local transposition */
1274: MatISGetLocalMat(A,&lA);
1275: MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
1276: MatISSetLocalMat(C,lC);
1277: MatDestroy(&lC);
1279: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1280: *B = C;
1281: } else {
1282: MatHeaderMerge(A,&C);
1283: }
1284: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1285: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1286: return(0);
1287: }
1289: PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
1290: {
1291: Mat_IS *is = (Mat_IS*)A->data;
1295: if (D) { /* MatShift_IS pass D = NULL */
1296: VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1297: VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1298: }
1299: VecPointwiseDivide(is->y,is->y,is->counter);
1300: MatDiagonalSet(is->A,is->y,insmode);
1301: return(0);
1302: }
1304: PetscErrorCode MatShift_IS(Mat A,PetscScalar a)
1305: {
1306: Mat_IS *is = (Mat_IS*)A->data;
1310: VecSet(is->y,a);
1311: MatDiagonalSet_IS(A,NULL,ADD_VALUES);
1312: return(0);
1313: }
1315: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1316: {
1318: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1321: if (PetscUnlikely(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);
1322: ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
1323: ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
1324: MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1325: return(0);
1326: }
1328: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1329: {
1331: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1334: if (PetscUnlikely(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);
1335: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
1336: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
1337: MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1338: return(0);
1339: }
1341: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
1342: {
1343: Mat locmat,newlocmat;
1344: Mat_IS *newmatis;
1345: const PetscInt *idxs;
1346: PetscInt i,m,n;
1347: PetscErrorCode ierr;
1350: if (scall == MAT_REUSE_MATRIX) {
1351: PetscBool ismatis;
1353: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
1354: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
1355: newmatis = (Mat_IS*)(*newmat)->data;
1356: if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
1357: if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
1358: }
1359: /* irow and icol may not have duplicate entries */
1360: if (PetscDefined(USE_DEBUG)) {
1361: Vec rtest,ltest;
1362: const PetscScalar *array;
1364: MatCreateVecs(mat,<est,&rtest);
1365: ISGetLocalSize(irow,&n);
1366: ISGetIndices(irow,&idxs);
1367: for (i=0;i<n;i++) {
1368: VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
1369: }
1370: VecAssemblyBegin(rtest);
1371: VecAssemblyEnd(rtest);
1372: VecGetLocalSize(rtest,&n);
1373: VecGetOwnershipRange(rtest,&m,NULL);
1374: VecGetArrayRead(rtest,&array);
1375: 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]));
1376: VecRestoreArrayRead(rtest,&array);
1377: ISRestoreIndices(irow,&idxs);
1378: ISGetLocalSize(icol,&n);
1379: ISGetIndices(icol,&idxs);
1380: for (i=0;i<n;i++) {
1381: VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
1382: }
1383: VecAssemblyBegin(ltest);
1384: VecAssemblyEnd(ltest);
1385: VecGetLocalSize(ltest,&n);
1386: VecGetOwnershipRange(ltest,&m,NULL);
1387: VecGetArrayRead(ltest,&array);
1388: 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]));
1389: VecRestoreArrayRead(ltest,&array);
1390: ISRestoreIndices(icol,&idxs);
1391: VecDestroy(&rtest);
1392: VecDestroy(<est);
1393: }
1394: if (scall == MAT_INITIAL_MATRIX) {
1395: Mat_IS *matis = (Mat_IS*)mat->data;
1396: ISLocalToGlobalMapping rl2g;
1397: IS is;
1398: PetscInt *lidxs,*lgidxs,*newgidxs;
1399: PetscInt ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
1400: PetscBool cong;
1401: MPI_Comm comm;
1403: PetscObjectGetComm((PetscObject)mat,&comm);
1404: MatGetBlockSizes(mat,&arbs,&acbs);
1405: ISGetBlockSize(irow,&irbs);
1406: ISGetBlockSize(icol,&icbs);
1407: rbs = arbs == irbs ? irbs : 1;
1408: cbs = acbs == icbs ? icbs : 1;
1409: ISGetLocalSize(irow,&m);
1410: ISGetLocalSize(icol,&n);
1411: MatCreate(comm,newmat);
1412: MatSetType(*newmat,MATIS);
1413: MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
1414: MatSetBlockSizes(*newmat,rbs,cbs);
1415: /* communicate irow to their owners in the layout */
1416: ISGetIndices(irow,&idxs);
1417: PetscLayoutMapLocal(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
1418: ISRestoreIndices(irow,&idxs);
1419: PetscArrayzero(matis->sf_rootdata,matis->sf->nroots);
1420: for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
1421: PetscFree(lidxs);
1422: PetscFree(lgidxs);
1423: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1424: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1425: for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
1426: PetscMalloc1(newloc,&newgidxs);
1427: PetscMalloc1(newloc,&lidxs);
1428: for (i=0,newloc=0;i<matis->sf->nleaves;i++)
1429: if (matis->sf_leafdata[i]) {
1430: lidxs[newloc] = i;
1431: newgidxs[newloc++] = matis->sf_leafdata[i]-1;
1432: }
1433: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1434: ISLocalToGlobalMappingCreateIS(is,&rl2g);
1435: ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);
1436: ISDestroy(&is);
1437: /* local is to extract local submatrix */
1438: newmatis = (Mat_IS*)(*newmat)->data;
1439: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
1440: MatHasCongruentLayouts(mat,&cong);
1441: if (cong && irow == icol && matis->csf == matis->sf) {
1442: MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
1443: PetscObjectReference((PetscObject)newmatis->getsub_ris);
1444: newmatis->getsub_cis = newmatis->getsub_ris;
1445: } else {
1446: ISLocalToGlobalMapping cl2g;
1448: /* communicate icol to their owners in the layout */
1449: ISGetIndices(icol,&idxs);
1450: PetscLayoutMapLocal(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
1451: ISRestoreIndices(icol,&idxs);
1452: PetscArrayzero(matis->csf_rootdata,matis->csf->nroots);
1453: for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
1454: PetscFree(lidxs);
1455: PetscFree(lgidxs);
1456: PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE);
1457: PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE);
1458: for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
1459: PetscMalloc1(newloc,&newgidxs);
1460: PetscMalloc1(newloc,&lidxs);
1461: for (i=0,newloc=0;i<matis->csf->nleaves;i++)
1462: if (matis->csf_leafdata[i]) {
1463: lidxs[newloc] = i;
1464: newgidxs[newloc++] = matis->csf_leafdata[i]-1;
1465: }
1466: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1467: ISLocalToGlobalMappingCreateIS(is,&cl2g);
1468: ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);
1469: ISDestroy(&is);
1470: /* local is to extract local submatrix */
1471: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
1472: MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
1473: ISLocalToGlobalMappingDestroy(&cl2g);
1474: }
1475: ISLocalToGlobalMappingDestroy(&rl2g);
1476: } else {
1477: MatISGetLocalMat(*newmat,&newlocmat);
1478: }
1479: MatISGetLocalMat(mat,&locmat);
1480: newmatis = (Mat_IS*)(*newmat)->data;
1481: MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
1482: if (scall == MAT_INITIAL_MATRIX) {
1483: MatISSetLocalMat(*newmat,newlocmat);
1484: MatDestroy(&newlocmat);
1485: }
1486: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1487: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1488: return(0);
1489: }
1491: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
1492: {
1493: Mat_IS *a = (Mat_IS*)A->data,*b;
1494: PetscBool ismatis;
1498: PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
1499: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
1500: b = (Mat_IS*)B->data;
1501: MatCopy(a->A,b->A,str);
1502: PetscObjectStateIncrease((PetscObject)B);
1503: return(0);
1504: }
1506: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d)
1507: {
1508: Vec v;
1509: const PetscScalar *array;
1510: PetscInt i,n;
1511: PetscErrorCode ierr;
1514: *missing = PETSC_FALSE;
1515: MatCreateVecs(A,NULL,&v);
1516: MatGetDiagonal(A,v);
1517: VecGetLocalSize(v,&n);
1518: VecGetArrayRead(v,&array);
1519: for (i=0;i<n;i++) if (array[i] == 0.) break;
1520: VecRestoreArrayRead(v,&array);
1521: VecDestroy(&v);
1522: if (i != n) *missing = PETSC_TRUE;
1523: if (d) {
1524: *d = -1;
1525: if (*missing) {
1526: PetscInt rstart;
1527: MatGetOwnershipRange(A,&rstart,NULL);
1528: *d = i+rstart;
1529: }
1530: }
1531: return(0);
1532: }
1534: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1535: {
1536: Mat_IS *matis = (Mat_IS*)(B->data);
1537: const PetscInt *gidxs;
1538: PetscInt nleaves;
1542: if (matis->sf) return(0);
1543: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
1544: ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
1545: ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
1546: PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1547: ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
1548: PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
1549: if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
1550: ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
1551: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
1552: ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
1553: PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1554: ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
1555: PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
1556: } else {
1557: matis->csf = matis->sf;
1558: matis->csf_leafdata = matis->sf_leafdata;
1559: matis->csf_rootdata = matis->sf_rootdata;
1560: }
1561: return(0);
1562: }
1564: /*@
1565: MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP.
1567: Collective
1569: Input Parameters:
1570: + A - the matrix
1571: - store - the boolean flag
1573: Level: advanced
1575: Notes:
1577: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
1578: @*/
1579: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1580: {
1587: PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
1588: return(0);
1589: }
1591: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1592: {
1593: Mat_IS *matis = (Mat_IS*)(A->data);
1597: matis->storel2l = store;
1598: if (!store) {
1599: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
1600: }
1601: return(0);
1602: }
1604: /*@
1605: MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1607: Collective
1609: Input Parameters:
1610: + A - the matrix
1611: - fix - the boolean flag
1613: Level: advanced
1615: Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process.
1617: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
1618: @*/
1619: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1620: {
1627: PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
1628: return(0);
1629: }
1631: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1632: {
1633: Mat_IS *matis = (Mat_IS*)(A->data);
1636: matis->locempty = fix;
1637: return(0);
1638: }
1640: /*@
1641: MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.
1643: Collective
1645: Input Parameters:
1646: + B - the matrix
1647: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1648: (same value is used for all local rows)
1649: . d_nnz - array containing the number of nonzeros in the various rows of the
1650: DIAGONAL portion of the local submatrix (possibly different for each row)
1651: or NULL, if d_nz is used to specify the nonzero structure.
1652: The size of this array is equal to the number of local rows, i.e 'm'.
1653: For matrices that will be factored, you must leave room for (and set)
1654: the diagonal entry even if it is zero.
1655: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1656: submatrix (same value is used for all local rows).
1657: - o_nnz - array containing the number of nonzeros in the various rows of the
1658: OFF-DIAGONAL portion of the local submatrix (possibly different for
1659: each row) or NULL, if o_nz is used to specify the nonzero
1660: structure. The size of this array is equal to the number
1661: of local rows, i.e 'm'.
1663: If the *_nnz parameter is given then the *_nz parameter is ignored
1665: Level: intermediate
1667: Notes:
1668: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
1669: from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
1670: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1672: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1673: @*/
1674: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1675: {
1681: PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1682: return(0);
1683: }
1685: /* this is used by DMDA */
1686: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1687: {
1688: Mat_IS *matis = (Mat_IS*)(B->data);
1689: PetscInt bs,i,nlocalcols;
1693: if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
1695: if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
1696: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
1698: if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz;
1699: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
1701: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1702: MatGetSize(matis->A,NULL,&nlocalcols);
1703: MatGetBlockSize(matis->A,&bs);
1704: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1706: for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1707: MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1708: #if defined(PETSC_HAVE_HYPRE)
1709: MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1710: #endif
1712: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
1713: MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1715: nlocalcols /= bs;
1716: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols - i);
1717: MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1719: /* for other matrix types */
1720: MatSetUp(matis->A);
1721: return(0);
1722: }
1724: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1725: {
1726: Mat_IS *matis = (Mat_IS*)(A->data);
1727: PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1728: const PetscInt *global_indices_r,*global_indices_c;
1729: PetscInt i,j,bs,rows,cols;
1730: PetscInt lrows,lcols;
1731: PetscInt local_rows,local_cols;
1732: PetscMPIInt size;
1733: PetscBool isdense,issbaij;
1734: PetscErrorCode ierr;
1737: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1738: MatGetSize(A,&rows,&cols);
1739: MatGetBlockSize(A,&bs);
1740: MatGetSize(matis->A,&local_rows,&local_cols);
1741: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1742: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1743: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
1744: if (A->rmap->mapping != A->cmap->mapping) {
1745: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
1746: } else {
1747: global_indices_c = global_indices_r;
1748: }
1750: if (issbaij) {
1751: MatGetRowUpperTriangular(matis->A);
1752: }
1753: /*
1754: An SF reduce is needed to sum up properly on shared rows.
1755: Note that generally preallocation is not exact, since it overestimates nonzeros
1756: */
1757: MatGetLocalSize(A,&lrows,&lcols);
1758: MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1759: /* All processes need to compute entire row ownership */
1760: PetscMalloc1(rows,&row_ownership);
1761: MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1762: for (i=0;i<size;i++) {
1763: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
1764: row_ownership[j] = i;
1765: }
1766: }
1767: MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);
1769: /*
1770: my_dnz and my_onz contains exact contribution to preallocation from each local mat
1771: then, they will be summed up properly. This way, preallocation is always sufficient
1772: */
1773: PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1774: /* preallocation as a MATAIJ */
1775: if (isdense) { /* special case for dense local matrices */
1776: for (i=0;i<local_rows;i++) {
1777: PetscInt owner = row_ownership[global_indices_r[i]];
1778: for (j=0;j<local_cols;j++) {
1779: PetscInt index_col = global_indices_c[j];
1780: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
1781: my_dnz[i] += 1;
1782: } else { /* offdiag block */
1783: my_onz[i] += 1;
1784: }
1785: }
1786: }
1787: } else if (matis->A->ops->getrowij) {
1788: const PetscInt *ii,*jj,*jptr;
1789: PetscBool done;
1790: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1791: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1792: jptr = jj;
1793: for (i=0;i<local_rows;i++) {
1794: PetscInt index_row = global_indices_r[i];
1795: for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1796: PetscInt owner = row_ownership[index_row];
1797: PetscInt index_col = global_indices_c[*jptr];
1798: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
1799: my_dnz[i] += 1;
1800: } else { /* offdiag block */
1801: my_onz[i] += 1;
1802: }
1803: /* same as before, interchanging rows and cols */
1804: if (issbaij && index_col != index_row) {
1805: owner = row_ownership[index_col];
1806: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) {
1807: my_dnz[*jptr] += 1;
1808: } else {
1809: my_onz[*jptr] += 1;
1810: }
1811: }
1812: }
1813: }
1814: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1815: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1816: } else { /* loop over rows and use MatGetRow */
1817: for (i=0;i<local_rows;i++) {
1818: const PetscInt *cols;
1819: PetscInt ncols,index_row = global_indices_r[i];
1820: MatGetRow(matis->A,i,&ncols,&cols,NULL);
1821: for (j=0;j<ncols;j++) {
1822: PetscInt owner = row_ownership[index_row];
1823: PetscInt index_col = global_indices_c[cols[j]];
1824: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
1825: my_dnz[i] += 1;
1826: } else { /* offdiag block */
1827: my_onz[i] += 1;
1828: }
1829: /* same as before, interchanging rows and cols */
1830: if (issbaij && index_col != index_row) {
1831: owner = row_ownership[index_col];
1832: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) {
1833: my_dnz[cols[j]] += 1;
1834: } else {
1835: my_onz[cols[j]] += 1;
1836: }
1837: }
1838: }
1839: MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1840: }
1841: }
1842: if (global_indices_c != global_indices_r) {
1843: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1844: }
1845: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1846: PetscFree(row_ownership);
1848: /* Reduce my_dnz and my_onz */
1849: if (maxreduce) {
1850: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1851: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1852: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1853: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1854: } else {
1855: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1856: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1857: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1858: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1859: }
1860: PetscFree2(my_dnz,my_onz);
1862: /* Resize preallocation if overestimated */
1863: for (i=0;i<lrows;i++) {
1864: dnz[i] = PetscMin(dnz[i],lcols);
1865: onz[i] = PetscMin(onz[i],cols-lcols);
1866: }
1868: /* Set preallocation */
1869: MatSeqAIJSetPreallocation(B,0,dnz);
1870: MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1871: for (i=0;i<lrows;i+=bs) {
1872: PetscInt b, d = dnz[i],o = onz[i];
1874: for (b=1;b<bs;b++) {
1875: d = PetscMax(d,dnz[i+b]);
1876: o = PetscMax(o,onz[i+b]);
1877: }
1878: dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
1879: onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
1880: }
1881: MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1882: MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1883: MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1884: MatPreallocateFinalize(dnz,onz);
1885: if (issbaij) {
1886: MatRestoreRowUpperTriangular(matis->A);
1887: }
1888: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1889: return(0);
1890: }
1892: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1893: {
1894: Mat_IS *matis = (Mat_IS*)(mat->data);
1895: Mat local_mat,MT;
1896: PetscInt rbs,cbs,rows,cols,lrows,lcols;
1897: PetscInt local_rows,local_cols;
1898: PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij;
1899: PetscMPIInt size;
1900: const PetscScalar *array;
1901: PetscErrorCode ierr;
1904: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
1905: if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1906: Mat B;
1907: IS irows = NULL,icols = NULL;
1908: PetscInt rbs,cbs;
1910: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
1911: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
1912: if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1913: IS rows,cols;
1914: const PetscInt *ridxs,*cidxs;
1915: PetscInt i,nw,*work;
1917: ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);
1918: ISLocalToGlobalMappingGetSize(mat->rmap->mapping,&nw);
1919: nw = nw/rbs;
1920: PetscCalloc1(nw,&work);
1921: for (i=0;i<nw;i++) work[ridxs[i]] += 1;
1922: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
1923: if (i == nw) {
1924: ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);
1925: ISSetPermutation(rows);
1926: ISInvertPermutation(rows,PETSC_DECIDE,&irows);
1927: ISDestroy(&rows);
1928: }
1929: ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);
1930: PetscFree(work);
1931: if (irows && mat->rmap->mapping != mat->cmap->mapping) {
1932: ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);
1933: ISLocalToGlobalMappingGetSize(mat->cmap->mapping,&nw);
1934: nw = nw/cbs;
1935: PetscCalloc1(nw,&work);
1936: for (i=0;i<nw;i++) work[cidxs[i]] += 1;
1937: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
1938: if (i == nw) {
1939: ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);
1940: ISSetPermutation(cols);
1941: ISInvertPermutation(cols,PETSC_DECIDE,&icols);
1942: ISDestroy(&cols);
1943: }
1944: ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);
1945: PetscFree(work);
1946: } else if (irows) {
1947: PetscObjectReference((PetscObject)irows);
1948: icols = irows;
1949: }
1950: } else {
1951: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);
1952: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);
1953: if (irows) { PetscObjectReference((PetscObject)irows); }
1954: if (icols) { PetscObjectReference((PetscObject)icols); }
1955: }
1956: if (!irows || !icols) {
1957: ISDestroy(&icols);
1958: ISDestroy(&irows);
1959: goto general_assembly;
1960: }
1961: MatConvert(matis->A,mtype,MAT_INITIAL_MATRIX,&B);
1962: if (reuse != MAT_INPLACE_MATRIX) {
1963: MatCreateSubMatrix(B,irows,icols,reuse,M);
1964: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);
1965: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);
1966: } else {
1967: Mat C;
1969: MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
1970: MatHeaderReplace(mat,&C);
1971: }
1972: MatDestroy(&B);
1973: ISDestroy(&icols);
1974: ISDestroy(&irows);
1975: return(0);
1976: }
1977: general_assembly:
1978: MatGetSize(mat,&rows,&cols);
1979: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
1980: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
1981: MatGetLocalSize(mat,&lrows,&lcols);
1982: MatGetSize(matis->A,&local_rows,&local_cols);
1983: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
1984: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
1985: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
1986: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
1987: if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
1988: if (PetscDefined (USE_DEBUG)) {
1989: PetscBool lb[4],bb[4];
1991: lb[0] = isseqdense;
1992: lb[1] = isseqaij;
1993: lb[2] = isseqbaij;
1994: lb[3] = isseqsbaij;
1995: MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
1996: if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
1997: }
1999: if (reuse != MAT_REUSE_MATRIX) {
2000: MatCreate(PetscObjectComm((PetscObject)mat),&MT);
2001: MatSetSizes(MT,lrows,lcols,rows,cols);
2002: MatSetType(MT,mtype);
2003: MatSetBlockSizes(MT,rbs,cbs);
2004: MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
2005: } else {
2006: PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;
2008: /* some checks */
2009: MT = *M;
2010: MatGetBlockSizes(MT,&mrbs,&mcbs);
2011: MatGetSize(MT,&mrows,&mcols);
2012: MatGetLocalSize(MT,&mlrows,&mlcols);
2013: if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
2014: if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
2015: if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
2016: if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
2017: if (mrbs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%d != %d)",rbs,mrbs);
2018: if (mcbs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%d != %d)",cbs,mcbs);
2019: MatZeroEntries(MT);
2020: }
2022: if (isseqsbaij || isseqbaij) {
2023: MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
2024: isseqaij = PETSC_TRUE;
2025: } else {
2026: PetscObjectReference((PetscObject)matis->A);
2027: local_mat = matis->A;
2028: }
2030: /* Set values */
2031: MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);
2032: if (isseqdense) { /* special case for dense local matrices */
2033: PetscInt i,*dummy;
2035: PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
2036: for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
2037: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
2038: MatDenseGetArrayRead(local_mat,&array);
2039: MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
2040: MatDenseRestoreArrayRead(local_mat,&array);
2041: PetscFree(dummy);
2042: } else if (isseqaij) {
2043: const PetscInt *blocks;
2044: PetscInt i,nvtxs,*xadj,*adjncy, nb;
2045: PetscBool done;
2046: PetscScalar *sarray;
2048: MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2049: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
2050: MatSeqAIJGetArray(local_mat,&sarray);
2051: MatGetVariableBlockSizes(local_mat,&nb,&blocks);
2052: if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2053: PetscInt sum;
2055: for (i=0,sum=0;i<nb;i++) sum += blocks[i];
2056: if (sum == nvtxs) {
2057: PetscInt r;
2059: for (i=0,r=0;i<nb;i++) {
2060: if (PetscUnlikelyDebug(blocks[i] != xadj[r+1] - xadj[r])) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %D: expected %D, got %D",i,blocks[i],xadj[r+1] - xadj[r]);
2061: MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],sarray+xadj[r],ADD_VALUES);
2062: r += blocks[i];
2063: }
2064: } else {
2065: for (i=0;i<nvtxs;i++) {
2066: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
2067: }
2068: }
2069: } else {
2070: for (i=0;i<nvtxs;i++) {
2071: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
2072: }
2073: }
2074: MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2075: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
2076: MatSeqAIJRestoreArray(local_mat,&sarray);
2077: } else { /* very basic values insertion for all other matrix types */
2078: PetscInt i;
2080: for (i=0;i<local_rows;i++) {
2081: PetscInt j;
2082: const PetscInt *local_indices_cols;
2084: MatGetRow(local_mat,i,&j,&local_indices_cols,&array);
2085: MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
2086: MatRestoreRow(local_mat,i,&j,&local_indices_cols,&array);
2087: }
2088: }
2089: MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
2090: MatDestroy(&local_mat);
2091: MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
2092: if (isseqdense) {
2093: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
2094: }
2095: if (reuse == MAT_INPLACE_MATRIX) {
2096: MatHeaderReplace(mat,&MT);
2097: } else if (reuse == MAT_INITIAL_MATRIX) {
2098: *M = MT;
2099: }
2100: return(0);
2101: }
2103: /*@
2104: MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format
2106: Input Parameters:
2107: + mat - the matrix (should be of type MATIS)
2108: - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
2110: Output Parameter:
2111: . newmat - the matrix in AIJ format
2113: Level: developer
2115: Notes:
2116: This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface.
2118: .seealso: MATIS, MatConvert()
2119: @*/
2120: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2121: {
2128: if (reuse == MAT_REUSE_MATRIX) {
2131: if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
2132: }
2133: PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
2134: return(0);
2135: }
2137: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
2138: {
2140: Mat_IS *matis = (Mat_IS*)(mat->data);
2141: PetscInt rbs,cbs,m,n,M,N;
2142: Mat B,localmat;
2145: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2146: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2147: MatGetSize(mat,&M,&N);
2148: MatGetLocalSize(mat,&m,&n);
2149: MatCreate(PetscObjectComm((PetscObject)mat),&B);
2150: MatSetSizes(B,m,n,M,N);
2151: MatSetBlockSize(B,rbs == cbs ? rbs : 1);
2152: MatSetType(B,MATIS);
2153: MatISSetLocalMatType(B,matis->lmattype);
2154: MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
2155: MatDuplicate(matis->A,op,&localmat);
2156: MatISSetLocalMat(B,localmat);
2157: MatDestroy(&localmat);
2158: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2159: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2160: *newmat = B;
2161: return(0);
2162: }
2164: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
2165: {
2167: Mat_IS *matis = (Mat_IS*)A->data;
2168: PetscBool local_sym;
2171: MatIsHermitian(matis->A,tol,&local_sym);
2172: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2173: return(0);
2174: }
2176: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
2177: {
2179: Mat_IS *matis = (Mat_IS*)A->data;
2180: PetscBool local_sym;
2183: MatIsSymmetric(matis->A,tol,&local_sym);
2184: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2185: return(0);
2186: }
2188: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg)
2189: {
2191: Mat_IS *matis = (Mat_IS*)A->data;
2192: PetscBool local_sym;
2195: if (A->rmap->mapping != A->cmap->mapping) {
2196: *flg = PETSC_FALSE;
2197: return(0);
2198: }
2199: MatIsStructurallySymmetric(matis->A,&local_sym);
2200: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2201: return(0);
2202: }
2204: static PetscErrorCode MatDestroy_IS(Mat A)
2205: {
2207: Mat_IS *b = (Mat_IS*)A->data;
2210: PetscFree(b->bdiag);
2211: PetscFree(b->lmattype);
2212: MatDestroy(&b->A);
2213: VecScatterDestroy(&b->cctx);
2214: VecScatterDestroy(&b->rctx);
2215: VecDestroy(&b->x);
2216: VecDestroy(&b->y);
2217: VecDestroy(&b->counter);
2218: ISDestroy(&b->getsub_ris);
2219: ISDestroy(&b->getsub_cis);
2220: if (b->sf != b->csf) {
2221: PetscSFDestroy(&b->csf);
2222: PetscFree2(b->csf_rootdata,b->csf_leafdata);
2223: } else b->csf = NULL;
2224: PetscSFDestroy(&b->sf);
2225: PetscFree2(b->sf_rootdata,b->sf_leafdata);
2226: PetscFree(A->data);
2227: PetscObjectChangeTypeName((PetscObject)A,NULL);
2228: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
2229: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
2230: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
2231: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
2232: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
2233: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
2234: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
2235: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
2236: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
2237: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
2238: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
2239: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
2240: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
2241: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
2242: return(0);
2243: }
2245: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
2246: {
2248: Mat_IS *is = (Mat_IS*)A->data;
2249: PetscScalar zero = 0.0;
2252: /* scatter the global vector x into the local work vector */
2253: VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2254: VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2256: /* multiply the local matrix */
2257: MatMult(is->A,is->x,is->y);
2259: /* scatter product back into global memory */
2260: VecSet(y,zero);
2261: VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2262: VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2263: return(0);
2264: }
2266: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2267: {
2268: Vec temp_vec;
2272: if (v3 != v2) {
2273: MatMult(A,v1,v3);
2274: VecAXPY(v3,1.0,v2);
2275: } else {
2276: VecDuplicate(v2,&temp_vec);
2277: MatMult(A,v1,temp_vec);
2278: VecAXPY(temp_vec,1.0,v2);
2279: VecCopy(temp_vec,v3);
2280: VecDestroy(&temp_vec);
2281: }
2282: return(0);
2283: }
2285: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
2286: {
2287: Mat_IS *is = (Mat_IS*)A->data;
2291: /* scatter the global vector x into the local work vector */
2292: VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2293: VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2295: /* multiply the local matrix */
2296: MatMultTranspose(is->A,is->y,is->x);
2298: /* scatter product back into global vector */
2299: VecSet(x,0);
2300: VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2301: VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2302: return(0);
2303: }
2305: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2306: {
2307: Vec temp_vec;
2311: if (v3 != v2) {
2312: MatMultTranspose(A,v1,v3);
2313: VecAXPY(v3,1.0,v2);
2314: } else {
2315: VecDuplicate(v2,&temp_vec);
2316: MatMultTranspose(A,v1,temp_vec);
2317: VecAXPY(temp_vec,1.0,v2);
2318: VecCopy(temp_vec,v3);
2319: VecDestroy(&temp_vec);
2320: }
2321: return(0);
2322: }
2324: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
2325: {
2326: Mat_IS *a = (Mat_IS*)A->data;
2328: PetscViewer sviewer;
2329: PetscBool isascii,view = PETSC_TRUE;
2332: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2333: if (isascii) {
2334: PetscViewerFormat format;
2336: PetscViewerGetFormat(viewer,&format);
2337: if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2338: }
2339: if (!view) return(0);
2340: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2341: MatView(a->A,sviewer);
2342: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2343: PetscViewerFlush(viewer);
2344: return(0);
2345: }
2347: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values)
2348: {
2349: Mat_IS *is = (Mat_IS*)mat->data;
2350: MPI_Datatype nodeType;
2351: const PetscScalar *lv;
2352: PetscInt bs;
2353: PetscErrorCode ierr;
2356: MatGetBlockSize(mat,&bs);
2357: MatSetBlockSize(is->A,bs);
2358: MatInvertBlockDiagonal(is->A,&lv);
2359: if (!is->bdiag) {
2360: PetscMalloc1(bs*mat->rmap->n,&is->bdiag);
2361: }
2362: MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType);
2363: MPI_Type_commit(&nodeType);
2364: PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
2365: PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
2366: MPI_Type_free(&nodeType);
2367: if (values) *values = is->bdiag;
2368: return(0);
2369: }
2371: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2372: {
2373: Vec cglobal,rglobal;
2374: IS from;
2375: Mat_IS *is = (Mat_IS*)A->data;
2376: PetscScalar sum;
2377: const PetscInt *garray;
2378: PetscInt nr,rbs,nc,cbs;
2379: PetscBool iscuda;
2383: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);
2384: ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);
2385: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);
2386: ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);
2387: VecDestroy(&is->x);
2388: VecDestroy(&is->y);
2389: VecDestroy(&is->counter);
2390: VecScatterDestroy(&is->rctx);
2391: VecScatterDestroy(&is->cctx);
2392: MatCreateVecs(is->A,&is->x,&is->y);
2393: VecBindToCPU(is->y,PETSC_TRUE);
2394: PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);
2395: if (iscuda) {
2396: PetscFree(A->defaultvectype);
2397: PetscStrallocpy(VECCUDA,&A->defaultvectype);
2398: }
2399: MatCreateVecs(A,&cglobal,&rglobal);
2400: VecBindToCPU(rglobal,PETSC_TRUE);
2401: ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);
2402: ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
2403: VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
2404: ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);
2405: ISDestroy(&from);
2406: if (A->rmap->mapping != A->cmap->mapping) {
2407: ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);
2408: ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
2409: VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
2410: ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);
2411: ISDestroy(&from);
2412: } else {
2413: PetscObjectReference((PetscObject)is->rctx);
2414: is->cctx = is->rctx;
2415: }
2416: VecDestroy(&cglobal);
2418: /* interface counter vector (local) */
2419: VecDuplicate(is->y,&is->counter);
2420: VecBindToCPU(is->counter,PETSC_TRUE);
2421: VecSet(is->y,1.);
2422: VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2423: VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2424: VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2425: VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2426: VecBindToCPU(is->y,PETSC_FALSE);
2427: VecBindToCPU(is->counter,PETSC_FALSE);
2429: /* special functions for block-diagonal matrices */
2430: VecSum(rglobal,&sum);
2431: if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && A->rmap->mapping == A->cmap->mapping) {
2432: A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2433: } else {
2434: A->ops->invertblockdiagonal = NULL;
2435: }
2436: VecDestroy(&rglobal);
2438: /* setup SF for general purpose shared indices based communications */
2439: MatISSetUpSF_IS(A);
2440: return(0);
2441: }
2443: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
2444: {
2446: PetscInt nr,rbs,nc,cbs;
2447: Mat_IS *is = (Mat_IS*)A->data;
2452: MatDestroy(&is->A);
2453: if (is->csf != is->sf) {
2454: PetscSFDestroy(&is->csf);
2455: PetscFree2(is->csf_rootdata,is->csf_leafdata);
2456: } else is->csf = NULL;
2457: PetscSFDestroy(&is->sf);
2458: PetscFree2(is->sf_rootdata,is->sf_leafdata);
2459: PetscFree(is->bdiag);
2461: /* Setup Layout and set local to global maps */
2462: PetscLayoutSetUp(A->rmap);
2463: PetscLayoutSetUp(A->cmap);
2464: ISLocalToGlobalMappingGetSize(rmapping,&nr);
2465: ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
2466: ISLocalToGlobalMappingGetSize(cmapping,&nc);
2467: ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
2468: /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */
2469: if (rmapping != cmapping && A->rmap->N == A->cmap->N) {
2470: PetscBool same,gsame;
2472: same = PETSC_FALSE;
2473: if (nr == nc && cbs == rbs) {
2474: const PetscInt *idxs1,*idxs2;
2476: ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);
2477: ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);
2478: PetscArraycmp(idxs1,idxs2,nr/rbs,&same);
2479: ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);
2480: ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);
2481: }
2482: MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2483: if (gsame) cmapping = rmapping;
2484: }
2485: PetscLayoutSetBlockSize(A->rmap,rbs);
2486: PetscLayoutSetBlockSize(A->cmap,cbs);
2487: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
2488: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
2490: /* Create the local matrix A */
2491: MatCreate(PETSC_COMM_SELF,&is->A);
2492: MatSetType(is->A,is->lmattype);
2493: MatSetSizes(is->A,nr,nc,nr,nc);
2494: MatSetBlockSizes(is->A,rbs,cbs);
2495: MatSetOptionsPrefix(is->A,"is_");
2496: MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);
2497: PetscLayoutSetUp(is->A->rmap);
2498: PetscLayoutSetUp(is->A->cmap);
2500: if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
2501: MatISSetUpScatters_Private(A);
2502: }
2503: MatSetUp(A);
2504: return(0);
2505: }
2507: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2508: {
2509: Mat_IS *is = (Mat_IS*)mat->data;
2511: PetscInt i,zm,zn;
2512: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2515: if (PetscDefined(USE_DEBUG)) {
2516: 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);
2517: /* count negative indices */
2518: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2519: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2520: }
2521: ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2522: ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2523: if (PetscDefined(USE_DEBUG)) {
2524: /* count negative indices (should be the same as before) */
2525: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2526: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2527: if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
2528: if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
2529: }
2530: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
2531: return(0);
2532: }
2534: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2535: {
2536: Mat_IS *is = (Mat_IS*)mat->data;
2538: PetscInt i,zm,zn;
2539: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2542: if (PetscDefined(USE_DEBUG)) {
2543: 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);
2544: /* count negative indices */
2545: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2546: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2547: }
2548: ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2549: ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2550: if (PetscDefined(USE_DEBUG)) {
2551: /* count negative indices (should be the same as before) */
2552: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2553: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2554: if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
2555: if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
2556: }
2557: MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
2558: return(0);
2559: }
2561: static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2562: {
2564: Mat_IS *is = (Mat_IS*)A->data;
2567: if (is->A->rmap->mapping) {
2568: MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
2569: } else {
2570: MatSetValues(is->A,m,rows,n,cols,values,addv);
2571: }
2572: return(0);
2573: }
2575: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2576: {
2578: Mat_IS *is = (Mat_IS*)A->data;
2581: if (is->A->rmap->mapping) {
2582: if (PetscDefined(USE_DEBUG)) {
2583: PetscInt ibs,bs;
2585: ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);
2586: MatGetBlockSize(is->A,&bs);
2587: if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs);
2588: }
2589: MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
2590: } else {
2591: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
2592: }
2593: return(0);
2594: }
2596: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns)
2597: {
2598: Mat_IS *is = (Mat_IS*)A->data;
2602: if (!n) {
2603: is->pure_neumann = PETSC_TRUE;
2604: } else {
2605: PetscInt i;
2606: is->pure_neumann = PETSC_FALSE;
2608: if (columns) {
2609: MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL);
2610: } else {
2611: MatZeroRows(is->A,n,rows,diag,NULL,NULL);
2612: }
2613: if (diag != 0.) {
2614: const PetscScalar *array;
2615: VecGetArrayRead(is->counter,&array);
2616: for (i=0; i<n; i++) {
2617: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
2618: }
2619: VecRestoreArrayRead(is->counter,&array);
2620: }
2621: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
2622: MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
2623: }
2624: return(0);
2625: }
2627: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
2628: {
2629: Mat_IS *matis = (Mat_IS*)A->data;
2630: PetscInt nr,nl,len,i;
2631: PetscInt *lrows;
2635: if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2636: PetscBool cong;
2638: PetscLayoutCompare(A->rmap,A->cmap,&cong);
2639: cong = (PetscBool)(cong && matis->sf == matis->csf);
2640: 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");
2641: 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");
2642: 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");
2643: }
2644: /* get locally owned rows */
2645: PetscLayoutMapLocal(A->rmap,n,rows,&len,&lrows,NULL);
2646: /* fix right hand side if needed */
2647: if (x && b) {
2648: const PetscScalar *xx;
2649: PetscScalar *bb;
2651: VecGetArrayRead(x, &xx);
2652: VecGetArray(b, &bb);
2653: for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
2654: VecRestoreArrayRead(x, &xx);
2655: VecRestoreArray(b, &bb);
2656: }
2657: /* get rows associated to the local matrices */
2658: MatGetSize(matis->A,&nl,NULL);
2659: PetscArrayzero(matis->sf_leafdata,nl);
2660: PetscArrayzero(matis->sf_rootdata,A->rmap->n);
2661: for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
2662: PetscFree(lrows);
2663: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
2664: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
2665: PetscMalloc1(nl,&lrows);
2666: for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
2667: MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
2668: PetscFree(lrows);
2669: return(0);
2670: }
2672: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2673: {
2677: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
2678: return(0);
2679: }
2681: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2682: {
2686: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
2687: return(0);
2688: }
2690: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
2691: {
2692: Mat_IS *is = (Mat_IS*)A->data;
2696: MatAssemblyBegin(is->A,type);
2697: return(0);
2698: }
2700: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
2701: {
2702: Mat_IS *is = (Mat_IS*)A->data;
2706: MatAssemblyEnd(is->A,type);
2707: /* fix for local empty rows/cols */
2708: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2709: Mat newlA;
2710: ISLocalToGlobalMapping rl2g,cl2g;
2711: IS nzr,nzc;
2712: PetscInt nr,nc,nnzr,nnzc;
2713: PetscBool lnewl2g,newl2g;
2715: MatGetSize(is->A,&nr,&nc);
2716: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
2717: if (!nzr) {
2718: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
2719: }
2720: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
2721: if (!nzc) {
2722: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
2723: }
2724: ISGetSize(nzr,&nnzr);
2725: ISGetSize(nzc,&nnzc);
2726: if (nnzr != nr || nnzc != nc) {
2727: ISLocalToGlobalMapping l2g;
2728: IS is1,is2;
2730: /* need new global l2g map */
2731: lnewl2g = PETSC_TRUE;
2732: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2734: /* extract valid submatrix */
2735: MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);
2737: /* attach local l2g maps for successive calls of MatSetValues on the local matrix */
2738: ISLocalToGlobalMappingCreateIS(nzr,&l2g);
2739: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);
2740: ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2741: ISLocalToGlobalMappingDestroy(&l2g);
2742: if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2743: const PetscInt *idxs1,*idxs2;
2744: PetscInt j,i,nl,*tidxs;
2746: ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);
2747: ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);
2748: PetscMalloc1(nl,&tidxs);
2749: ISGetIndices(is2,&idxs2);
2750: for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2751: if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr);
2752: ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);
2753: ISRestoreIndices(is2,&idxs2);
2754: ISDestroy(&is2);
2755: ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2756: }
2757: ISLocalToGlobalMappingCreateIS(is2,&rl2g);
2758: ISDestroy(&is1);
2759: ISDestroy(&is2);
2761: ISLocalToGlobalMappingCreateIS(nzc,&l2g);
2762: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);
2763: ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2764: ISLocalToGlobalMappingDestroy(&l2g);
2765: if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2766: const PetscInt *idxs1,*idxs2;
2767: PetscInt j,i,nl,*tidxs;
2769: ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);
2770: ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);
2771: PetscMalloc1(nl,&tidxs);
2772: ISGetIndices(is2,&idxs2);
2773: for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2774: if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc);
2775: ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);
2776: ISRestoreIndices(is2,&idxs2);
2777: ISDestroy(&is2);
2778: ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2779: }
2780: ISLocalToGlobalMappingCreateIS(is2,&cl2g);
2781: ISDestroy(&is1);
2782: ISDestroy(&is2);
2784: MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);
2786: ISLocalToGlobalMappingDestroy(&rl2g);
2787: ISLocalToGlobalMappingDestroy(&cl2g);
2788: } else { /* local matrix fully populated */
2789: lnewl2g = PETSC_FALSE;
2790: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2791: PetscObjectReference((PetscObject)is->A);
2792: newlA = is->A;
2793: }
2794: /* attach new global l2g map if needed */
2795: if (newl2g) {
2796: IS gnzr,gnzc;
2797: const PetscInt *grid,*gcid;
2799: ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);
2800: ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);
2801: ISGetIndices(gnzr,&grid);
2802: ISGetIndices(gnzc,&gcid);
2803: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);
2804: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);
2805: ISRestoreIndices(gnzr,&grid);
2806: ISRestoreIndices(gnzc,&gcid);
2807: ISDestroy(&gnzr);
2808: ISDestroy(&gnzc);
2809: MatSetLocalToGlobalMapping(A,rl2g,cl2g);
2810: ISLocalToGlobalMappingDestroy(&rl2g);
2811: ISLocalToGlobalMappingDestroy(&cl2g);
2812: }
2813: MatISSetLocalMat(A,newlA);
2814: MatDestroy(&newlA);
2815: ISDestroy(&nzr);
2816: ISDestroy(&nzc);
2817: is->locempty = PETSC_FALSE;
2818: }
2819: return(0);
2820: }
2822: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
2823: {
2824: Mat_IS *is = (Mat_IS*)mat->data;
2827: *local = is->A;
2828: return(0);
2829: }
2831: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
2832: {
2834: *local = NULL;
2835: return(0);
2836: }
2838: /*@
2839: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
2841: Input Parameter:
2842: . mat - the matrix
2844: Output Parameter:
2845: . local - the local matrix
2847: Level: advanced
2849: Notes:
2850: This can be called if you have precomputed the nonzero structure of the
2851: matrix and want to provide it to the inner matrix object to improve the performance
2852: of the MatSetValues() operation.
2854: Call MatISRestoreLocalMat() when finished with the local matrix.
2856: .seealso: MATIS
2857: @*/
2858: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
2859: {
2865: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
2866: return(0);
2867: }
2869: /*@
2870: MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()
2872: Input Parameter:
2873: . mat - the matrix
2875: Output Parameter:
2876: . local - the local matrix
2878: Level: advanced
2880: .seealso: MATIS
2881: @*/
2882: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
2883: {
2889: PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
2890: return(0);
2891: }
2893: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
2894: {
2895: Mat_IS *is = (Mat_IS*)mat->data;
2899: if (is->A) {
2900: MatSetType(is->A,mtype);
2901: }
2902: PetscFree(is->lmattype);
2903: PetscStrallocpy(mtype,&is->lmattype);
2904: return(0);
2905: }
2907: /*@
2908: MatISSetLocalMatType - Specifies the type of local matrix
2910: Input Parameters:
2911: + mat - the matrix
2912: - mtype - the local matrix type
2914: Output Parameter:
2916: Level: advanced
2918: .seealso: MATIS, MatSetType(), MatType
2919: @*/
2920: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
2921: {
2926: PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
2927: return(0);
2928: }
2930: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
2931: {
2932: Mat_IS *is = (Mat_IS*)mat->data;
2933: PetscInt nrows,ncols,orows,ocols;
2935: MatType mtype,otype;
2936: PetscBool sametype = PETSC_TRUE;
2939: if (is->A) {
2940: MatGetSize(is->A,&orows,&ocols);
2941: MatGetSize(local,&nrows,&ncols);
2942: 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);
2943: MatGetType(local,&mtype);
2944: MatGetType(is->A,&otype);
2945: PetscStrcmp(mtype,otype,&sametype);
2946: }
2947: PetscObjectReference((PetscObject)local);
2948: MatDestroy(&is->A);
2949: is->A = local;
2950: MatGetType(is->A,&mtype);
2951: MatISSetLocalMatType(mat,mtype);
2952: if (!sametype && !is->islocalref) {
2953: MatISSetUpScatters_Private(mat);
2954: }
2955: return(0);
2956: }
2958: /*@
2959: MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.
2961: Collective on Mat
2963: Input Parameters:
2964: + mat - the matrix
2965: - local - the local matrix
2967: Output Parameter:
2969: Level: advanced
2971: Notes:
2972: This can be called if you have precomputed the local matrix and
2973: want to provide it to the matrix object MATIS.
2975: .seealso: MATIS
2976: @*/
2977: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
2978: {
2984: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
2985: return(0);
2986: }
2988: static PetscErrorCode MatZeroEntries_IS(Mat A)
2989: {
2990: Mat_IS *a = (Mat_IS*)A->data;
2994: MatZeroEntries(a->A);
2995: return(0);
2996: }
2998: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
2999: {
3000: Mat_IS *is = (Mat_IS*)A->data;
3004: MatScale(is->A,a);
3005: return(0);
3006: }
3008: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3009: {
3010: Mat_IS *is = (Mat_IS*)A->data;
3014: /* get diagonal of the local matrix */
3015: MatGetDiagonal(is->A,is->y);
3017: /* scatter diagonal back into global vector */
3018: VecSet(v,0);
3019: VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3020: VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3021: return(0);
3022: }
3024: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
3025: {
3026: Mat_IS *a = (Mat_IS*)A->data;
3030: MatSetOption(a->A,op,flg);
3031: return(0);
3032: }
3034: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
3035: {
3036: Mat_IS *y = (Mat_IS*)Y->data;
3037: Mat_IS *x;
3041: if (PetscDefined(USE_DEBUG)) {
3042: PetscBool ismatis;
3043: PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
3044: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3045: }
3046: x = (Mat_IS*)X->data;
3047: MatAXPY(y->A,a,x->A,str);
3048: return(0);
3049: }
3051: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
3052: {
3053: Mat lA;
3054: Mat_IS *matis;
3055: ISLocalToGlobalMapping rl2g,cl2g;
3056: IS is;
3057: const PetscInt *rg,*rl;
3058: PetscInt nrg;
3059: PetscInt N,M,nrl,i,*idxs;
3060: PetscErrorCode ierr;
3063: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
3064: ISGetLocalSize(row,&nrl);
3065: ISGetIndices(row,&rl);
3066: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
3067: if (PetscDefined(USE_DEBUG)) {
3068: for (i=0; i<nrl; i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater then maximum possible %D",i,rl[i],nrg);
3069: }
3070: PetscMalloc1(nrg,&idxs);
3071: /* map from [0,nrl) to row */
3072: for (i=0;i<nrl;i++) idxs[i] = rl[i];
3073: for (i=nrl;i<nrg;i++) idxs[i] = -1;
3074: ISRestoreIndices(row,&rl);
3075: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
3076: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
3077: ISLocalToGlobalMappingCreateIS(is,&rl2g);
3078: ISDestroy(&is);
3079: /* compute new l2g map for columns */
3080: if (col != row || A->rmap->mapping != A->cmap->mapping) {
3081: const PetscInt *cg,*cl;
3082: PetscInt ncg;
3083: PetscInt ncl;
3085: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
3086: ISGetLocalSize(col,&ncl);
3087: ISGetIndices(col,&cl);
3088: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
3089: if (PetscDefined(USE_DEBUG)) {
3090: for (i=0; i<ncl; i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater then maximum possible %D",i,cl[i],ncg);
3091: }
3092: PetscMalloc1(ncg,&idxs);
3093: /* map from [0,ncl) to col */
3094: for (i=0;i<ncl;i++) idxs[i] = cl[i];
3095: for (i=ncl;i<ncg;i++) idxs[i] = -1;
3096: ISRestoreIndices(col,&cl);
3097: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
3098: ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
3099: ISLocalToGlobalMappingCreateIS(is,&cl2g);
3100: ISDestroy(&is);
3101: } else {
3102: PetscObjectReference((PetscObject)rl2g);
3103: cl2g = rl2g;
3104: }
3105: /* create the MATIS submatrix */
3106: MatGetSize(A,&M,&N);
3107: MatCreate(PetscObjectComm((PetscObject)A),submat);
3108: MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
3109: MatSetType(*submat,MATIS);
3110: matis = (Mat_IS*)((*submat)->data);
3111: matis->islocalref = PETSC_TRUE;
3112: MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
3113: MatISGetLocalMat(A,&lA);
3114: MatISSetLocalMat(*submat,lA);
3115: MatSetUp(*submat);
3116: MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
3117: MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
3118: ISLocalToGlobalMappingDestroy(&rl2g);
3119: ISLocalToGlobalMappingDestroy(&cl2g);
3120: /* remove unsupported ops */
3121: PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
3122: (*submat)->ops->destroy = MatDestroy_IS;
3123: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3124: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3125: (*submat)->ops->assemblybegin = MatAssemblyBegin_IS;
3126: (*submat)->ops->assemblyend = MatAssemblyEnd_IS;
3127: return(0);
3128: }
3130: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
3131: {
3132: Mat_IS *a = (Mat_IS*)A->data;
3133: char type[256];
3134: PetscBool flg;
3138: PetscOptionsHead(PetscOptionsObject,"MATIS options");
3139: PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
3140: PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
3141: PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
3142: if (flg) {
3143: MatISSetLocalMatType(A,type);
3144: }
3145: if (a->A) {
3146: MatSetFromOptions(a->A);
3147: }
3148: PetscOptionsTail();
3149: return(0);
3150: }
3152: /*@
3153: MatCreateIS - Creates a "process" unassembled matrix, assembled on each
3154: process but not across processes.
3156: Input Parameters:
3157: + comm - MPI communicator that will share the matrix
3158: . bs - block size of the matrix
3159: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3160: . rmap - local to global map for rows
3161: - cmap - local to global map for cols
3163: Output Parameter:
3164: . A - the resulting matrix
3166: Level: advanced
3168: Notes:
3169: See MATIS for more details.
3170: m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors
3171: used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
3172: If either rmap or cmap are NULL, then the matrix is assumed to be square.
3174: .seealso: MATIS, MatSetLocalToGlobalMapping()
3175: @*/
3176: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
3177: {
3181: if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
3182: MatCreate(comm,A);
3183: MatSetSizes(*A,m,n,M,N);
3184: if (bs > 0) {
3185: MatSetBlockSize(*A,bs);
3186: }
3187: MatSetType(*A,MATIS);
3188: if (rmap && cmap) {
3189: MatSetLocalToGlobalMapping(*A,rmap,cmap);
3190: } else if (!rmap) {
3191: MatSetLocalToGlobalMapping(*A,cmap,cmap);
3192: } else {
3193: MatSetLocalToGlobalMapping(*A,rmap,rmap);
3194: }
3195: return(0);
3196: }
3198: /*MC
3199: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP).
3200: This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector
3201: product is handled "implicitly".
3203: Options Database Keys:
3204: + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
3205: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3206: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP().
3208: Notes:
3209: Options prefix for the inner matrix are given by -is_mat_xxx
3211: You must call MatSetLocalToGlobalMapping() before using this matrix type.
3213: You can do matrix preallocation on the local matrix after you obtain it with
3214: MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation()
3216: Level: advanced
3218: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP
3220: M*/
3222: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3223: {
3225: Mat_IS *b;
3228: PetscNewLog(A,&b);
3229: PetscStrallocpy(MATAIJ,&b->lmattype);
3230: A->data = (void*)b;
3232: /* matrix ops */
3233: PetscMemzero(A->ops,sizeof(struct _MatOps));
3234: A->ops->mult = MatMult_IS;
3235: A->ops->multadd = MatMultAdd_IS;
3236: A->ops->multtranspose = MatMultTranspose_IS;
3237: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3238: A->ops->destroy = MatDestroy_IS;
3239: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3240: A->ops->setvalues = MatSetValues_IS;
3241: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3242: A->ops->setvalueslocal = MatSetValuesLocal_IS;
3243: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3244: A->ops->zerorows = MatZeroRows_IS;
3245: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3246: A->ops->assemblybegin = MatAssemblyBegin_IS;
3247: A->ops->assemblyend = MatAssemblyEnd_IS;
3248: A->ops->view = MatView_IS;
3249: A->ops->zeroentries = MatZeroEntries_IS;
3250: A->ops->scale = MatScale_IS;
3251: A->ops->getdiagonal = MatGetDiagonal_IS;
3252: A->ops->setoption = MatSetOption_IS;
3253: A->ops->ishermitian = MatIsHermitian_IS;
3254: A->ops->issymmetric = MatIsSymmetric_IS;
3255: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3256: A->ops->duplicate = MatDuplicate_IS;
3257: A->ops->missingdiagonal = MatMissingDiagonal_IS;
3258: A->ops->copy = MatCopy_IS;
3259: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3260: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3261: A->ops->axpy = MatAXPY_IS;
3262: A->ops->diagonalset = MatDiagonalSet_IS;
3263: A->ops->shift = MatShift_IS;
3264: A->ops->transpose = MatTranspose_IS;
3265: A->ops->getinfo = MatGetInfo_IS;
3266: A->ops->diagonalscale = MatDiagonalScale_IS;
3267: A->ops->setfromoptions = MatSetFromOptions_IS;
3269: /* special MATIS functions */
3270: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
3271: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
3272: PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
3273: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
3274: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
3275: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
3276: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
3277: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
3278: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
3279: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
3280: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
3281: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
3282: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
3283: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
3284: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
3285: PetscObjectChangeTypeName((PetscObject)A,MATIS);
3286: return(0);
3287: }