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>
13: #include <petsc/private/hashseti.h>
15: #define MATIS_MAX_ENTRIES_INSERTION 2048
16: static PetscErrorCode MatSetValuesLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode);
17: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode);
18: static PetscErrorCode MatISSetUpScatters_Private(Mat);
20: static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr)
21: {
22: MatISPtAP ptap = (MatISPtAP)ptr;
24: MatDestroySubMatrices(ptap->ris1 ? 2 : 1,&ptap->lP);
25: ISDestroy(&ptap->cis0);
26: ISDestroy(&ptap->cis1);
27: ISDestroy(&ptap->ris0);
28: ISDestroy(&ptap->ris1);
29: PetscFree(ptap);
30: return 0;
31: }
33: static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
34: {
35: MatISPtAP ptap;
36: Mat_IS *matis = (Mat_IS*)A->data;
37: Mat lA,lC;
38: MatReuse reuse;
39: IS ris[2],cis[2];
40: PetscContainer c;
41: PetscInt n;
43: PetscObjectQuery((PetscObject)C,"_MatIS_PtAP",(PetscObject*)&c);
45: PetscContainerGetPointer(c,(void**)&ptap);
46: ris[0] = ptap->ris0;
47: ris[1] = ptap->ris1;
48: cis[0] = ptap->cis0;
49: cis[1] = ptap->cis1;
50: n = ptap->ris1 ? 2 : 1;
51: reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
52: MatCreateSubMatrices(P,n,ris,cis,reuse,&ptap->lP);
54: MatISGetLocalMat(A,&lA);
55: MatISGetLocalMat(C,&lC);
56: if (ptap->ris1) { /* unsymmetric A mapping */
57: Mat lPt;
59: MatTranspose(ptap->lP[1],MAT_INITIAL_MATRIX,&lPt);
60: MatMatMatMult(lPt,lA,ptap->lP[0],reuse,ptap->fill,&lC);
61: if (matis->storel2l) {
62: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",(PetscObject)lPt);
63: }
64: MatDestroy(&lPt);
65: } else {
66: MatPtAP(lA,ptap->lP[0],reuse,ptap->fill,&lC);
67: if (matis->storel2l) {
68: PetscObjectCompose((PetscObject)C,"_MatIS_PtAP_l2l",(PetscObject)ptap->lP[0]);
69: }
70: }
71: if (reuse == MAT_INITIAL_MATRIX) {
72: MatISSetLocalMat(C,lC);
73: MatDestroy(&lC);
74: }
75: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
77: return 0;
78: }
80: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT,IS *cis)
81: {
82: Mat Po,Pd;
83: IS zd,zo;
84: const PetscInt *garray;
85: PetscInt *aux,i,bs;
86: PetscInt dc,stc,oc,ctd,cto;
87: PetscBool ismpiaij,ismpibaij,isseqaij,isseqbaij;
88: MPI_Comm comm;
92: PetscObjectGetComm((PetscObject)PT,&comm);
93: bs = 1;
94: PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIAIJ,&ismpiaij);
95: PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIBAIJ,&ismpibaij);
96: PetscObjectBaseTypeCompare((PetscObject)PT,MATSEQAIJ,&isseqaij);
97: PetscObjectTypeCompare((PetscObject)PT,MATSEQBAIJ,&isseqbaij);
98: if (isseqaij || isseqbaij) {
99: Pd = PT;
100: Po = NULL;
101: garray = NULL;
102: } else if (ismpiaij) {
103: MatMPIAIJGetSeqAIJ(PT,&Pd,&Po,&garray);
104: } else if (ismpibaij) {
105: MatMPIBAIJGetSeqBAIJ(PT,&Pd,&Po,&garray);
106: MatGetBlockSize(PT,&bs);
107: } else SETERRQ(comm,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(PT))->type_name);
109: /* identify any null columns in Pd or Po */
110: /* We use a tolerance comparison since it may happen that, with geometric multigrid,
111: some of the columns are not really zero, but very close to */
112: zo = zd = NULL;
113: if (Po) {
114: MatFindNonzeroRowsOrCols_Basic(Po,PETSC_TRUE,PETSC_SMALL,&zo);
115: }
116: MatFindNonzeroRowsOrCols_Basic(Pd,PETSC_TRUE,PETSC_SMALL,&zd);
118: MatGetLocalSize(PT,NULL,&dc);
119: MatGetOwnershipRangeColumn(PT,&stc,NULL);
120: if (Po) MatGetLocalSize(Po,NULL,&oc);
121: else oc = 0;
122: PetscMalloc1((dc+oc)/bs,&aux);
123: if (zd) {
124: const PetscInt *idxs;
125: PetscInt nz;
127: /* this will throw an error if bs is not valid */
128: ISSetBlockSize(zd,bs);
129: ISGetLocalSize(zd,&nz);
130: ISGetIndices(zd,&idxs);
131: ctd = nz/bs;
132: for (i=0; i<ctd; i++) aux[i] = (idxs[bs*i]+stc)/bs;
133: ISRestoreIndices(zd,&idxs);
134: } else {
135: ctd = dc/bs;
136: for (i=0; i<ctd; i++) aux[i] = i+stc/bs;
137: }
138: if (zo) {
139: const PetscInt *idxs;
140: PetscInt nz;
142: /* this will throw an error if bs is not valid */
143: ISSetBlockSize(zo,bs);
144: ISGetLocalSize(zo,&nz);
145: ISGetIndices(zo,&idxs);
146: cto = nz/bs;
147: for (i=0; i<cto; i++) aux[i+ctd] = garray[idxs[bs*i]/bs];
148: ISRestoreIndices(zo,&idxs);
149: } else {
150: cto = oc/bs;
151: for (i=0; i<cto; i++) aux[i+ctd] = garray[i];
152: }
153: ISCreateBlock(comm,bs,ctd+cto,aux,PETSC_OWN_POINTER,cis);
154: ISDestroy(&zd);
155: ISDestroy(&zo);
156: return 0;
157: }
159: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A,Mat P,PetscReal fill,Mat C)
160: {
161: Mat PT,lA;
162: MatISPtAP ptap;
163: ISLocalToGlobalMapping Crl2g,Ccl2g,rl2g,cl2g;
164: PetscContainer c;
165: MatType lmtype;
166: const PetscInt *garray;
167: PetscInt ibs,N,dc;
168: MPI_Comm comm;
170: PetscObjectGetComm((PetscObject)A,&comm);
171: MatSetType(C,MATIS);
172: MatISGetLocalMat(A,&lA);
173: MatGetType(lA,&lmtype);
174: MatISSetLocalMatType(C,lmtype);
175: MatGetSize(P,NULL,&N);
176: MatGetLocalSize(P,NULL,&dc);
177: MatSetSizes(C,dc,dc,N,N);
178: /* Not sure about this
179: MatGetBlockSizes(P,NULL,&ibs);
180: MatSetBlockSize(*C,ibs);
181: */
183: PetscNew(&ptap);
184: PetscContainerCreate(PETSC_COMM_SELF,&c);
185: PetscContainerSetPointer(c,ptap);
186: PetscContainerSetUserDestroy(c,MatISContainerDestroyPtAP_Private);
187: PetscObjectCompose((PetscObject)C,"_MatIS_PtAP",(PetscObject)c);
188: PetscContainerDestroy(&c);
189: ptap->fill = fill;
191: MatISGetLocalToGlobalMapping(A,&rl2g,&cl2g);
193: ISLocalToGlobalMappingGetBlockSize(cl2g,&ibs);
194: ISLocalToGlobalMappingGetSize(cl2g,&N);
195: ISLocalToGlobalMappingGetBlockIndices(cl2g,&garray);
196: ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris0);
197: ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&garray);
199: MatCreateSubMatrix(P,ptap->ris0,NULL,MAT_INITIAL_MATRIX,&PT);
200: MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis0);
201: ISLocalToGlobalMappingCreateIS(ptap->cis0,&Ccl2g);
202: MatDestroy(&PT);
204: Crl2g = NULL;
205: if (rl2g != cl2g) { /* unsymmetric A mapping */
206: PetscBool same,lsame = PETSC_FALSE;
207: PetscInt N1,ibs1;
209: ISLocalToGlobalMappingGetSize(rl2g,&N1);
210: ISLocalToGlobalMappingGetBlockSize(rl2g,&ibs1);
211: ISLocalToGlobalMappingGetBlockIndices(rl2g,&garray);
212: ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris1);
213: ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&garray);
214: if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
215: const PetscInt *i1,*i2;
217: ISBlockGetIndices(ptap->ris0,&i1);
218: ISBlockGetIndices(ptap->ris1,&i2);
219: PetscArraycmp(i1,i2,N,&lsame);
220: }
221: MPIU_Allreduce(&lsame,&same,1,MPIU_BOOL,MPI_LAND,comm);
222: if (same) {
223: ISDestroy(&ptap->ris1);
224: } else {
225: MatCreateSubMatrix(P,ptap->ris1,NULL,MAT_INITIAL_MATRIX,&PT);
226: MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis1);
227: ISLocalToGlobalMappingCreateIS(ptap->cis1,&Crl2g);
228: MatDestroy(&PT);
229: }
230: }
231: /* Not sure about this
232: if (!Crl2g) {
233: MatGetBlockSize(C,&ibs);
234: ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs);
235: }
236: */
237: MatSetLocalToGlobalMapping(C,Crl2g ? Crl2g : Ccl2g,Ccl2g);
238: ISLocalToGlobalMappingDestroy(&Crl2g);
239: ISLocalToGlobalMappingDestroy(&Ccl2g);
241: C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
242: return 0;
243: }
245: /* ----------------------------------------- */
246: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
247: {
248: Mat_Product *product = C->product;
249: Mat A=product->A,P=product->B;
250: PetscReal fill=product->fill;
252: MatPtAPSymbolic_IS_XAIJ(A,P,fill,C);
253: C->ops->productnumeric = MatProductNumeric_PtAP;
254: return 0;
255: }
257: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
258: {
259: C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
260: return 0;
261: }
263: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
264: {
265: Mat_Product *product = C->product;
267: if (product->type == MATPRODUCT_PtAP) {
268: MatProductSetFromOptions_IS_XAIJ_PtAP(C);
269: }
270: return 0;
271: }
273: /* ----------------------------------------- */
274: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
275: {
276: MatISLocalFields lf = (MatISLocalFields)ptr;
277: PetscInt i;
279: for (i=0;i<lf->nr;i++) {
280: ISDestroy(&lf->rf[i]);
281: }
282: for (i=0;i<lf->nc;i++) {
283: ISDestroy(&lf->cf[i]);
284: }
285: PetscFree2(lf->rf,lf->cf);
286: PetscFree(lf);
287: return 0;
288: }
290: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
291: {
292: Mat B,lB;
294: if (reuse != MAT_REUSE_MATRIX) {
295: ISLocalToGlobalMapping rl2g,cl2g;
296: PetscInt bs;
297: IS is;
299: MatGetBlockSize(A,&bs);
300: ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is);
301: if (bs > 1) {
302: IS is2;
303: PetscInt i,*aux;
305: ISGetLocalSize(is,&i);
306: ISGetIndices(is,(const PetscInt**)&aux);
307: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2);
308: ISRestoreIndices(is,(const PetscInt**)&aux);
309: ISDestroy(&is);
310: is = is2;
311: }
312: ISSetIdentity(is);
313: ISLocalToGlobalMappingCreateIS(is,&rl2g);
314: ISDestroy(&is);
315: ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->n/bs,0,1,&is);
316: if (bs > 1) {
317: IS is2;
318: PetscInt i,*aux;
320: ISGetLocalSize(is,&i);
321: ISGetIndices(is,(const PetscInt**)&aux);
322: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2);
323: ISRestoreIndices(is,(const PetscInt**)&aux);
324: ISDestroy(&is);
325: is = is2;
326: }
327: ISSetIdentity(is);
328: ISLocalToGlobalMappingCreateIS(is,&cl2g);
329: ISDestroy(&is);
330: MatCreateIS(PetscObjectComm((PetscObject)A),bs,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,rl2g,cl2g,&B);
331: ISLocalToGlobalMappingDestroy(&rl2g);
332: ISLocalToGlobalMappingDestroy(&cl2g);
333: MatDuplicate(A,MAT_COPY_VALUES,&lB);
334: if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
335: } else {
336: B = *newmat;
337: PetscObjectReference((PetscObject)A);
338: lB = A;
339: }
340: MatISSetLocalMat(B,lB);
341: MatDestroy(&lB);
342: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
343: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
344: if (reuse == MAT_INPLACE_MATRIX) {
345: MatHeaderReplace(A,&B);
346: }
347: return 0;
348: }
350: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
351: {
352: Mat_IS *matis = (Mat_IS*)(A->data);
353: PetscScalar *aa;
354: const PetscInt *ii,*jj;
355: PetscInt i,n,m;
356: PetscInt *ecount,**eneighs;
357: PetscBool flg;
359: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
361: ISLocalToGlobalMappingGetNodeInfo(matis->rmapping,&n,&ecount,&eneighs);
363: MatSeqAIJGetArray(matis->A,&aa);
364: for (i=0;i<n;i++) {
365: if (ecount[i] > 1) {
366: PetscInt j;
368: for (j=ii[i];j<ii[i+1];j++) {
369: PetscInt i2 = jj[j],p,p2;
370: PetscReal scal = 0.0;
372: for (p=0;p<ecount[i];p++) {
373: for (p2=0;p2<ecount[i2];p2++) {
374: if (eneighs[i][p] == eneighs[i2][p2]) { scal += 1.0; break; }
375: }
376: }
377: if (scal) aa[j] /= scal;
378: }
379: }
380: }
381: ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping,&n,&ecount,&eneighs);
382: MatSeqAIJRestoreArray(matis->A,&aa);
383: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
385: return 0;
386: }
388: typedef enum {MAT_IS_DISASSEMBLE_L2G_NATURAL,MAT_IS_DISASSEMBLE_L2G_MAT, MAT_IS_DISASSEMBLE_L2G_ND} MatISDisassemblel2gType;
390: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
391: {
392: Mat Ad,Ao;
393: IS is,ndmap,ndsub;
394: MPI_Comm comm;
395: const PetscInt *garray,*ndmapi;
396: PetscInt bs,i,cnt,nl,*ncount,*ndmapc;
397: PetscBool ismpiaij,ismpibaij;
398: const char *const MatISDisassemblel2gTypes[] = {"NATURAL","MAT","ND","MatISDisassemblel2gType","MAT_IS_DISASSEMBLE_L2G_",NULL};
399: MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
400: MatPartitioning part;
401: PetscSF sf;
402: PetscErrorCode ierr;
404: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatIS l2g disassembling options","Mat");
405: PetscOptionsEnum("-mat_is_disassemble_l2g_type","Type of local-to-global mapping to be used for disassembling","MatISDisassemblel2gType",MatISDisassemblel2gTypes,(PetscEnum)mode,(PetscEnum*)&mode,NULL);
406: PetscOptionsEnd();
407: if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
408: MatGetLocalToGlobalMapping(A,l2g,NULL);
409: return 0;
410: }
411: PetscObjectGetComm((PetscObject)A,&comm);
412: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij);
413: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij);
414: MatGetBlockSize(A,&bs);
415: switch (mode) {
416: case MAT_IS_DISASSEMBLE_L2G_ND:
417: MatPartitioningCreate(comm,&part);
418: MatPartitioningSetAdjacency(part,A);
419: PetscObjectSetOptionsPrefix((PetscObject)part,((PetscObject)A)->prefix);
420: MatPartitioningSetFromOptions(part);
421: MatPartitioningApplyND(part,&ndmap);
422: MatPartitioningDestroy(&part);
423: ISBuildTwoSided(ndmap,NULL,&ndsub);
424: MatMPIAIJSetUseScalableIncreaseOverlap(A,PETSC_TRUE);
425: MatIncreaseOverlap(A,1,&ndsub,1);
426: ISLocalToGlobalMappingCreateIS(ndsub,l2g);
428: /* it may happen that a separator node is not properly shared */
429: ISLocalToGlobalMappingGetNodeInfo(*l2g,&nl,&ncount,NULL);
430: PetscSFCreate(comm,&sf);
431: ISLocalToGlobalMappingGetIndices(*l2g,&garray);
432: PetscSFSetGraphLayout(sf,A->rmap,nl,NULL,PETSC_OWN_POINTER,garray);
433: ISLocalToGlobalMappingRestoreIndices(*l2g,&garray);
434: PetscCalloc1(A->rmap->n,&ndmapc);
435: PetscSFReduceBegin(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE);
436: PetscSFReduceEnd(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE);
437: ISLocalToGlobalMappingRestoreNodeInfo(*l2g,NULL,&ncount,NULL);
438: ISGetIndices(ndmap,&ndmapi);
439: for (i = 0, cnt = 0; i < A->rmap->n; i++)
440: if (ndmapi[i] < 0 && ndmapc[i] < 2)
441: cnt++;
443: MPIU_Allreduce(&cnt,&i,1,MPIU_INT,MPI_MAX,comm);
444: if (i) { /* we detected isolated separator nodes */
445: Mat A2,A3;
446: IS *workis,is2;
447: PetscScalar *vals;
448: PetscInt gcnt = i,*dnz,*onz,j,*lndmapi;
449: ISLocalToGlobalMapping ll2g;
450: PetscBool flg;
451: const PetscInt *ii,*jj;
453: /* communicate global id of separators */
454: MatPreallocateInitialize(comm,A->rmap->n,A->cmap->n,dnz,onz);
455: for (i = 0, cnt = 0; i < A->rmap->n; i++)
456: dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;
458: PetscMalloc1(nl,&lndmapi);
459: PetscSFBcastBegin(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE);
461: /* compute adjacency of isolated separators node */
462: PetscMalloc1(gcnt,&workis);
463: for (i = 0, cnt = 0; i < A->rmap->n; i++) {
464: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
465: ISCreateStride(comm,1,i+A->rmap->rstart,1,&workis[cnt++]);
466: }
467: }
468: for (i = cnt; i < gcnt; i++) {
469: ISCreateStride(comm,0,0,1,&workis[i]);
470: }
471: for (i = 0; i < gcnt; i++) {
472: PetscObjectSetName((PetscObject)workis[i],"ISOLATED");
473: ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
474: }
476: /* no communications since all the ISes correspond to locally owned rows */
477: MatIncreaseOverlap(A,gcnt,workis,1);
479: /* end communicate global id of separators */
480: PetscSFBcastEnd(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE);
482: /* communicate new layers : create a matrix and transpose it */
483: PetscArrayzero(dnz,A->rmap->n);
484: PetscArrayzero(onz,A->rmap->n);
485: for (i = 0, j = 0; i < A->rmap->n; i++) {
486: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
487: const PetscInt* idxs;
488: PetscInt s;
490: ISGetLocalSize(workis[j],&s);
491: ISGetIndices(workis[j],&idxs);
492: MatPreallocateSet(i+A->rmap->rstart,s,idxs,dnz,onz);
493: j++;
494: }
495: }
498: for (i = 0; i < gcnt; i++) {
499: PetscObjectSetName((PetscObject)workis[i],"EXTENDED");
500: ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
501: }
503: for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j,dnz[i]+onz[i]);
504: PetscMalloc1(j,&vals);
505: for (i = 0; i < j; i++) vals[i] = 1.0;
507: MatCreate(comm,&A2);
508: MatSetType(A2,MATMPIAIJ);
509: MatSetSizes(A2,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
510: MatMPIAIJSetPreallocation(A2,0,dnz,0,onz);
511: MatSetOption(A2,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
512: for (i = 0, j = 0; i < A2->rmap->n; i++) {
513: PetscInt row = i+A2->rmap->rstart,s = dnz[i] + onz[i];
514: const PetscInt* idxs;
516: if (s) {
517: ISGetIndices(workis[j],&idxs);
518: MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES);
519: ISRestoreIndices(workis[j],&idxs);
520: j++;
521: }
522: }
524: PetscFree(vals);
525: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
526: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
527: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
529: /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
530: for (i = 0, j = 0; i < nl; i++)
531: if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
532: ISCreateGeneral(comm,j,lndmapi,PETSC_USE_POINTER,&is);
533: MatMPIAIJGetLocalMatCondensed(A2,MAT_INITIAL_MATRIX,&is,NULL,&A3);
534: ISDestroy(&is);
535: MatDestroy(&A2);
537: /* extend local to global map to include connected isolated separators */
538: PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is);
540: ISLocalToGlobalMappingCreateIS(is,&ll2g);
541: MatGetRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
543: ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is);
544: MatRestoreRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
546: ISLocalToGlobalMappingApplyIS(ll2g,is,&is2);
547: ISDestroy(&is);
548: ISLocalToGlobalMappingDestroy(&ll2g);
550: /* add new nodes to the local-to-global map */
551: ISLocalToGlobalMappingDestroy(l2g);
552: ISExpand(ndsub,is2,&is);
553: ISDestroy(&is2);
554: ISLocalToGlobalMappingCreateIS(is,l2g);
555: ISDestroy(&is);
557: MatDestroy(&A3);
558: PetscFree(lndmapi);
559: MatPreallocateFinalize(dnz,onz);
560: for (i = 0; i < gcnt; i++) {
561: ISDestroy(&workis[i]);
562: }
563: PetscFree(workis);
564: }
565: ISRestoreIndices(ndmap,&ndmapi);
566: PetscSFDestroy(&sf);
567: PetscFree(ndmapc);
568: ISDestroy(&ndmap);
569: ISDestroy(&ndsub);
570: ISLocalToGlobalMappingSetBlockSize(*l2g,bs);
571: ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view");
572: break;
573: case MAT_IS_DISASSEMBLE_L2G_NATURAL:
574: if (ismpiaij) {
575: MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
576: } else if (ismpibaij) {
577: MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
578: } else SETERRQ(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
580: if (A->rmap->n) {
581: PetscInt dc,oc,stc,*aux;
583: MatGetLocalSize(Ad,NULL,&dc);
584: MatGetLocalSize(Ao,NULL,&oc);
585: MatGetOwnershipRangeColumn(A,&stc,NULL);
586: PetscMalloc1((dc+oc)/bs,&aux);
587: for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs;
588: for (i=0; i<oc/bs; i++) aux[i+dc/bs] = (ismpiaij ? garray[i*bs]/bs : garray[i]);
589: ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
590: } else {
591: ISCreateBlock(comm,1,0,NULL,PETSC_OWN_POINTER,&is);
592: }
593: ISLocalToGlobalMappingCreateIS(is,l2g);
594: ISDestroy(&is);
595: break;
596: default:
597: SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Unsupported l2g disassembling type %d",mode);
598: }
599: return 0;
600: }
602: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
603: {
604: Mat lA,Ad,Ao,B = NULL;
605: ISLocalToGlobalMapping rl2g,cl2g;
606: IS is;
607: MPI_Comm comm;
608: void *ptrs[2];
609: const char *names[2] = {"_convert_csr_aux","_convert_csr_data"};
610: const PetscInt *garray;
611: PetscScalar *dd,*od,*aa,*data;
612: const PetscInt *di,*dj,*oi,*oj;
613: const PetscInt *odi,*odj,*ooi,*ooj;
614: PetscInt *aux,*ii,*jj;
615: PetscInt bs,lc,dr,dc,oc,str,stc,nnz,i,jd,jo,cum;
616: PetscBool flg,ismpiaij,ismpibaij,was_inplace = PETSC_FALSE;
617: PetscMPIInt size;
619: PetscObjectGetComm((PetscObject)A,&comm);
620: MPI_Comm_size(comm,&size);
621: if (size == 1) {
622: MatConvert_SeqXAIJ_IS(A,type,reuse,newmat);
623: return 0;
624: }
625: if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
626: MatMPIXAIJComputeLocalToGlobalMapping_Private(A,&rl2g);
627: MatCreate(comm,&B);
628: MatSetType(B,MATIS);
629: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
630: MatSetLocalToGlobalMapping(B,rl2g,rl2g);
631: MatGetBlockSize(A,&bs);
632: MatSetBlockSize(B,bs);
633: ISLocalToGlobalMappingDestroy(&rl2g);
634: if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
635: reuse = MAT_REUSE_MATRIX;
636: }
637: if (reuse == MAT_REUSE_MATRIX) {
638: Mat *newlA, lA;
639: IS rows, cols;
640: const PetscInt *ridx, *cidx;
641: PetscInt rbs, cbs, nr, nc;
643: if (!B) B = *newmat;
644: MatISGetLocalToGlobalMapping(B,&rl2g,&cl2g);
645: ISLocalToGlobalMappingGetBlockIndices(rl2g,&ridx);
646: ISLocalToGlobalMappingGetBlockIndices(cl2g,&cidx);
647: ISLocalToGlobalMappingGetSize(rl2g,&nr);
648: ISLocalToGlobalMappingGetSize(cl2g,&nc);
649: ISLocalToGlobalMappingGetBlockSize(rl2g,&rbs);
650: ISLocalToGlobalMappingGetBlockSize(cl2g,&cbs);
651: ISCreateBlock(comm,rbs,nr/rbs,ridx,PETSC_USE_POINTER,&rows);
652: if (rl2g != cl2g) {
653: ISCreateBlock(comm,cbs,nc/cbs,cidx,PETSC_USE_POINTER,&cols);
654: } else {
655: PetscObjectReference((PetscObject)rows);
656: cols = rows;
657: }
658: MatISGetLocalMat(B,&lA);
659: MatCreateSubMatrices(A,1,&rows,&cols,MAT_INITIAL_MATRIX,&newlA);
660: MatConvert(newlA[0],MATSEQAIJ,MAT_INPLACE_MATRIX,&newlA[0]);
661: ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&ridx);
662: ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&cidx);
663: ISDestroy(&rows);
664: ISDestroy(&cols);
665: if (!lA->preallocated) { /* first time */
666: MatDuplicate(newlA[0],MAT_COPY_VALUES,&lA);
667: MatISSetLocalMat(B,lA);
668: PetscObjectDereference((PetscObject)lA);
669: }
670: MatCopy(newlA[0],lA,SAME_NONZERO_PATTERN);
671: MatDestroySubMatrices(1,&newlA);
672: MatISScaleDisassembling_Private(B);
673: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
674: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
675: if (was_inplace) MatHeaderReplace(A,&B);
676: else *newmat = B;
677: return 0;
678: }
679: /* rectangular case, just compress out the column space */
680: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij);
681: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij);
682: if (ismpiaij) {
683: bs = 1;
684: MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
685: } else if (ismpibaij) {
686: MatGetBlockSize(A,&bs);
687: MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
688: MatConvert(Ad,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ad);
689: MatConvert(Ao,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ao);
690: } else SETERRQ(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
691: MatSeqAIJGetArray(Ad,&dd);
692: MatSeqAIJGetArray(Ao,&od);
695: /* access relevant information from MPIAIJ */
696: MatGetOwnershipRange(A,&str,NULL);
697: MatGetOwnershipRangeColumn(A,&stc,NULL);
698: MatGetLocalSize(A,&dr,&dc);
699: MatGetLocalSize(Ao,NULL,&oc);
700: MatGetRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&di,&dj,&flg);
702: MatGetRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&oi,&oj,&flg);
704: nnz = di[dr] + oi[dr];
705: /* store original pointers to be restored later */
706: odi = di; odj = dj; ooi = oi; ooj = oj;
708: /* generate l2g maps for rows and cols */
709: ISCreateStride(comm,dr/bs,str/bs,1,&is);
710: if (bs > 1) {
711: IS is2;
713: ISGetLocalSize(is,&i);
714: ISGetIndices(is,(const PetscInt**)&aux);
715: ISCreateBlock(comm,bs,i,aux,PETSC_COPY_VALUES,&is2);
716: ISRestoreIndices(is,(const PetscInt**)&aux);
717: ISDestroy(&is);
718: is = is2;
719: }
720: ISLocalToGlobalMappingCreateIS(is,&rl2g);
721: ISDestroy(&is);
722: if (dr) {
723: PetscMalloc1((dc+oc)/bs,&aux);
724: for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs;
725: for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
726: ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
727: lc = dc+oc;
728: } else {
729: ISCreateBlock(comm,bs,0,NULL,PETSC_OWN_POINTER,&is);
730: lc = 0;
731: }
732: ISLocalToGlobalMappingCreateIS(is,&cl2g);
733: ISDestroy(&is);
735: /* create MATIS object */
736: MatCreate(comm,&B);
737: MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
738: MatSetType(B,MATIS);
739: MatSetBlockSize(B,bs);
740: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
741: ISLocalToGlobalMappingDestroy(&rl2g);
742: ISLocalToGlobalMappingDestroy(&cl2g);
744: /* merge local matrices */
745: PetscMalloc1(nnz+dr+1,&aux);
746: PetscMalloc1(nnz,&data);
747: ii = aux;
748: jj = aux+dr+1;
749: aa = data;
750: *ii = *(di++) + *(oi++);
751: for (jd=0,jo=0,cum=0;*ii<nnz;cum++)
752: {
753: for (;jd<*di;jd++) { *jj++ = *dj++; *aa++ = *dd++; }
754: for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; }
755: *(++ii) = *(di++) + *(oi++);
756: }
757: for (;cum<dr;cum++) *(++ii) = nnz;
759: MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg);
761: MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg);
763: MatSeqAIJRestoreArray(Ad,&dd);
764: MatSeqAIJRestoreArray(Ao,&od);
766: ii = aux;
767: jj = aux+dr+1;
768: aa = data;
769: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);
771: /* create containers to destroy the data */
772: ptrs[0] = aux;
773: ptrs[1] = data;
774: for (i=0; i<2; i++) {
775: PetscContainer c;
777: PetscContainerCreate(PETSC_COMM_SELF,&c);
778: PetscContainerSetPointer(c,ptrs[i]);
779: PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);
780: PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
781: PetscContainerDestroy(&c);
782: }
783: if (ismpibaij) { /* destroy converted local matrices */
784: MatDestroy(&Ad);
785: MatDestroy(&Ao);
786: }
788: /* finalize matrix */
789: MatISSetLocalMat(B,lA);
790: MatDestroy(&lA);
791: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
792: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
793: if (reuse == MAT_INPLACE_MATRIX) {
794: MatHeaderReplace(A,&B);
795: } else *newmat = B;
796: return 0;
797: }
799: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
800: {
801: Mat **nest,*snest,**rnest,lA,B;
802: IS *iscol,*isrow,*islrow,*islcol;
803: ISLocalToGlobalMapping rl2g,cl2g;
804: MPI_Comm comm;
805: PetscInt *lr,*lc,*l2gidxs;
806: PetscInt i,j,nr,nc,rbs,cbs;
807: PetscBool convert,lreuse,*istrans;
809: MatNestGetSubMats(A,&nr,&nc,&nest);
810: lreuse = PETSC_FALSE;
811: rnest = NULL;
812: if (reuse == MAT_REUSE_MATRIX) {
813: PetscBool ismatis,isnest;
815: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
817: MatISGetLocalMat(*newmat,&lA);
818: PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
819: if (isnest) {
820: MatNestGetSubMats(lA,&i,&j,&rnest);
821: lreuse = (PetscBool)(i == nr && j == nc);
822: if (!lreuse) rnest = NULL;
823: }
824: }
825: PetscObjectGetComm((PetscObject)A,&comm);
826: PetscCalloc2(nr,&lr,nc,&lc);
827: PetscCalloc6(nr,&isrow,nc,&iscol,nr,&islrow,nc,&islcol,nr*nc,&snest,nr*nc,&istrans);
828: MatNestGetISs(A,isrow,iscol);
829: for (i=0;i<nr;i++) {
830: for (j=0;j<nc;j++) {
831: PetscBool ismatis;
832: PetscInt l1,l2,lb1,lb2,ij=i*nc+j;
834: /* Null matrix pointers are allowed in MATNEST */
835: if (!nest[i][j]) continue;
837: /* Nested matrices should be of type MATIS */
838: PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
839: if (istrans[ij]) {
840: Mat T,lT;
841: MatTransposeGetMat(nest[i][j],&T);
842: PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
844: MatISGetLocalMat(T,&lT);
845: MatCreateTranspose(lT,&snest[ij]);
846: } else {
847: PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
849: MatISGetLocalMat(nest[i][j],&snest[ij]);
850: }
852: /* Check compatibility of local sizes */
853: MatGetSize(snest[ij],&l1,&l2);
854: MatGetBlockSizes(snest[ij],&lb1,&lb2);
855: if (!l1 || !l2) continue;
858: lr[i] = l1;
859: lc[j] = l2;
861: /* check compatibilty for local matrix reusage */
862: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
863: }
864: }
866: if (PetscDefined (USE_DEBUG)) {
867: /* Check compatibility of l2g maps for rows */
868: for (i=0;i<nr;i++) {
869: rl2g = NULL;
870: for (j=0;j<nc;j++) {
871: PetscInt n1,n2;
873: if (!nest[i][j]) continue;
874: if (istrans[i*nc+j]) {
875: Mat T;
877: MatTransposeGetMat(nest[i][j],&T);
878: MatISGetLocalToGlobalMapping(T,NULL,&cl2g);
879: } else {
880: MatISGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
881: }
882: ISLocalToGlobalMappingGetSize(cl2g,&n1);
883: if (!n1) continue;
884: if (!rl2g) {
885: rl2g = cl2g;
886: } else {
887: const PetscInt *idxs1,*idxs2;
888: PetscBool same;
890: ISLocalToGlobalMappingGetSize(rl2g,&n2);
892: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
893: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
894: PetscArraycmp(idxs1,idxs2,n1,&same);
895: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
896: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
898: }
899: }
900: }
901: /* Check compatibility of l2g maps for columns */
902: for (i=0;i<nc;i++) {
903: rl2g = NULL;
904: for (j=0;j<nr;j++) {
905: PetscInt n1,n2;
907: if (!nest[j][i]) continue;
908: if (istrans[j*nc+i]) {
909: Mat T;
911: MatTransposeGetMat(nest[j][i],&T);
912: MatISGetLocalToGlobalMapping(T,&cl2g,NULL);
913: } else {
914: MatISGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
915: }
916: ISLocalToGlobalMappingGetSize(cl2g,&n1);
917: if (!n1) continue;
918: if (!rl2g) {
919: rl2g = cl2g;
920: } else {
921: const PetscInt *idxs1,*idxs2;
922: PetscBool same;
924: ISLocalToGlobalMappingGetSize(rl2g,&n2);
926: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
927: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
928: PetscArraycmp(idxs1,idxs2,n1,&same);
929: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
930: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
932: }
933: }
934: }
935: }
937: B = NULL;
938: if (reuse != MAT_REUSE_MATRIX) {
939: PetscInt stl;
941: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
942: for (i=0,stl=0;i<nr;i++) stl += lr[i];
943: PetscMalloc1(stl,&l2gidxs);
944: for (i=0,stl=0;i<nr;i++) {
945: Mat usedmat;
946: Mat_IS *matis;
947: const PetscInt *idxs;
949: /* local IS for local NEST */
950: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
952: /* l2gmap */
953: j = 0;
954: usedmat = nest[i][j];
955: while (!usedmat && j < nc-1) usedmat = nest[i][++j];
958: if (istrans[i*nc+j]) {
959: Mat T;
960: MatTransposeGetMat(usedmat,&T);
961: usedmat = T;
962: }
963: matis = (Mat_IS*)(usedmat->data);
964: ISGetIndices(isrow[i],&idxs);
965: if (istrans[i*nc+j]) {
966: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
967: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
968: } else {
969: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
970: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
971: }
972: ISRestoreIndices(isrow[i],&idxs);
973: stl += lr[i];
974: }
975: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);
977: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
978: for (i=0,stl=0;i<nc;i++) stl += lc[i];
979: PetscMalloc1(stl,&l2gidxs);
980: for (i=0,stl=0;i<nc;i++) {
981: Mat usedmat;
982: Mat_IS *matis;
983: const PetscInt *idxs;
985: /* local IS for local NEST */
986: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
988: /* l2gmap */
989: j = 0;
990: usedmat = nest[j][i];
991: while (!usedmat && j < nr-1) usedmat = nest[++j][i];
993: if (istrans[j*nc+i]) {
994: Mat T;
995: MatTransposeGetMat(usedmat,&T);
996: usedmat = T;
997: }
998: matis = (Mat_IS*)(usedmat->data);
999: ISGetIndices(iscol[i],&idxs);
1000: if (istrans[j*nc+i]) {
1001: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1002: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1003: } else {
1004: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1005: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE);
1006: }
1007: ISRestoreIndices(iscol[i],&idxs);
1008: stl += lc[i];
1009: }
1010: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);
1012: /* Create MATIS */
1013: MatCreate(comm,&B);
1014: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1015: MatGetBlockSizes(A,&rbs,&cbs);
1016: MatSetBlockSizes(B,rbs,cbs);
1017: MatSetType(B,MATIS);
1018: MatISSetLocalMatType(B,MATNEST);
1019: { /* hack : avoid setup of scatters */
1020: Mat_IS *matis = (Mat_IS*)(B->data);
1021: matis->islocalref = PETSC_TRUE;
1022: }
1023: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
1024: ISLocalToGlobalMappingDestroy(&rl2g);
1025: ISLocalToGlobalMappingDestroy(&cl2g);
1026: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1027: MatNestSetVecType(lA,VECNEST);
1028: for (i=0;i<nr*nc;i++) {
1029: if (istrans[i]) {
1030: MatDestroy(&snest[i]);
1031: }
1032: }
1033: MatISSetLocalMat(B,lA);
1034: MatDestroy(&lA);
1035: { /* hack : setup of scatters done here */
1036: Mat_IS *matis = (Mat_IS*)(B->data);
1038: matis->islocalref = PETSC_FALSE;
1039: MatISSetUpScatters_Private(B);
1040: }
1041: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1042: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1043: if (reuse == MAT_INPLACE_MATRIX) {
1044: MatHeaderReplace(A,&B);
1045: } else {
1046: *newmat = B;
1047: }
1048: } else {
1049: if (lreuse) {
1050: MatISGetLocalMat(*newmat,&lA);
1051: for (i=0;i<nr;i++) {
1052: for (j=0;j<nc;j++) {
1053: if (snest[i*nc+j]) {
1054: MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
1055: if (istrans[i*nc+j]) {
1056: MatDestroy(&snest[i*nc+j]);
1057: }
1058: }
1059: }
1060: }
1061: } else {
1062: PetscInt stl;
1063: for (i=0,stl=0;i<nr;i++) {
1064: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
1065: stl += lr[i];
1066: }
1067: for (i=0,stl=0;i<nc;i++) {
1068: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
1069: stl += lc[i];
1070: }
1071: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1072: for (i=0;i<nr*nc;i++) {
1073: if (istrans[i]) {
1074: MatDestroy(&snest[i]);
1075: }
1076: }
1077: MatISSetLocalMat(*newmat,lA);
1078: MatDestroy(&lA);
1079: }
1080: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1081: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1082: }
1084: /* Create local matrix in MATNEST format */
1085: convert = PETSC_FALSE;
1086: PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
1087: if (convert) {
1088: Mat M;
1089: MatISLocalFields lf;
1090: PetscContainer c;
1092: MatISGetLocalMat(*newmat,&lA);
1093: MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
1094: MatISSetLocalMat(*newmat,M);
1095: MatDestroy(&M);
1097: /* attach local fields to the matrix */
1098: PetscNew(&lf);
1099: PetscMalloc2(nr,&lf->rf,nc,&lf->cf);
1100: for (i=0;i<nr;i++) {
1101: PetscInt n,st;
1103: ISGetLocalSize(islrow[i],&n);
1104: ISStrideGetInfo(islrow[i],&st,NULL);
1105: ISCreateStride(comm,n,st,1,&lf->rf[i]);
1106: }
1107: for (i=0;i<nc;i++) {
1108: PetscInt n,st;
1110: ISGetLocalSize(islcol[i],&n);
1111: ISStrideGetInfo(islcol[i],&st,NULL);
1112: ISCreateStride(comm,n,st,1,&lf->cf[i]);
1113: }
1114: lf->nr = nr;
1115: lf->nc = nc;
1116: PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
1117: PetscContainerSetPointer(c,lf);
1118: PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
1119: PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
1120: PetscContainerDestroy(&c);
1121: }
1123: /* Free workspace */
1124: for (i=0;i<nr;i++) {
1125: ISDestroy(&islrow[i]);
1126: }
1127: for (i=0;i<nc;i++) {
1128: ISDestroy(&islcol[i]);
1129: }
1130: PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
1131: PetscFree2(lr,lc);
1132: return 0;
1133: }
1135: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1136: {
1137: Mat_IS *matis = (Mat_IS*)A->data;
1138: Vec ll,rr;
1139: const PetscScalar *Y,*X;
1140: PetscScalar *x,*y;
1142: if (l) {
1143: ll = matis->y;
1144: VecGetArrayRead(l,&Y);
1145: VecGetArray(ll,&y);
1146: PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE);
1147: } else {
1148: ll = NULL;
1149: }
1150: if (r) {
1151: rr = matis->x;
1152: VecGetArrayRead(r,&X);
1153: VecGetArray(rr,&x);
1154: PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE);
1155: } else {
1156: rr = NULL;
1157: }
1158: if (ll) {
1159: PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE);
1160: VecRestoreArrayRead(l,&Y);
1161: VecRestoreArray(ll,&y);
1162: }
1163: if (rr) {
1164: PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE);
1165: VecRestoreArrayRead(r,&X);
1166: VecRestoreArray(rr,&x);
1167: }
1168: MatDiagonalScale(matis->A,ll,rr);
1169: return 0;
1170: }
1172: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
1173: {
1174: Mat_IS *matis = (Mat_IS*)A->data;
1175: MatInfo info;
1176: PetscLogDouble isend[6],irecv[6];
1177: PetscInt bs;
1179: MatGetBlockSize(A,&bs);
1180: if (matis->A->ops->getinfo) {
1181: MatGetInfo(matis->A,MAT_LOCAL,&info);
1182: isend[0] = info.nz_used;
1183: isend[1] = info.nz_allocated;
1184: isend[2] = info.nz_unneeded;
1185: isend[3] = info.memory;
1186: isend[4] = info.mallocs;
1187: } else {
1188: isend[0] = 0.;
1189: isend[1] = 0.;
1190: isend[2] = 0.;
1191: isend[3] = 0.;
1192: isend[4] = 0.;
1193: }
1194: isend[5] = matis->A->num_ass;
1195: if (flag == MAT_LOCAL) {
1196: ginfo->nz_used = isend[0];
1197: ginfo->nz_allocated = isend[1];
1198: ginfo->nz_unneeded = isend[2];
1199: ginfo->memory = isend[3];
1200: ginfo->mallocs = isend[4];
1201: ginfo->assemblies = isend[5];
1202: } else if (flag == MAT_GLOBAL_MAX) {
1203: MPIU_Allreduce(isend,irecv,6,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)A));
1205: ginfo->nz_used = irecv[0];
1206: ginfo->nz_allocated = irecv[1];
1207: ginfo->nz_unneeded = irecv[2];
1208: ginfo->memory = irecv[3];
1209: ginfo->mallocs = irecv[4];
1210: ginfo->assemblies = irecv[5];
1211: } else if (flag == MAT_GLOBAL_SUM) {
1212: MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)A));
1214: ginfo->nz_used = irecv[0];
1215: ginfo->nz_allocated = irecv[1];
1216: ginfo->nz_unneeded = irecv[2];
1217: ginfo->memory = irecv[3];
1218: ginfo->mallocs = irecv[4];
1219: ginfo->assemblies = A->num_ass;
1220: }
1221: ginfo->block_size = bs;
1222: ginfo->fill_ratio_given = 0;
1223: ginfo->fill_ratio_needed = 0;
1224: ginfo->factor_mallocs = 0;
1225: return 0;
1226: }
1228: static PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
1229: {
1230: Mat C,lC,lA;
1232: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1233: ISLocalToGlobalMapping rl2g,cl2g;
1234: MatCreate(PetscObjectComm((PetscObject)A),&C);
1235: MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
1236: MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
1237: MatSetType(C,MATIS);
1238: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
1239: MatSetLocalToGlobalMapping(C,cl2g,rl2g);
1240: } else C = *B;
1242: /* perform local transposition */
1243: MatISGetLocalMat(A,&lA);
1244: MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
1245: MatSetLocalToGlobalMapping(lC,lA->cmap->mapping,lA->rmap->mapping);
1246: MatISSetLocalMat(C,lC);
1247: MatDestroy(&lC);
1249: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1250: *B = C;
1251: } else {
1252: MatHeaderMerge(A,&C);
1253: }
1254: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1255: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1256: return 0;
1257: }
1259: static PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
1260: {
1261: Mat_IS *is = (Mat_IS*)A->data;
1263: if (D) { /* MatShift_IS pass D = NULL */
1264: VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1265: VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1266: }
1267: VecPointwiseDivide(is->y,is->y,is->counter);
1268: MatDiagonalSet(is->A,is->y,insmode);
1269: return 0;
1270: }
1272: static PetscErrorCode MatShift_IS(Mat A,PetscScalar a)
1273: {
1274: Mat_IS *is = (Mat_IS*)A->data;
1276: VecSet(is->y,a);
1277: MatDiagonalSet_IS(A,NULL,ADD_VALUES);
1278: return 0;
1279: }
1281: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1282: {
1283: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1286: ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
1287: ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
1288: MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1289: return 0;
1290: }
1292: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1293: {
1294: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1297: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
1298: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
1299: MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1300: return 0;
1301: }
1303: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
1304: {
1305: Mat locmat,newlocmat;
1306: Mat_IS *newmatis;
1307: const PetscInt *idxs;
1308: PetscInt i,m,n;
1310: if (scall == MAT_REUSE_MATRIX) {
1311: PetscBool ismatis;
1313: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
1315: newmatis = (Mat_IS*)(*newmat)->data;
1318: }
1319: /* irow and icol may not have duplicate entries */
1320: if (PetscDefined(USE_DEBUG)) {
1321: Vec rtest,ltest;
1322: const PetscScalar *array;
1324: MatCreateVecs(mat,<est,&rtest);
1325: ISGetLocalSize(irow,&n);
1326: ISGetIndices(irow,&idxs);
1327: for (i=0;i<n;i++) {
1328: VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
1329: }
1330: VecAssemblyBegin(rtest);
1331: VecAssemblyEnd(rtest);
1332: VecGetLocalSize(rtest,&n);
1333: VecGetOwnershipRange(rtest,&m,NULL);
1334: VecGetArrayRead(rtest,&array);
1336: VecRestoreArrayRead(rtest,&array);
1337: ISRestoreIndices(irow,&idxs);
1338: ISGetLocalSize(icol,&n);
1339: ISGetIndices(icol,&idxs);
1340: for (i=0;i<n;i++) {
1341: VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
1342: }
1343: VecAssemblyBegin(ltest);
1344: VecAssemblyEnd(ltest);
1345: VecGetLocalSize(ltest,&n);
1346: VecGetOwnershipRange(ltest,&m,NULL);
1347: VecGetArrayRead(ltest,&array);
1349: VecRestoreArrayRead(ltest,&array);
1350: ISRestoreIndices(icol,&idxs);
1351: VecDestroy(&rtest);
1352: VecDestroy(<est);
1353: }
1354: if (scall == MAT_INITIAL_MATRIX) {
1355: Mat_IS *matis = (Mat_IS*)mat->data;
1356: ISLocalToGlobalMapping rl2g;
1357: IS is;
1358: PetscInt *lidxs,*lgidxs,*newgidxs;
1359: PetscInt ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
1360: PetscBool cong;
1361: MPI_Comm comm;
1363: PetscObjectGetComm((PetscObject)mat,&comm);
1364: MatGetBlockSizes(mat,&arbs,&acbs);
1365: ISGetBlockSize(irow,&irbs);
1366: ISGetBlockSize(icol,&icbs);
1367: rbs = arbs == irbs ? irbs : 1;
1368: cbs = acbs == icbs ? icbs : 1;
1369: ISGetLocalSize(irow,&m);
1370: ISGetLocalSize(icol,&n);
1371: MatCreate(comm,newmat);
1372: MatSetType(*newmat,MATIS);
1373: MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
1374: MatSetBlockSizes(*newmat,rbs,cbs);
1375: /* communicate irow to their owners in the layout */
1376: ISGetIndices(irow,&idxs);
1377: PetscLayoutMapLocal(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
1378: ISRestoreIndices(irow,&idxs);
1379: PetscArrayzero(matis->sf_rootdata,matis->sf->nroots);
1380: for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
1381: PetscFree(lidxs);
1382: PetscFree(lgidxs);
1383: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1384: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1385: for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
1386: PetscMalloc1(newloc,&newgidxs);
1387: PetscMalloc1(newloc,&lidxs);
1388: for (i=0,newloc=0;i<matis->sf->nleaves;i++)
1389: if (matis->sf_leafdata[i]) {
1390: lidxs[newloc] = i;
1391: newgidxs[newloc++] = matis->sf_leafdata[i]-1;
1392: }
1393: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1394: ISLocalToGlobalMappingCreateIS(is,&rl2g);
1395: ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);
1396: ISDestroy(&is);
1397: /* local is to extract local submatrix */
1398: newmatis = (Mat_IS*)(*newmat)->data;
1399: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
1400: MatHasCongruentLayouts(mat,&cong);
1401: if (cong && irow == icol && matis->csf == matis->sf) {
1402: MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
1403: PetscObjectReference((PetscObject)newmatis->getsub_ris);
1404: newmatis->getsub_cis = newmatis->getsub_ris;
1405: } else {
1406: ISLocalToGlobalMapping cl2g;
1408: /* communicate icol to their owners in the layout */
1409: ISGetIndices(icol,&idxs);
1410: PetscLayoutMapLocal(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
1411: ISRestoreIndices(icol,&idxs);
1412: PetscArrayzero(matis->csf_rootdata,matis->csf->nroots);
1413: for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
1414: PetscFree(lidxs);
1415: PetscFree(lgidxs);
1416: PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE);
1417: PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE);
1418: for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
1419: PetscMalloc1(newloc,&newgidxs);
1420: PetscMalloc1(newloc,&lidxs);
1421: for (i=0,newloc=0;i<matis->csf->nleaves;i++)
1422: if (matis->csf_leafdata[i]) {
1423: lidxs[newloc] = i;
1424: newgidxs[newloc++] = matis->csf_leafdata[i]-1;
1425: }
1426: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1427: ISLocalToGlobalMappingCreateIS(is,&cl2g);
1428: ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);
1429: ISDestroy(&is);
1430: /* local is to extract local submatrix */
1431: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
1432: MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
1433: ISLocalToGlobalMappingDestroy(&cl2g);
1434: }
1435: ISLocalToGlobalMappingDestroy(&rl2g);
1436: } else {
1437: MatISGetLocalMat(*newmat,&newlocmat);
1438: }
1439: MatISGetLocalMat(mat,&locmat);
1440: newmatis = (Mat_IS*)(*newmat)->data;
1441: MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
1442: if (scall == MAT_INITIAL_MATRIX) {
1443: MatISSetLocalMat(*newmat,newlocmat);
1444: MatDestroy(&newlocmat);
1445: }
1446: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1447: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1448: return 0;
1449: }
1451: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
1452: {
1453: Mat_IS *a = (Mat_IS*)A->data,*b;
1454: PetscBool ismatis;
1456: PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
1458: b = (Mat_IS*)B->data;
1459: MatCopy(a->A,b->A,str);
1460: PetscObjectStateIncrease((PetscObject)B);
1461: return 0;
1462: }
1464: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d)
1465: {
1466: Vec v;
1467: const PetscScalar *array;
1468: PetscInt i,n;
1470: *missing = PETSC_FALSE;
1471: MatCreateVecs(A,NULL,&v);
1472: MatGetDiagonal(A,v);
1473: VecGetLocalSize(v,&n);
1474: VecGetArrayRead(v,&array);
1475: for (i=0;i<n;i++) if (array[i] == 0.) break;
1476: VecRestoreArrayRead(v,&array);
1477: VecDestroy(&v);
1478: if (i != n) *missing = PETSC_TRUE;
1479: if (d) {
1480: *d = -1;
1481: if (*missing) {
1482: PetscInt rstart;
1483: MatGetOwnershipRange(A,&rstart,NULL);
1484: *d = i+rstart;
1485: }
1486: }
1487: return 0;
1488: }
1490: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1491: {
1492: Mat_IS *matis = (Mat_IS*)(B->data);
1493: const PetscInt *gidxs;
1494: PetscInt nleaves;
1496: if (matis->sf) return 0;
1497: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
1498: ISLocalToGlobalMappingGetIndices(matis->rmapping,&gidxs);
1499: ISLocalToGlobalMappingGetSize(matis->rmapping,&nleaves);
1500: PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1501: ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&gidxs);
1502: PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
1503: if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1504: ISLocalToGlobalMappingGetSize(matis->cmapping,&nleaves);
1505: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
1506: ISLocalToGlobalMappingGetIndices(matis->cmapping,&gidxs);
1507: PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1508: ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&gidxs);
1509: PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
1510: } else {
1511: matis->csf = matis->sf;
1512: matis->csf_leafdata = matis->sf_leafdata;
1513: matis->csf_rootdata = matis->sf_rootdata;
1514: }
1515: return 0;
1516: }
1518: /*@
1519: MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP.
1521: Collective
1523: Input Parameters:
1524: + A - the matrix
1525: - store - the boolean flag
1527: Level: advanced
1529: Notes:
1531: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
1532: @*/
1533: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1534: {
1538: PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
1539: return 0;
1540: }
1542: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1543: {
1544: Mat_IS *matis = (Mat_IS*)(A->data);
1546: matis->storel2l = store;
1547: if (!store) {
1548: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
1549: }
1550: return 0;
1551: }
1553: /*@
1554: MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1556: Collective
1558: Input Parameters:
1559: + A - the matrix
1560: - fix - the boolean flag
1562: Level: advanced
1564: Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process.
1566: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
1567: @*/
1568: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1569: {
1573: PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
1574: return 0;
1575: }
1577: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1578: {
1579: Mat_IS *matis = (Mat_IS*)(A->data);
1581: matis->locempty = fix;
1582: return 0;
1583: }
1585: /*@
1586: MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.
1588: Collective
1590: Input Parameters:
1591: + B - the matrix
1592: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1593: (same value is used for all local rows)
1594: . d_nnz - array containing the number of nonzeros in the various rows of the
1595: DIAGONAL portion of the local submatrix (possibly different for each row)
1596: or NULL, if d_nz is used to specify the nonzero structure.
1597: The size of this array is equal to the number of local rows, i.e 'm'.
1598: For matrices that will be factored, you must leave room for (and set)
1599: the diagonal entry even if it is zero.
1600: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1601: submatrix (same value is used for all local rows).
1602: - o_nnz - array containing the number of nonzeros in the various rows of the
1603: OFF-DIAGONAL portion of the local submatrix (possibly different for
1604: each row) or NULL, if o_nz is used to specify the nonzero
1605: structure. The size of this array is equal to the number
1606: of local rows, i.e 'm'.
1608: If the *_nnz parameter is given then the *_nz parameter is ignored
1610: Level: intermediate
1612: Notes:
1613: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
1614: from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
1615: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1617: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1618: @*/
1619: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1620: {
1623: PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1624: return 0;
1625: }
1627: /* this is used by DMDA */
1628: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1629: {
1630: Mat_IS *matis = (Mat_IS*)(B->data);
1631: PetscInt bs,i,nlocalcols;
1633: MatSetUp(B);
1634: if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
1635: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
1637: if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz;
1638: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
1640: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1641: MatGetSize(matis->A,NULL,&nlocalcols);
1642: MatGetBlockSize(matis->A,&bs);
1643: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1645: for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1646: MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1647: #if defined(PETSC_HAVE_HYPRE)
1648: MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1649: #endif
1651: for (i=0;i<matis->sf->nleaves/bs;i++) {
1652: PetscInt b;
1654: matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
1655: for (b=1;b<bs;b++) {
1656: matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i],matis->sf_leafdata[i*bs+b]/bs);
1657: }
1658: }
1659: MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1661: nlocalcols /= bs;
1662: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols - i);
1663: MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1665: /* for other matrix types */
1666: MatSetUp(matis->A);
1667: return 0;
1668: }
1670: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1671: {
1672: Mat_IS *matis = (Mat_IS*)(A->data);
1673: PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1674: const PetscInt *global_indices_r,*global_indices_c;
1675: PetscInt i,j,bs,rows,cols;
1676: PetscInt lrows,lcols;
1677: PetscInt local_rows,local_cols;
1678: PetscMPIInt size;
1679: PetscBool isdense,issbaij;
1680: PetscErrorCode ierr;
1682: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1683: MatGetSize(A,&rows,&cols);
1684: MatGetBlockSize(A,&bs);
1685: MatGetSize(matis->A,&local_rows,&local_cols);
1686: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1687: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1688: ISLocalToGlobalMappingGetIndices(matis->rmapping,&global_indices_r);
1689: if (matis->rmapping != matis->cmapping) {
1690: ISLocalToGlobalMappingGetIndices(matis->cmapping,&global_indices_c);
1691: } else global_indices_c = global_indices_r;
1693: if (issbaij) MatGetRowUpperTriangular(matis->A);
1694: /*
1695: An SF reduce is needed to sum up properly on shared rows.
1696: Note that generally preallocation is not exact, since it overestimates nonzeros
1697: */
1698: MatGetLocalSize(A,&lrows,&lcols);
1699: MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1700: /* All processes need to compute entire row ownership */
1701: PetscMalloc1(rows,&row_ownership);
1702: MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1703: for (i=0;i<size;i++) {
1704: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) row_ownership[j] = i;
1705: }
1706: MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);
1708: /*
1709: my_dnz and my_onz contains exact contribution to preallocation from each local mat
1710: then, they will be summed up properly. This way, preallocation is always sufficient
1711: */
1712: PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1713: /* preallocation as a MATAIJ */
1714: if (isdense) { /* special case for dense local matrices */
1715: for (i=0;i<local_rows;i++) {
1716: PetscInt owner = row_ownership[global_indices_r[i]];
1717: for (j=0;j<local_cols;j++) {
1718: PetscInt index_col = global_indices_c[j];
1719: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
1720: my_dnz[i] += 1;
1721: } else { /* offdiag block */
1722: my_onz[i] += 1;
1723: }
1724: }
1725: }
1726: } else if (matis->A->ops->getrowij) {
1727: const PetscInt *ii,*jj,*jptr;
1728: PetscBool done;
1729: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1731: jptr = jj;
1732: for (i=0;i<local_rows;i++) {
1733: PetscInt index_row = global_indices_r[i];
1734: for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1735: PetscInt owner = row_ownership[index_row];
1736: PetscInt index_col = global_indices_c[*jptr];
1737: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
1738: my_dnz[i] += 1;
1739: } else { /* offdiag block */
1740: my_onz[i] += 1;
1741: }
1742: /* same as before, interchanging rows and cols */
1743: if (issbaij && index_col != index_row) {
1744: owner = row_ownership[index_col];
1745: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) {
1746: my_dnz[*jptr] += 1;
1747: } else {
1748: my_onz[*jptr] += 1;
1749: }
1750: }
1751: }
1752: }
1753: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1755: } else { /* loop over rows and use MatGetRow */
1756: for (i=0;i<local_rows;i++) {
1757: const PetscInt *cols;
1758: PetscInt ncols,index_row = global_indices_r[i];
1759: MatGetRow(matis->A,i,&ncols,&cols,NULL);
1760: for (j=0;j<ncols;j++) {
1761: PetscInt owner = row_ownership[index_row];
1762: PetscInt index_col = global_indices_c[cols[j]];
1763: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
1764: my_dnz[i] += 1;
1765: } else { /* offdiag block */
1766: my_onz[i] += 1;
1767: }
1768: /* same as before, interchanging rows and cols */
1769: if (issbaij && index_col != index_row) {
1770: owner = row_ownership[index_col];
1771: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) {
1772: my_dnz[cols[j]] += 1;
1773: } else {
1774: my_onz[cols[j]] += 1;
1775: }
1776: }
1777: }
1778: MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1779: }
1780: }
1781: if (global_indices_c != global_indices_r) {
1782: ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&global_indices_c);
1783: }
1784: ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&global_indices_r);
1785: PetscFree(row_ownership);
1787: /* Reduce my_dnz and my_onz */
1788: if (maxreduce) {
1789: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1790: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1791: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1792: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1793: } else {
1794: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1795: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1796: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1797: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1798: }
1799: PetscFree2(my_dnz,my_onz);
1801: /* Resize preallocation if overestimated */
1802: for (i=0;i<lrows;i++) {
1803: dnz[i] = PetscMin(dnz[i],lcols);
1804: onz[i] = PetscMin(onz[i],cols-lcols);
1805: }
1807: /* Set preallocation */
1808: MatSetBlockSizesFromMats(B,A,A);
1809: MatSeqAIJSetPreallocation(B,0,dnz);
1810: MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1811: for (i=0;i<lrows;i+=bs) {
1812: PetscInt b, d = dnz[i],o = onz[i];
1814: for (b=1;b<bs;b++) {
1815: d = PetscMax(d,dnz[i+b]);
1816: o = PetscMax(o,onz[i+b]);
1817: }
1818: dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
1819: onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
1820: }
1821: MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1822: MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1823: MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1824: MatPreallocateFinalize(dnz,onz);
1825: if (issbaij) MatRestoreRowUpperTriangular(matis->A);
1826: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1827: return 0;
1828: }
1830: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1831: {
1832: Mat_IS *matis = (Mat_IS*)(mat->data);
1833: Mat local_mat,MT;
1834: PetscInt rbs,cbs,rows,cols,lrows,lcols;
1835: PetscInt local_rows,local_cols;
1836: PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij;
1837: PetscMPIInt size;
1838: const PetscScalar *array;
1840: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
1841: if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1842: Mat B;
1843: IS irows = NULL,icols = NULL;
1844: PetscInt rbs,cbs;
1846: ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs);
1847: ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs);
1848: if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1849: IS rows,cols;
1850: const PetscInt *ridxs,*cidxs;
1851: PetscInt i,nw,*work;
1853: ISLocalToGlobalMappingGetBlockIndices(matis->rmapping,&ridxs);
1854: ISLocalToGlobalMappingGetSize(matis->rmapping,&nw);
1855: nw = nw/rbs;
1856: PetscCalloc1(nw,&work);
1857: for (i=0;i<nw;i++) work[ridxs[i]] += 1;
1858: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
1859: if (i == nw) {
1860: ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);
1861: ISSetPermutation(rows);
1862: ISInvertPermutation(rows,PETSC_DECIDE,&irows);
1863: ISDestroy(&rows);
1864: }
1865: ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping,&ridxs);
1866: PetscFree(work);
1867: if (irows && matis->rmapping != matis->cmapping) {
1868: ISLocalToGlobalMappingGetBlockIndices(matis->cmapping,&cidxs);
1869: ISLocalToGlobalMappingGetSize(matis->cmapping,&nw);
1870: nw = nw/cbs;
1871: PetscCalloc1(nw,&work);
1872: for (i=0;i<nw;i++) work[cidxs[i]] += 1;
1873: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
1874: if (i == nw) {
1875: ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);
1876: ISSetPermutation(cols);
1877: ISInvertPermutation(cols,PETSC_DECIDE,&icols);
1878: ISDestroy(&cols);
1879: }
1880: ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping,&cidxs);
1881: PetscFree(work);
1882: } else if (irows) {
1883: PetscObjectReference((PetscObject)irows);
1884: icols = irows;
1885: }
1886: } else {
1887: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);
1888: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);
1889: if (irows) PetscObjectReference((PetscObject)irows);
1890: if (icols) PetscObjectReference((PetscObject)icols);
1891: }
1892: if (!irows || !icols) {
1893: ISDestroy(&icols);
1894: ISDestroy(&irows);
1895: goto general_assembly;
1896: }
1897: MatConvert(matis->A,mtype,MAT_INITIAL_MATRIX,&B);
1898: if (reuse != MAT_INPLACE_MATRIX) {
1899: MatCreateSubMatrix(B,irows,icols,reuse,M);
1900: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);
1901: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);
1902: } else {
1903: Mat C;
1905: MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
1906: MatHeaderReplace(mat,&C);
1907: }
1908: MatDestroy(&B);
1909: ISDestroy(&icols);
1910: ISDestroy(&irows);
1911: return 0;
1912: }
1913: general_assembly:
1914: MatGetSize(mat,&rows,&cols);
1915: ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs);
1916: ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs);
1917: MatGetLocalSize(mat,&lrows,&lcols);
1918: MatGetSize(matis->A,&local_rows,&local_cols);
1919: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
1920: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
1921: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
1922: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
1924: if (PetscDefined (USE_DEBUG)) {
1925: PetscBool lb[4],bb[4];
1927: lb[0] = isseqdense;
1928: lb[1] = isseqaij;
1929: lb[2] = isseqbaij;
1930: lb[3] = isseqsbaij;
1931: MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
1933: }
1935: if (reuse != MAT_REUSE_MATRIX) {
1936: MatCreate(PetscObjectComm((PetscObject)mat),&MT);
1937: MatSetSizes(MT,lrows,lcols,rows,cols);
1938: MatSetType(MT,mtype);
1939: MatSetBlockSizes(MT,rbs,cbs);
1940: MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
1941: } else {
1942: PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;
1944: /* some checks */
1945: MT = *M;
1946: MatGetBlockSizes(MT,&mrbs,&mcbs);
1947: MatGetSize(MT,&mrows,&mcols);
1948: MatGetLocalSize(MT,&mlrows,&mlcols);
1955: MatZeroEntries(MT);
1956: }
1958: if (isseqsbaij || isseqbaij) {
1959: MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
1960: isseqaij = PETSC_TRUE;
1961: } else {
1962: PetscObjectReference((PetscObject)matis->A);
1963: local_mat = matis->A;
1964: }
1966: /* Set values */
1967: MatSetLocalToGlobalMapping(MT,matis->rmapping,matis->cmapping);
1968: if (isseqdense) { /* special case for dense local matrices */
1969: PetscInt i,*dummy;
1971: PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
1972: for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
1973: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
1974: MatDenseGetArrayRead(local_mat,&array);
1975: MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
1976: MatDenseRestoreArrayRead(local_mat,&array);
1977: PetscFree(dummy);
1978: } else if (isseqaij) {
1979: const PetscInt *blocks;
1980: PetscInt i,nvtxs,*xadj,*adjncy, nb;
1981: PetscBool done;
1982: PetscScalar *sarray;
1984: MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1986: MatSeqAIJGetArray(local_mat,&sarray);
1987: MatGetVariableBlockSizes(local_mat,&nb,&blocks);
1988: if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
1989: PetscInt sum;
1991: for (i=0,sum=0;i<nb;i++) sum += blocks[i];
1992: if (sum == nvtxs) {
1993: PetscInt r;
1995: for (i=0,r=0;i<nb;i++) {
1996: PetscAssert(blocks[i] == xadj[r+1] - xadj[r],PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT,i,blocks[i],xadj[r+1] - xadj[r]);
1997: MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],sarray+xadj[r],ADD_VALUES);
1998: r += blocks[i];
1999: }
2000: } else {
2001: for (i=0;i<nvtxs;i++) {
2002: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
2003: }
2004: }
2005: } else {
2006: for (i=0;i<nvtxs;i++) {
2007: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
2008: }
2009: }
2010: MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2012: MatSeqAIJRestoreArray(local_mat,&sarray);
2013: } else { /* very basic values insertion for all other matrix types */
2014: PetscInt i;
2016: for (i=0;i<local_rows;i++) {
2017: PetscInt j;
2018: const PetscInt *local_indices_cols;
2020: MatGetRow(local_mat,i,&j,&local_indices_cols,&array);
2021: MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
2022: MatRestoreRow(local_mat,i,&j,&local_indices_cols,&array);
2023: }
2024: }
2025: MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
2026: MatDestroy(&local_mat);
2027: MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
2028: if (isseqdense) {
2029: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
2030: }
2031: if (reuse == MAT_INPLACE_MATRIX) {
2032: MatHeaderReplace(mat,&MT);
2033: } else if (reuse == MAT_INITIAL_MATRIX) {
2034: *M = MT;
2035: }
2036: return 0;
2037: }
2039: /*@
2040: MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format
2042: Input Parameters:
2043: + mat - the matrix (should be of type MATIS)
2044: - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
2046: Output Parameter:
2047: . newmat - the matrix in AIJ format
2049: Level: developer
2051: Notes:
2052: This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface.
2054: .seealso: MATIS, MatConvert()
2055: @*/
2056: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2057: {
2061: if (reuse == MAT_REUSE_MATRIX) {
2065: }
2066: PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
2067: return 0;
2068: }
2070: static PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
2071: {
2072: Mat_IS *matis = (Mat_IS*)(mat->data);
2073: PetscInt rbs,cbs,m,n,M,N;
2074: Mat B,localmat;
2076: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2077: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2078: MatGetSize(mat,&M,&N);
2079: MatGetLocalSize(mat,&m,&n);
2080: MatCreate(PetscObjectComm((PetscObject)mat),&B);
2081: MatSetSizes(B,m,n,M,N);
2082: MatSetBlockSize(B,rbs == cbs ? rbs : 1);
2083: MatSetType(B,MATIS);
2084: MatISSetLocalMatType(B,matis->lmattype);
2085: MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
2086: MatDuplicate(matis->A,op,&localmat);
2087: MatSetLocalToGlobalMapping(localmat,matis->A->rmap->mapping,matis->A->cmap->mapping);
2088: MatISSetLocalMat(B,localmat);
2089: MatDestroy(&localmat);
2090: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2091: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2092: *newmat = B;
2093: return 0;
2094: }
2096: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
2097: {
2098: Mat_IS *matis = (Mat_IS*)A->data;
2099: PetscBool local_sym;
2101: MatIsHermitian(matis->A,tol,&local_sym);
2102: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2103: return 0;
2104: }
2106: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
2107: {
2108: Mat_IS *matis = (Mat_IS*)A->data;
2109: PetscBool local_sym;
2111: if (matis->rmapping != matis->cmapping) {
2112: *flg = PETSC_FALSE;
2113: return 0;
2114: }
2115: MatIsSymmetric(matis->A,tol,&local_sym);
2116: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2117: return 0;
2118: }
2120: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg)
2121: {
2122: Mat_IS *matis = (Mat_IS*)A->data;
2123: PetscBool local_sym;
2125: if (matis->rmapping != matis->cmapping) {
2126: *flg = PETSC_FALSE;
2127: return 0;
2128: }
2129: MatIsStructurallySymmetric(matis->A,&local_sym);
2130: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2131: return 0;
2132: }
2134: static PetscErrorCode MatDestroy_IS(Mat A)
2135: {
2136: Mat_IS *b = (Mat_IS*)A->data;
2138: PetscFree(b->bdiag);
2139: PetscFree(b->lmattype);
2140: MatDestroy(&b->A);
2141: VecScatterDestroy(&b->cctx);
2142: VecScatterDestroy(&b->rctx);
2143: VecDestroy(&b->x);
2144: VecDestroy(&b->y);
2145: VecDestroy(&b->counter);
2146: ISDestroy(&b->getsub_ris);
2147: ISDestroy(&b->getsub_cis);
2148: if (b->sf != b->csf) {
2149: PetscSFDestroy(&b->csf);
2150: PetscFree2(b->csf_rootdata,b->csf_leafdata);
2151: } else b->csf = NULL;
2152: PetscSFDestroy(&b->sf);
2153: PetscFree2(b->sf_rootdata,b->sf_leafdata);
2154: ISLocalToGlobalMappingDestroy(&b->rmapping);
2155: ISLocalToGlobalMappingDestroy(&b->cmapping);
2156: PetscFree(A->data);
2157: PetscObjectChangeTypeName((PetscObject)A,NULL);
2158: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
2159: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
2160: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
2161: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
2162: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
2163: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
2164: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
2165: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",NULL);
2166: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
2167: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
2168: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
2169: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
2170: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
2171: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
2172: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
2173: PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",NULL);
2174: PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL);
2175: PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL);
2176: return 0;
2177: }
2179: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
2180: {
2181: Mat_IS *is = (Mat_IS*)A->data;
2182: PetscScalar zero = 0.0;
2184: /* scatter the global vector x into the local work vector */
2185: VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2186: VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2188: /* multiply the local matrix */
2189: MatMult(is->A,is->x,is->y);
2191: /* scatter product back into global memory */
2192: VecSet(y,zero);
2193: VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2194: VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2195: return 0;
2196: }
2198: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2199: {
2200: Vec temp_vec;
2202: if (v3 != v2) {
2203: MatMult(A,v1,v3);
2204: VecAXPY(v3,1.0,v2);
2205: } else {
2206: VecDuplicate(v2,&temp_vec);
2207: MatMult(A,v1,temp_vec);
2208: VecAXPY(temp_vec,1.0,v2);
2209: VecCopy(temp_vec,v3);
2210: VecDestroy(&temp_vec);
2211: }
2212: return 0;
2213: }
2215: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
2216: {
2217: Mat_IS *is = (Mat_IS*)A->data;
2219: /* scatter the global vector x into the local work vector */
2220: VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2221: VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2223: /* multiply the local matrix */
2224: MatMultTranspose(is->A,is->y,is->x);
2226: /* scatter product back into global vector */
2227: VecSet(x,0);
2228: VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2229: VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2230: return 0;
2231: }
2233: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2234: {
2235: Vec temp_vec;
2237: if (v3 != v2) {
2238: MatMultTranspose(A,v1,v3);
2239: VecAXPY(v3,1.0,v2);
2240: } else {
2241: VecDuplicate(v2,&temp_vec);
2242: MatMultTranspose(A,v1,temp_vec);
2243: VecAXPY(temp_vec,1.0,v2);
2244: VecCopy(temp_vec,v3);
2245: VecDestroy(&temp_vec);
2246: }
2247: return 0;
2248: }
2250: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
2251: {
2252: Mat_IS *a = (Mat_IS*)A->data;
2253: PetscViewer sviewer;
2254: PetscBool isascii,view = PETSC_TRUE;
2256: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2257: if (isascii) {
2258: PetscViewerFormat format;
2260: PetscViewerGetFormat(viewer,&format);
2261: if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2262: }
2263: if (!view) return 0;
2264: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2265: MatView(a->A,sviewer);
2266: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2267: PetscViewerFlush(viewer);
2268: return 0;
2269: }
2271: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values)
2272: {
2273: Mat_IS *is = (Mat_IS*)mat->data;
2274: MPI_Datatype nodeType;
2275: const PetscScalar *lv;
2276: PetscInt bs;
2278: MatGetBlockSize(mat,&bs);
2279: MatSetBlockSize(is->A,bs);
2280: MatInvertBlockDiagonal(is->A,&lv);
2281: if (!is->bdiag) {
2282: PetscMalloc1(bs*mat->rmap->n,&is->bdiag);
2283: }
2284: MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType);
2285: MPI_Type_commit(&nodeType);
2286: PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
2287: PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
2288: MPI_Type_free(&nodeType);
2289: if (values) *values = is->bdiag;
2290: return 0;
2291: }
2293: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2294: {
2295: Vec cglobal,rglobal;
2296: IS from;
2297: Mat_IS *is = (Mat_IS*)A->data;
2298: PetscScalar sum;
2299: const PetscInt *garray;
2300: PetscInt nr,rbs,nc,cbs;
2301: VecType rtype;
2303: ISLocalToGlobalMappingGetSize(is->rmapping,&nr);
2304: ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs);
2305: ISLocalToGlobalMappingGetSize(is->cmapping,&nc);
2306: ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs);
2307: VecDestroy(&is->x);
2308: VecDestroy(&is->y);
2309: VecDestroy(&is->counter);
2310: VecScatterDestroy(&is->rctx);
2311: VecScatterDestroy(&is->cctx);
2312: MatCreateVecs(is->A,&is->x,&is->y);
2313: VecBindToCPU(is->y,PETSC_TRUE);
2314: VecGetRootType_Private(is->y,&rtype);
2315: PetscFree(A->defaultvectype);
2316: PetscStrallocpy(rtype,&A->defaultvectype);
2317: MatCreateVecs(A,&cglobal,&rglobal);
2318: VecBindToCPU(rglobal,PETSC_TRUE);
2319: ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&garray);
2320: ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
2321: VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
2322: ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&garray);
2323: ISDestroy(&from);
2324: if (is->rmapping != is->cmapping) {
2325: ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&garray);
2326: ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
2327: VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
2328: ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&garray);
2329: ISDestroy(&from);
2330: } else {
2331: PetscObjectReference((PetscObject)is->rctx);
2332: is->cctx = is->rctx;
2333: }
2334: VecDestroy(&cglobal);
2336: /* interface counter vector (local) */
2337: VecDuplicate(is->y,&is->counter);
2338: VecBindToCPU(is->counter,PETSC_TRUE);
2339: VecSet(is->y,1.);
2340: VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2341: VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2342: VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2343: VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2344: VecBindToCPU(is->y,PETSC_FALSE);
2345: VecBindToCPU(is->counter,PETSC_FALSE);
2347: /* special functions for block-diagonal matrices */
2348: VecSum(rglobal,&sum);
2349: A->ops->invertblockdiagonal = NULL;
2350: if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2351: VecDestroy(&rglobal);
2353: /* setup SF for general purpose shared indices based communications */
2354: MatISSetUpSF_IS(A);
2355: return 0;
2356: }
2358: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2359: {
2360: IS is;
2361: ISLocalToGlobalMappingType l2gtype;
2362: const PetscInt *idxs;
2363: PetscHSetI ht;
2364: PetscInt *nidxs;
2365: PetscInt i,n,bs,c;
2366: PetscBool flg[] = {PETSC_FALSE,PETSC_FALSE};
2368: ISLocalToGlobalMappingGetSize(map,&n);
2369: ISLocalToGlobalMappingGetBlockSize(map,&bs);
2370: ISLocalToGlobalMappingGetBlockIndices(map,&idxs);
2371: PetscHSetICreate(&ht);
2372: PetscMalloc1(n/bs,&nidxs);
2373: for (i=0,c=0;i<n/bs;i++) {
2374: PetscBool missing;
2375: if (idxs[i] < 0) { flg[0] = PETSC_TRUE; continue; }
2376: PetscHSetIQueryAdd(ht,idxs[i],&missing);
2377: if (!missing) flg[1] = PETSC_TRUE;
2378: else nidxs[c++] = idxs[i];
2379: }
2380: PetscHSetIDestroy(&ht);
2381: MPIU_Allreduce(MPI_IN_PLACE,flg,2,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2382: if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2383: *nmap = NULL;
2384: *lmap = NULL;
2385: PetscFree(nidxs);
2386: ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs);
2387: return 0;
2388: }
2390: /* New l2g map without negative or repeated indices */
2391: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,c,nidxs,PETSC_USE_POINTER,&is);
2392: ISLocalToGlobalMappingCreateIS(is,nmap);
2393: ISDestroy(&is);
2394: ISLocalToGlobalMappingGetType(map,&l2gtype);
2395: ISLocalToGlobalMappingSetType(*nmap,l2gtype);
2397: /* New local l2g map for repeated indices */
2398: ISGlobalToLocalMappingApplyBlock(*nmap,IS_GTOLM_MASK,n/bs,idxs,NULL,nidxs);
2399: ISCreateBlock(PETSC_COMM_SELF,bs,n/bs,nidxs,PETSC_USE_POINTER,&is);
2400: ISLocalToGlobalMappingCreateIS(is,lmap);
2401: ISDestroy(&is);
2403: PetscFree(nidxs);
2404: ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs);
2405: return 0;
2406: }
2408: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
2409: {
2410: Mat_IS *is = (Mat_IS*)A->data;
2411: ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2412: PetscBool cong, freem[] = { PETSC_FALSE, PETSC_FALSE };
2413: PetscInt nr,rbs,nc,cbs;
2418: ISLocalToGlobalMappingDestroy(&is->rmapping);
2419: ISLocalToGlobalMappingDestroy(&is->cmapping);
2420: PetscLayoutSetUp(A->rmap);
2421: PetscLayoutSetUp(A->cmap);
2422: MatHasCongruentLayouts(A,&cong);
2424: /* If NULL, local space matches global space */
2425: if (!rmapping) {
2426: IS is;
2428: ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->N,0,1,&is);
2429: ISLocalToGlobalMappingCreateIS(is,&rmapping);
2430: if (A->rmap->bs > 0) ISLocalToGlobalMappingSetBlockSize(rmapping,A->rmap->bs);
2431: ISDestroy(&is);
2432: freem[0] = PETSC_TRUE;
2433: if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2434: } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2435: MatISFilterL2GMap(A,rmapping,&is->rmapping,&localrmapping);
2436: if (rmapping == cmapping) {
2437: PetscObjectReference((PetscObject)is->rmapping);
2438: is->cmapping = is->rmapping;
2439: PetscObjectReference((PetscObject)localrmapping);
2440: localcmapping = localrmapping;
2441: }
2442: }
2443: if (!cmapping) {
2444: IS is;
2446: ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->N,0,1,&is);
2447: ISLocalToGlobalMappingCreateIS(is,&cmapping);
2448: if (A->cmap->bs > 0) ISLocalToGlobalMappingSetBlockSize(cmapping,A->cmap->bs);
2449: ISDestroy(&is);
2450: freem[1] = PETSC_TRUE;
2451: } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2452: MatISFilterL2GMap(A,cmapping,&is->cmapping,&localcmapping);
2453: }
2454: if (!is->rmapping) {
2455: PetscObjectReference((PetscObject)rmapping);
2456: is->rmapping = rmapping;
2457: }
2458: if (!is->cmapping) {
2459: PetscObjectReference((PetscObject)cmapping);
2460: is->cmapping = cmapping;
2461: }
2463: /* Clean up */
2464: MatDestroy(&is->A);
2465: if (is->csf != is->sf) {
2466: PetscSFDestroy(&is->csf);
2467: PetscFree2(is->csf_rootdata,is->csf_leafdata);
2468: } else is->csf = NULL;
2469: PetscSFDestroy(&is->sf);
2470: PetscFree2(is->sf_rootdata,is->sf_leafdata);
2471: PetscFree(is->bdiag);
2473: /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2474: (DOLFIN passes 2 different objects) */
2475: ISLocalToGlobalMappingGetSize(is->rmapping,&nr);
2476: ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs);
2477: ISLocalToGlobalMappingGetSize(is->cmapping,&nc);
2478: ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs);
2479: if (is->rmapping != is->cmapping && cong) {
2480: PetscBool same = PETSC_FALSE;
2481: if (nr == nc && cbs == rbs) {
2482: const PetscInt *idxs1,*idxs2;
2484: ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&idxs1);
2485: ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&idxs2);
2486: PetscArraycmp(idxs1,idxs2,nr/rbs,&same);
2487: ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&idxs1);
2488: ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&idxs2);
2489: }
2490: MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2491: if (same) {
2492: ISLocalToGlobalMappingDestroy(&is->cmapping);
2493: PetscObjectReference((PetscObject)is->rmapping);
2494: is->cmapping = is->rmapping;
2495: }
2496: }
2497: PetscLayoutSetBlockSize(A->rmap,rbs);
2498: PetscLayoutSetBlockSize(A->cmap,cbs);
2499: /* Pass the user defined maps to the layout */
2500: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
2501: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
2502: if (freem[0]) ISLocalToGlobalMappingDestroy(&rmapping);
2503: if (freem[1]) ISLocalToGlobalMappingDestroy(&cmapping);
2505: /* Create the local matrix A */
2506: MatCreate(PETSC_COMM_SELF,&is->A);
2507: MatSetType(is->A,is->lmattype);
2508: MatSetSizes(is->A,nr,nc,nr,nc);
2509: MatSetBlockSizes(is->A,rbs,cbs);
2510: MatSetOptionsPrefix(is->A,"is_");
2511: MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);
2512: PetscLayoutSetUp(is->A->rmap);
2513: PetscLayoutSetUp(is->A->cmap);
2514: MatSetLocalToGlobalMapping(is->A,localrmapping,localcmapping);
2515: ISLocalToGlobalMappingDestroy(&localrmapping);
2516: ISLocalToGlobalMappingDestroy(&localcmapping);
2518: /* setup scatters and local vectors for MatMult */
2519: if (!is->islocalref) MatISSetUpScatters_Private(A);
2520: A->preallocated = PETSC_TRUE;
2521: return 0;
2522: }
2524: static PetscErrorCode MatSetUp_IS(Mat A)
2525: {
2526: ISLocalToGlobalMapping rmap, cmap;
2528: MatGetLocalToGlobalMapping(A,&rmap,&cmap);
2529: if (!rmap && !cmap) {
2530: MatSetLocalToGlobalMapping(A,NULL,NULL);
2531: }
2532: return 0;
2533: }
2535: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2536: {
2537: Mat_IS *is = (Mat_IS*)mat->data;
2538: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2540: ISGlobalToLocalMappingApply(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2541: if (m != n || rows != cols || is->cmapping != is->rmapping) {
2542: ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2543: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
2544: } else {
2545: MatSetValues(is->A,m,rows_l,m,rows_l,values,addv);
2546: }
2547: return 0;
2548: }
2550: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2551: {
2552: Mat_IS *is = (Mat_IS*)mat->data;
2553: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2555: ISGlobalToLocalMappingApplyBlock(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2556: if (m != n || rows != cols || is->cmapping != is->rmapping) {
2557: ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2558: MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
2559: } else {
2560: MatSetValuesBlocked(is->A,m,rows_l,n,rows_l,values,addv);
2561: }
2562: return 0;
2563: }
2565: static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2566: {
2567: Mat_IS *is = (Mat_IS*)A->data;
2569: if (is->A->rmap->mapping || is->A->cmap->mapping) {
2570: MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
2571: } else {
2572: MatSetValues(is->A,m,rows,n,cols,values,addv);
2573: }
2574: return 0;
2575: }
2577: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2578: {
2579: Mat_IS *is = (Mat_IS*)A->data;
2581: if (is->A->rmap->mapping || is->A->cmap->mapping) {
2582: MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
2583: } else {
2584: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
2585: }
2586: return 0;
2587: }
2589: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns)
2590: {
2591: Mat_IS *is = (Mat_IS*)A->data;
2593: if (!n) {
2594: is->pure_neumann = PETSC_TRUE;
2595: } else {
2596: PetscInt i;
2597: is->pure_neumann = PETSC_FALSE;
2599: if (columns) {
2600: MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL);
2601: } else {
2602: MatZeroRows(is->A,n,rows,diag,NULL,NULL);
2603: }
2604: if (diag != 0.) {
2605: const PetscScalar *array;
2606: VecGetArrayRead(is->counter,&array);
2607: for (i=0; i<n; i++) {
2608: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
2609: }
2610: VecRestoreArrayRead(is->counter,&array);
2611: }
2612: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
2613: MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
2614: }
2615: return 0;
2616: }
2618: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
2619: {
2620: Mat_IS *matis = (Mat_IS*)A->data;
2621: PetscInt nr,nl,len,i;
2622: PetscInt *lrows;
2624: if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2625: PetscBool cong;
2627: PetscLayoutCompare(A->rmap,A->cmap,&cong);
2628: cong = (PetscBool)(cong && matis->sf == matis->csf);
2632: }
2633: /* get locally owned rows */
2634: PetscLayoutMapLocal(A->rmap,n,rows,&len,&lrows,NULL);
2635: /* fix right hand side if needed */
2636: if (x && b) {
2637: const PetscScalar *xx;
2638: PetscScalar *bb;
2640: VecGetArrayRead(x, &xx);
2641: VecGetArray(b, &bb);
2642: for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
2643: VecRestoreArrayRead(x, &xx);
2644: VecRestoreArray(b, &bb);
2645: }
2646: /* get rows associated to the local matrices */
2647: MatGetSize(matis->A,&nl,NULL);
2648: PetscArrayzero(matis->sf_leafdata,nl);
2649: PetscArrayzero(matis->sf_rootdata,A->rmap->n);
2650: for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
2651: PetscFree(lrows);
2652: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
2653: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
2654: PetscMalloc1(nl,&lrows);
2655: for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
2656: MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
2657: PetscFree(lrows);
2658: return 0;
2659: }
2661: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2662: {
2663: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
2664: return 0;
2665: }
2667: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2668: {
2669: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
2670: return 0;
2671: }
2673: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
2674: {
2675: Mat_IS *is = (Mat_IS*)A->data;
2677: MatAssemblyBegin(is->A,type);
2678: return 0;
2679: }
2681: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
2682: {
2683: Mat_IS *is = (Mat_IS*)A->data;
2685: MatAssemblyEnd(is->A,type);
2686: /* fix for local empty rows/cols */
2687: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2688: Mat newlA;
2689: ISLocalToGlobalMapping rl2g,cl2g;
2690: IS nzr,nzc;
2691: PetscInt nr,nc,nnzr,nnzc;
2692: PetscBool lnewl2g,newl2g;
2694: MatGetSize(is->A,&nr,&nc);
2695: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
2696: if (!nzr) {
2697: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
2698: }
2699: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
2700: if (!nzc) {
2701: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
2702: }
2703: ISGetSize(nzr,&nnzr);
2704: ISGetSize(nzc,&nnzc);
2705: if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
2706: lnewl2g = PETSC_TRUE;
2707: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2709: /* extract valid submatrix */
2710: MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);
2711: } else { /* local matrix fully populated */
2712: lnewl2g = PETSC_FALSE;
2713: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2714: PetscObjectReference((PetscObject)is->A);
2715: newlA = is->A;
2716: }
2718: /* attach new global l2g map if needed */
2719: if (newl2g) {
2720: IS zr,zc;
2721: const PetscInt *ridxs,*cidxs,*zridxs,*zcidxs;
2722: PetscInt *nidxs,i;
2724: ISComplement(nzr,0,nr,&zr);
2725: ISComplement(nzc,0,nc,&zc);
2726: PetscMalloc1(PetscMax(nr,nc),&nidxs);
2727: ISLocalToGlobalMappingGetIndices(is->rmapping,&ridxs);
2728: ISLocalToGlobalMappingGetIndices(is->cmapping,&cidxs);
2729: ISGetIndices(zr,&zridxs);
2730: ISGetIndices(zc,&zcidxs);
2731: ISGetLocalSize(zr,&nnzr);
2732: ISGetLocalSize(zc,&nnzc);
2734: PetscArraycpy(nidxs,ridxs,nr);
2735: for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
2736: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nr,nidxs,PETSC_COPY_VALUES,&rl2g);
2737: PetscArraycpy(nidxs,cidxs,nc);
2738: for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
2739: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nc,nidxs,PETSC_COPY_VALUES,&cl2g);
2741: ISRestoreIndices(zr,&zridxs);
2742: ISRestoreIndices(zc,&zcidxs);
2743: ISLocalToGlobalMappingRestoreIndices(is->rmapping,&ridxs);
2744: ISLocalToGlobalMappingRestoreIndices(is->cmapping,&cidxs);
2745: ISDestroy(&nzr);
2746: ISDestroy(&nzc);
2747: ISDestroy(&zr);
2748: ISDestroy(&zc);
2749: PetscFree(nidxs);
2750: MatSetLocalToGlobalMapping(A,rl2g,cl2g);
2751: ISLocalToGlobalMappingDestroy(&rl2g);
2752: ISLocalToGlobalMappingDestroy(&cl2g);
2753: }
2754: MatISSetLocalMat(A,newlA);
2755: MatDestroy(&newlA);
2756: ISDestroy(&nzr);
2757: ISDestroy(&nzc);
2758: is->locempty = PETSC_FALSE;
2759: }
2760: return 0;
2761: }
2763: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
2764: {
2765: Mat_IS *is = (Mat_IS*)mat->data;
2767: *local = is->A;
2768: return 0;
2769: }
2771: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
2772: {
2773: *local = NULL;
2774: return 0;
2775: }
2777: /*@
2778: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
2780: Input Parameter:
2781: . mat - the matrix
2783: Output Parameter:
2784: . local - the local matrix
2786: Level: advanced
2788: Notes:
2789: This can be called if you have precomputed the nonzero structure of the
2790: matrix and want to provide it to the inner matrix object to improve the performance
2791: of the MatSetValues() operation.
2793: Call MatISRestoreLocalMat() when finished with the local matrix.
2795: .seealso: MATIS
2796: @*/
2797: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
2798: {
2801: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
2802: return 0;
2803: }
2805: /*@
2806: MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()
2808: Input Parameter:
2809: . mat - the matrix
2811: Output Parameter:
2812: . local - the local matrix
2814: Level: advanced
2816: .seealso: MATIS
2817: @*/
2818: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
2819: {
2822: PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
2823: return 0;
2824: }
2826: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
2827: {
2828: Mat_IS *is = (Mat_IS*)mat->data;
2830: if (is->A) {
2831: MatSetType(is->A,mtype);
2832: }
2833: PetscFree(is->lmattype);
2834: PetscStrallocpy(mtype,&is->lmattype);
2835: return 0;
2836: }
2838: /*@
2839: MatISSetLocalMatType - Specifies the type of local matrix
2841: Input Parameters:
2842: + mat - the matrix
2843: - mtype - the local matrix type
2845: Output Parameter:
2847: Level: advanced
2849: .seealso: MATIS, MatSetType(), MatType
2850: @*/
2851: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
2852: {
2854: PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
2855: return 0;
2856: }
2858: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
2859: {
2860: Mat_IS *is = (Mat_IS*)mat->data;
2861: PetscInt nrows,ncols,orows,ocols;
2862: MatType mtype,otype;
2863: PetscBool sametype = PETSC_TRUE;
2865: if (is->A && !is->islocalref) {
2866: MatGetSize(is->A,&orows,&ocols);
2867: MatGetSize(local,&nrows,&ncols);
2869: MatGetType(local,&mtype);
2870: MatGetType(is->A,&otype);
2871: PetscStrcmp(mtype,otype,&sametype);
2872: }
2873: PetscObjectReference((PetscObject)local);
2874: MatDestroy(&is->A);
2875: is->A = local;
2876: MatGetType(is->A,&mtype);
2877: MatISSetLocalMatType(mat,mtype);
2878: if (!sametype && !is->islocalref) {
2879: MatISSetUpScatters_Private(mat);
2880: }
2881: return 0;
2882: }
2884: /*@
2885: MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.
2887: Collective on Mat
2889: Input Parameters:
2890: + mat - the matrix
2891: - local - the local matrix
2893: Output Parameter:
2895: Level: advanced
2897: Notes:
2898: This can be called if you have precomputed the local matrix and
2899: want to provide it to the matrix object MATIS.
2901: .seealso: MATIS
2902: @*/
2903: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
2904: {
2907: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
2908: return 0;
2909: }
2911: static PetscErrorCode MatZeroEntries_IS(Mat A)
2912: {
2913: Mat_IS *a = (Mat_IS*)A->data;
2915: MatZeroEntries(a->A);
2916: return 0;
2917: }
2919: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
2920: {
2921: Mat_IS *is = (Mat_IS*)A->data;
2923: MatScale(is->A,a);
2924: return 0;
2925: }
2927: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2928: {
2929: Mat_IS *is = (Mat_IS*)A->data;
2931: /* get diagonal of the local matrix */
2932: MatGetDiagonal(is->A,is->y);
2934: /* scatter diagonal back into global vector */
2935: VecSet(v,0);
2936: VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2937: VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2938: return 0;
2939: }
2941: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
2942: {
2943: Mat_IS *a = (Mat_IS*)A->data;
2945: MatSetOption(a->A,op,flg);
2946: return 0;
2947: }
2949: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
2950: {
2951: Mat_IS *y = (Mat_IS*)Y->data;
2952: Mat_IS *x;
2954: if (PetscDefined(USE_DEBUG)) {
2955: PetscBool ismatis;
2956: PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
2958: }
2959: x = (Mat_IS*)X->data;
2960: MatAXPY(y->A,a,x->A,str);
2961: return 0;
2962: }
2964: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
2965: {
2966: Mat lA;
2967: Mat_IS *matis = (Mat_IS*)(A->data);
2968: ISLocalToGlobalMapping rl2g,cl2g;
2969: IS is;
2970: const PetscInt *rg,*rl;
2971: PetscInt nrg;
2972: PetscInt N,M,nrl,i,*idxs;
2974: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
2975: ISGetLocalSize(row,&nrl);
2976: ISGetIndices(row,&rl);
2977: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
2978: if (PetscDefined(USE_DEBUG)) {
2980: }
2981: PetscMalloc1(nrg,&idxs);
2982: /* map from [0,nrl) to row */
2983: for (i=0;i<nrl;i++) idxs[i] = rl[i];
2984: for (i=nrl;i<nrg;i++) idxs[i] = -1;
2985: ISRestoreIndices(row,&rl);
2986: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
2987: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
2988: ISLocalToGlobalMappingCreateIS(is,&rl2g);
2989: ISDestroy(&is);
2990: /* compute new l2g map for columns */
2991: if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
2992: const PetscInt *cg,*cl;
2993: PetscInt ncg;
2994: PetscInt ncl;
2996: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
2997: ISGetLocalSize(col,&ncl);
2998: ISGetIndices(col,&cl);
2999: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
3000: if (PetscDefined(USE_DEBUG)) {
3002: }
3003: PetscMalloc1(ncg,&idxs);
3004: /* map from [0,ncl) to col */
3005: for (i=0;i<ncl;i++) idxs[i] = cl[i];
3006: for (i=ncl;i<ncg;i++) idxs[i] = -1;
3007: ISRestoreIndices(col,&cl);
3008: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
3009: ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
3010: ISLocalToGlobalMappingCreateIS(is,&cl2g);
3011: ISDestroy(&is);
3012: } else {
3013: PetscObjectReference((PetscObject)rl2g);
3014: cl2g = rl2g;
3015: }
3016: /* create the MATIS submatrix */
3017: MatGetSize(A,&M,&N);
3018: MatCreate(PetscObjectComm((PetscObject)A),submat);
3019: MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
3020: MatSetType(*submat,MATIS);
3021: matis = (Mat_IS*)((*submat)->data);
3022: matis->islocalref = PETSC_TRUE;
3023: MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
3024: MatISGetLocalMat(A,&lA);
3025: MatISSetLocalMat(*submat,lA);
3026: MatSetUp(*submat);
3027: MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
3028: MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
3029: ISLocalToGlobalMappingDestroy(&rl2g);
3030: ISLocalToGlobalMappingDestroy(&cl2g);
3032: /* remove unsupported ops */
3033: PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
3034: (*submat)->ops->destroy = MatDestroy_IS;
3035: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3036: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3037: (*submat)->ops->assemblybegin = MatAssemblyBegin_IS;
3038: (*submat)->ops->assemblyend = MatAssemblyEnd_IS;
3039: return 0;
3040: }
3042: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
3043: {
3044: Mat_IS *a = (Mat_IS*)A->data;
3045: char type[256];
3046: PetscBool flg;
3048: PetscOptionsHead(PetscOptionsObject,"MATIS options");
3049: PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
3050: PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
3051: PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
3052: if (flg) {
3053: MatISSetLocalMatType(A,type);
3054: }
3055: if (a->A) {
3056: MatSetFromOptions(a->A);
3057: }
3058: PetscOptionsTail();
3059: return 0;
3060: }
3062: /*@
3063: MatCreateIS - Creates a "process" unassembled matrix, assembled on each
3064: process but not across processes.
3066: Input Parameters:
3067: + comm - MPI communicator that will share the matrix
3068: . bs - block size of the matrix
3069: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3070: . rmap - local to global map for rows
3071: - cmap - local to global map for cols
3073: Output Parameter:
3074: . A - the resulting matrix
3076: Level: advanced
3078: Notes:
3079: See MATIS for more details.
3080: m and n are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors
3081: used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
3082: If rmap (cmap) is NULL, then the local row (column) spaces matches the global space.
3084: .seealso: MATIS, MatSetLocalToGlobalMapping()
3085: @*/
3086: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
3087: {
3088: MatCreate(comm,A);
3089: MatSetSizes(*A,m,n,M,N);
3090: if (bs > 0) {
3091: MatSetBlockSize(*A,bs);
3092: }
3093: MatSetType(*A,MATIS);
3094: MatSetLocalToGlobalMapping(*A,rmap,cmap);
3095: return 0;
3096: }
3098: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3099: {
3100: Mat_IS *a = (Mat_IS*)A->data;
3102: *has = PETSC_FALSE;
3103: if (!((void**)A->ops)[op]) return 0;
3104: MatHasOperation(a->A,op,has);
3105: return 0;
3106: }
3108: static PetscErrorCode MatSetValuesCOO_IS(Mat A,const PetscScalar v[],InsertMode imode)
3109: {
3110: Mat_IS *a = (Mat_IS*)A->data;
3112: MatSetValuesCOO(a->A,v,imode);
3113: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
3114: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
3115: return 0;
3116: }
3118: static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A,PetscCount ncoo,PetscInt coo_i[],PetscInt coo_j[])
3119: {
3120: Mat_IS *a = (Mat_IS*)A->data;
3123: if (a->A->rmap->mapping || a->A->cmap->mapping) {
3124: MatSetPreallocationCOOLocal(a->A,ncoo,coo_i,coo_j);
3125: } else {
3126: MatSetPreallocationCOO(a->A,ncoo,coo_i,coo_j);
3127: }
3128: PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS);
3129: A->preallocated = PETSC_TRUE;
3130: return 0;
3131: }
3133: static PetscErrorCode MatSetPreallocationCOO_IS(Mat A,PetscCount ncoo,const PetscInt coo_i[],const PetscInt coo_j[])
3134: {
3135: Mat_IS *a = (Mat_IS*)A->data;
3136: PetscInt *coo_il, *coo_jl, incoo;
3140: PetscMalloc2(ncoo,&coo_il,ncoo,&coo_jl);
3141: ISGlobalToLocalMappingApply(a->rmapping,IS_GTOLM_MASK,ncoo,coo_i,&incoo,coo_il);
3142: ISGlobalToLocalMappingApply(a->cmapping,IS_GTOLM_MASK,ncoo,coo_j,&incoo,coo_jl);
3143: MatSetPreallocationCOO(a->A,ncoo,coo_il,coo_jl);
3144: PetscFree2(coo_il,coo_jl);
3145: PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS);
3146: A->preallocated = PETSC_TRUE;
3147: return 0;
3148: }
3150: /*@
3151: MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the MATIS object
3153: Not Collective
3155: Input Parameter:
3156: . A - the matrix
3158: Output Parameters:
3159: + rmapping - row mapping
3160: - cmapping - column mapping
3162: Notes: The returned map can be different from the one used to construct the MATIS object, since it will not contain negative or repeated indices.
3164: Level: advanced
3166: .seealso: MatSetLocalToGlobalMapping()
3167: @*/
3168: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping)
3169: {
3174: PetscUseMethod(A,"MatISGetLocalToGlobalMapping_C",(Mat,ISLocalToGlobalMapping*,ISLocalToGlobalMapping*),(A,rmapping,cmapping));
3175: return 0;
3176: }
3178: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3179: {
3180: Mat_IS *a = (Mat_IS*)A->data;
3182: if (r) *r = a->rmapping;
3183: if (c) *c = a->cmapping;
3184: return 0;
3185: }
3187: /*MC
3188: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP).
3189: This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector
3190: product is handled "implicitly".
3192: Options Database Keys:
3193: + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
3194: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3195: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP().
3197: Notes:
3198: Options prefix for the inner matrix are given by -is_mat_xxx
3200: You must call MatSetLocalToGlobalMapping() before using this matrix type.
3202: You can do matrix preallocation on the local matrix after you obtain it with
3203: MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation()
3205: Level: advanced
3207: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP
3209: M*/
3210: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3211: {
3212: Mat_IS *a;
3214: PetscNewLog(A,&a);
3215: PetscStrallocpy(MATAIJ,&a->lmattype);
3216: A->data = (void*)a;
3218: /* matrix ops */
3219: PetscMemzero(A->ops,sizeof(struct _MatOps));
3220: A->ops->mult = MatMult_IS;
3221: A->ops->multadd = MatMultAdd_IS;
3222: A->ops->multtranspose = MatMultTranspose_IS;
3223: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3224: A->ops->destroy = MatDestroy_IS;
3225: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3226: A->ops->setvalues = MatSetValues_IS;
3227: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3228: A->ops->setvalueslocal = MatSetValuesLocal_IS;
3229: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3230: A->ops->zerorows = MatZeroRows_IS;
3231: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3232: A->ops->assemblybegin = MatAssemblyBegin_IS;
3233: A->ops->assemblyend = MatAssemblyEnd_IS;
3234: A->ops->view = MatView_IS;
3235: A->ops->zeroentries = MatZeroEntries_IS;
3236: A->ops->scale = MatScale_IS;
3237: A->ops->getdiagonal = MatGetDiagonal_IS;
3238: A->ops->setoption = MatSetOption_IS;
3239: A->ops->ishermitian = MatIsHermitian_IS;
3240: A->ops->issymmetric = MatIsSymmetric_IS;
3241: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3242: A->ops->duplicate = MatDuplicate_IS;
3243: A->ops->missingdiagonal = MatMissingDiagonal_IS;
3244: A->ops->copy = MatCopy_IS;
3245: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3246: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3247: A->ops->axpy = MatAXPY_IS;
3248: A->ops->diagonalset = MatDiagonalSet_IS;
3249: A->ops->shift = MatShift_IS;
3250: A->ops->transpose = MatTranspose_IS;
3251: A->ops->getinfo = MatGetInfo_IS;
3252: A->ops->diagonalscale = MatDiagonalScale_IS;
3253: A->ops->setfromoptions = MatSetFromOptions_IS;
3254: A->ops->setup = MatSetUp_IS;
3255: A->ops->hasoperation = MatHasOperation_IS;
3257: /* special MATIS functions */
3258: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
3259: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
3260: PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
3261: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
3262: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
3263: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
3264: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
3265: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
3266: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",MatISGetLocalToGlobalMapping_IS);
3267: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
3268: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
3269: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
3270: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
3271: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
3272: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
3273: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
3274: PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",MatSetPreallocationCOOLocal_IS);
3275: PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_IS);
3276: PetscObjectChangeTypeName((PetscObject)A,MATIS);
3277: return 0;
3278: }