Actual source code: matis.c
petsc-3.11.4 2019-09-28
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>
13: #define MATIS_MAX_ENTRIES_INSERTION 2048
14: static PetscErrorCode MatSetValuesLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode);
15: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode);
16: static PetscErrorCode MatISSetUpScatters_Private(Mat);
18: static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr)
19: {
20: 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;
45: PetscObjectQuery((PetscObject)C,"_MatIS_PtAP",(PetscObject*)&c);
46: if (!c) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing PtAP information");
47: PetscContainerGetPointer(c,(void**)&ptap);
48: ris[0] = ptap->ris0;
49: ris[1] = ptap->ris1;
50: cis[0] = ptap->cis0;
51: cis[1] = ptap->cis1;
52: n = ptap->ris1 ? 2 : 1;
53: reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
54: MatCreateSubMatrices(P,n,ris,cis,reuse,&ptap->lP);
56: MatISGetLocalMat(A,&lA);
57: MatISGetLocalMat(C,&lC);
58: if (ptap->ris1) { /* unsymmetric A mapping */
59: Mat lPt;
61: MatTranspose(ptap->lP[1],MAT_INITIAL_MATRIX,&lPt);
62: MatMatMatMult(lPt,lA,ptap->lP[0],reuse,ptap->fill,&lC);
63: if (matis->storel2l) {
64: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",(PetscObject)lPt);
65: }
66: MatDestroy(&lPt);
67: } else {
68: MatPtAP(lA,ptap->lP[0],reuse,ptap->fill,&lC);
69: if (matis->storel2l) {
70: PetscObjectCompose((PetscObject)(C),"_MatIS_PtAP_l2l",(PetscObject)ptap->lP[0]);
71: }
72: }
73: if (reuse == MAT_INITIAL_MATRIX) {
74: MatISSetLocalMat(C,lC);
75: MatDestroy(&lC);
76: }
77: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
79: return(0);
80: }
82: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT,IS *cis)
83: {
84: Mat Po,Pd;
85: IS zd,zo;
86: const PetscInt *garray;
87: PetscInt *aux,i,bs;
88: PetscInt dc,stc,oc,ctd,cto;
89: PetscBool ismpiaij,ismpibaij,isseqaij,isseqbaij;
90: MPI_Comm comm;
96: PetscObjectGetComm((PetscObject)PT,&comm);
97: bs = 1;
98: PetscObjectTypeCompare((PetscObject)PT,MATMPIAIJ,&ismpiaij);
99: PetscObjectTypeCompare((PetscObject)PT,MATMPIBAIJ,&ismpibaij);
100: PetscObjectBaseTypeCompare((PetscObject)PT,MATSEQAIJ,&isseqaij);
101: PetscObjectTypeCompare((PetscObject)PT,MATSEQBAIJ,&isseqbaij);
102: if (isseqaij || isseqbaij) {
103: Pd = PT;
104: Po = NULL;
105: garray = NULL;
106: } else if (ismpiaij) {
107: MatMPIAIJGetSeqAIJ(PT,&Pd,&Po,&garray);
108: } else if (ismpibaij) {
109: MatMPIBAIJGetSeqBAIJ(PT,&Pd,&Po,&garray);
110: MatGetBlockSize(PT,&bs);
111: } else SETERRQ1(comm,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(PT))->type_name);
113: /* identify any null columns in Pd or Po */
114: /* We use a tolerance comparison since it may happen that, with geometric multigrid,
115: some of the columns are not really zero, but very close to */
116: zo = zd = NULL;
117: if (Po) {
118: MatFindNonzeroRowsOrCols_Basic(Po,PETSC_TRUE,PETSC_SMALL,&zo);
119: }
120: MatFindNonzeroRowsOrCols_Basic(Pd,PETSC_TRUE,PETSC_SMALL,&zd);
122: MatGetLocalSize(PT,NULL,&dc);
123: MatGetOwnershipRangeColumn(PT,&stc,NULL);
124: if (Po) { MatGetLocalSize(Po,NULL,&oc); }
125: else oc = 0;
126: PetscMalloc1((dc+oc)/bs,&aux);
127: if (zd) {
128: const PetscInt *idxs;
129: PetscInt nz;
131: /* this will throw an error if bs is not valid */
132: ISSetBlockSize(zd,bs);
133: ISGetLocalSize(zd,&nz);
134: ISGetIndices(zd,&idxs);
135: ctd = nz/bs;
136: for (i=0; i<ctd; i++) aux[i] = (idxs[bs*i]+stc)/bs;
137: ISRestoreIndices(zd,&idxs);
138: } else {
139: ctd = dc/bs;
140: for (i=0; i<ctd; i++) aux[i] = i+stc/bs;
141: }
142: if (zo) {
143: const PetscInt *idxs;
144: PetscInt nz;
146: /* this will throw an error if bs is not valid */
147: ISSetBlockSize(zo,bs);
148: ISGetLocalSize(zo,&nz);
149: ISGetIndices(zo,&idxs);
150: cto = nz/bs;
151: for (i=0; i<cto; i++) aux[i+ctd] = garray[idxs[bs*i]/bs];
152: ISRestoreIndices(zo,&idxs);
153: } else {
154: cto = oc/bs;
155: for (i=0; i<cto; i++) aux[i+ctd] = garray[i];
156: }
157: ISCreateBlock(comm,bs,ctd+cto,aux,PETSC_OWN_POINTER,cis);
158: ISDestroy(&zd);
159: ISDestroy(&zo);
160: return(0);
161: }
163: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
164: {
165: Mat PT,lA;
166: MatISPtAP ptap;
167: ISLocalToGlobalMapping Crl2g,Ccl2g,rl2g,cl2g;
168: PetscContainer c;
169: MatType lmtype;
170: const PetscInt *garray;
171: PetscInt ibs,N,dc;
172: MPI_Comm comm;
173: PetscErrorCode ierr;
176: PetscObjectGetComm((PetscObject)A,&comm);
177: MatCreate(comm,C);
178: MatSetType(*C,MATIS);
179: MatISGetLocalMat(A,&lA);
180: MatGetType(lA,&lmtype);
181: MatISSetLocalMatType(*C,lmtype);
182: MatGetSize(P,NULL,&N);
183: MatGetLocalSize(P,NULL,&dc);
184: MatSetSizes(*C,dc,dc,N,N);
185: /* Not sure about this
186: MatGetBlockSizes(P,NULL,&ibs);
187: MatSetBlockSize(*C,ibs);
188: */
190: PetscNew(&ptap);
191: PetscContainerCreate(PETSC_COMM_SELF,&c);
192: PetscContainerSetPointer(c,ptap);
193: PetscContainerSetUserDestroy(c,MatISContainerDestroyPtAP_Private);
194: PetscObjectCompose((PetscObject)(*C),"_MatIS_PtAP",(PetscObject)c);
195: PetscContainerDestroy(&c);
196: ptap->fill = fill;
198: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
200: ISLocalToGlobalMappingGetBlockSize(cl2g,&ibs);
201: ISLocalToGlobalMappingGetSize(cl2g,&N);
202: ISLocalToGlobalMappingGetBlockIndices(cl2g,&garray);
203: ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris0);
204: ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&garray);
206: MatCreateSubMatrix(P,ptap->ris0,NULL,MAT_INITIAL_MATRIX,&PT);
207: MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis0);
208: ISLocalToGlobalMappingCreateIS(ptap->cis0,&Ccl2g);
209: MatDestroy(&PT);
211: Crl2g = NULL;
212: if (rl2g != cl2g) { /* unsymmetric A mapping */
213: PetscBool same,lsame = PETSC_FALSE;
214: PetscInt N1,ibs1;
216: ISLocalToGlobalMappingGetSize(rl2g,&N1);
217: ISLocalToGlobalMappingGetBlockSize(rl2g,&ibs1);
218: ISLocalToGlobalMappingGetBlockIndices(rl2g,&garray);
219: ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris1);
220: ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&garray);
221: if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
222: const PetscInt *i1,*i2;
224: ISBlockGetIndices(ptap->ris0,&i1);
225: ISBlockGetIndices(ptap->ris1,&i2);
226: PetscMemcmp(i1,i2,N*sizeof(*i1),&lsame);
227: }
228: MPIU_Allreduce(&lsame,&same,1,MPIU_BOOL,MPI_LAND,comm);
229: if (same) {
230: ISDestroy(&ptap->ris1);
231: } else {
232: MatCreateSubMatrix(P,ptap->ris1,NULL,MAT_INITIAL_MATRIX,&PT);
233: MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis1);
234: ISLocalToGlobalMappingCreateIS(ptap->cis1,&Crl2g);
235: MatDestroy(&PT);
236: }
237: }
238: /* Not sure about this
239: if (!Crl2g) {
240: MatGetBlockSize(*C,&ibs);
241: ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs);
242: }
243: */
244: MatSetLocalToGlobalMapping(*C,Crl2g ? Crl2g : Ccl2g,Ccl2g);
245: ISLocalToGlobalMappingDestroy(&Crl2g);
246: ISLocalToGlobalMappingDestroy(&Ccl2g);
248: (*C)->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
249: return(0);
250: }
252: PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
253: {
257: if (scall == MAT_INITIAL_MATRIX) {
258: PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
259: MatPtAPSymbolic_IS_XAIJ(A,P,fill,C);
260: PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
261: }
262: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
263: ((*C)->ops->ptapnumeric)(A,P,*C);
264: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
265: return(0);
266: }
268: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
269: {
270: MatISLocalFields lf = (MatISLocalFields)ptr;
271: PetscInt i;
272: PetscErrorCode ierr;
275: for (i=0;i<lf->nr;i++) {
276: ISDestroy(&lf->rf[i]);
277: }
278: for (i=0;i<lf->nc;i++) {
279: ISDestroy(&lf->cf[i]);
280: }
281: PetscFree2(lf->rf,lf->cf);
282: PetscFree(lf);
283: return(0);
284: }
286: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
287: {
288: Mat B,lB;
292: if (reuse != MAT_REUSE_MATRIX) {
293: ISLocalToGlobalMapping rl2g,cl2g;
294: PetscInt bs;
295: IS is;
297: MatGetBlockSize(A,&bs);
298: ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is);
299: if (bs > 1) {
300: IS is2;
301: PetscInt i,*aux;
303: ISGetLocalSize(is,&i);
304: ISGetIndices(is,(const PetscInt**)&aux);
305: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2);
306: ISRestoreIndices(is,(const PetscInt**)&aux);
307: ISDestroy(&is);
308: is = is2;
309: }
310: ISSetIdentity(is);
311: ISLocalToGlobalMappingCreateIS(is,&rl2g);
312: ISDestroy(&is);
313: ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->n/bs,0,1,&is);
314: if (bs > 1) {
315: IS is2;
316: PetscInt i,*aux;
318: ISGetLocalSize(is,&i);
319: ISGetIndices(is,(const PetscInt**)&aux);
320: ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2);
321: ISRestoreIndices(is,(const PetscInt**)&aux);
322: ISDestroy(&is);
323: is = is2;
324: }
325: ISSetIdentity(is);
326: ISLocalToGlobalMappingCreateIS(is,&cl2g);
327: ISDestroy(&is);
328: MatCreateIS(PetscObjectComm((PetscObject)A),bs,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,rl2g,cl2g,&B);
329: ISLocalToGlobalMappingDestroy(&rl2g);
330: ISLocalToGlobalMappingDestroy(&cl2g);
331: MatDuplicate(A,MAT_COPY_VALUES,&lB);
332: if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
333: } else {
334: B = *newmat;
335: PetscObjectReference((PetscObject)A);
336: lB = A;
337: }
338: MatISSetLocalMat(B,lB);
339: MatDestroy(&lB);
340: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
341: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
342: if (reuse == MAT_INPLACE_MATRIX) {
343: MatHeaderReplace(A,&B);
344: }
345: return(0);
346: }
348: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
349: {
350: Mat_IS *matis = (Mat_IS*)(A->data);
351: PetscScalar *aa;
352: const PetscInt *ii,*jj;
353: PetscInt i,n,m;
354: PetscInt *ecount,**eneighs;
355: PetscBool flg;
359: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
360: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
361: ISLocalToGlobalMappingGetNodeInfo(A->rmap->mapping,&n,&ecount,&eneighs);
362: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D != %D",m,n);
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(A->rmap->mapping,&n,&ecount,&eneighs);
382: MatSeqAIJRestoreArray(matis->A,&aa);
383: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
384: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
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_",0};
399: MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
400: MatPartitioning part;
401: PetscSF sf;
402: PetscErrorCode ierr;
405: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatIS l2g disassembling options","Mat");
406: PetscOptionsEnum("-mat_is_disassemble_l2g_type","Type of local-to-global mapping to be used for disassembling","MatISDisassemblel2gType",MatISDisassemblel2gTypes,(PetscEnum)mode,(PetscEnum*)&mode,NULL);
407: PetscOptionsEnd();
408: if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
409: MatGetLocalToGlobalMapping(A,l2g,NULL);
410: return(0);
411: }
412: PetscObjectGetComm((PetscObject)A,&comm);
413: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij);
414: PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij);
415: MatGetBlockSize(A,&bs);
416: switch (mode) {
417: case MAT_IS_DISASSEMBLE_L2G_ND:
418: MatPartitioningCreate(comm,&part);
419: MatPartitioningSetAdjacency(part,A);
420: PetscObjectSetOptionsPrefix((PetscObject)part,((PetscObject)A)->prefix);
421: MatPartitioningSetFromOptions(part);
422: MatPartitioningApplyND(part,&ndmap);
423: MatPartitioningDestroy(&part);
424: ISBuildTwoSided(ndmap,NULL,&ndsub);
425: MatMPIAIJSetUseScalableIncreaseOverlap(A,PETSC_TRUE);
426: MatIncreaseOverlap(A,1,&ndsub,1);
427: ISLocalToGlobalMappingCreateIS(ndsub,l2g);
429: /* it may happen that a separator node is not properly shared */
430: ISLocalToGlobalMappingGetNodeInfo(*l2g,&nl,&ncount,NULL);
431: PetscSFCreate(comm,&sf);
432: ISLocalToGlobalMappingGetIndices(*l2g,&garray);
433: PetscSFSetGraphLayout(sf,A->rmap,nl,NULL,PETSC_OWN_POINTER,garray);
434: ISLocalToGlobalMappingRestoreIndices(*l2g,&garray);
435: PetscCalloc1(A->rmap->n,&ndmapc);
436: PetscSFReduceBegin(sf,MPIU_INT,ncount,ndmapc,MPIU_REPLACE);
437: PetscSFReduceEnd(sf,MPIU_INT,ncount,ndmapc,MPIU_REPLACE);
438: ISLocalToGlobalMappingRestoreNodeInfo(*l2g,NULL,&ncount,NULL);
439: ISGetIndices(ndmap,&ndmapi);
440: for (i = 0, cnt = 0; i < A->rmap->n; i++)
441: if (ndmapi[i] < 0 && ndmapc[i] < 2)
442: cnt++;
444: MPIU_Allreduce(&cnt,&i,1,MPIU_INT,MPI_MAX,comm);
445: if (i) { /* we detected isolated separator nodes */
446: Mat A2,A3;
447: IS *workis,is2;
448: PetscScalar *vals;
449: PetscInt gcnt = i,*dnz,*onz,j,*lndmapi;
450: ISLocalToGlobalMapping ll2g;
451: PetscBool flg;
452: const PetscInt *ii,*jj;
454: /* communicate global id of separators */
455: MatPreallocateInitialize(comm,A->rmap->n,A->cmap->n,dnz,onz);
456: for (i = 0, cnt = 0; i < A->rmap->n; i++)
457: dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;
459: PetscMalloc1(nl,&lndmapi);
460: PetscSFBcastBegin(sf,MPIU_INT,dnz,lndmapi);
462: /* compute adjacency of isolated separators node */
463: PetscMalloc1(gcnt,&workis);
464: for (i = 0, cnt = 0; i < A->rmap->n; i++) {
465: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
466: ISCreateStride(comm,1,i+A->rmap->rstart,1,&workis[cnt++]);
467: }
468: }
469: for (i = cnt; i < gcnt; i++) {
470: ISCreateStride(comm,0,0,1,&workis[i]);
471: }
472: for (i = 0; i < gcnt; i++) {
473: PetscObjectSetName((PetscObject)workis[i],"ISOLATED");
474: ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
475: }
477: /* no communications since all the ISes correspond to locally owned rows */
478: MatIncreaseOverlap(A,gcnt,workis,1);
480: /* end communicate global id of separators */
481: PetscSFBcastEnd(sf,MPIU_INT,dnz,lndmapi);
483: /* communicate new layers : create a matrix and transpose it */
484: PetscMemzero(dnz,A->rmap->n*sizeof(*dnz));
485: PetscMemzero(onz,A->rmap->n*sizeof(*onz));
486: for (i = 0, j = 0; i < A->rmap->n; i++) {
487: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
488: const PetscInt* idxs;
489: PetscInt s;
491: ISGetLocalSize(workis[j],&s);
492: ISGetIndices(workis[j],&idxs);
493: MatPreallocateSet(i+A->rmap->rstart,s,idxs,dnz,onz);
494: j++;
495: }
496: }
497: if (j != cnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %D != %D",j,cnt);
499: for (i = 0; i < gcnt; i++) {
500: PetscObjectSetName((PetscObject)workis[i],"EXTENDED");
501: ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
502: }
504: for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j,dnz[i]+onz[i]);
505: PetscMalloc1(j,&vals);
506: for (i = 0; i < j; i++) vals[i] = 1.0;
508: MatCreate(comm,&A2);
509: MatSetType(A2,MATMPIAIJ);
510: MatSetSizes(A2,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
511: MatMPIAIJSetPreallocation(A2,0,dnz,0,onz);
512: MatSetOption(A2,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
513: for (i = 0, j = 0; i < A2->rmap->n; i++) {
514: PetscInt row = i+A2->rmap->rstart,s = dnz[i] + onz[i];
515: const PetscInt* idxs;
517: if (s) {
518: ISGetIndices(workis[j],&idxs);
519: MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES);
520: ISRestoreIndices(workis[j],&idxs);
521: j++;
522: }
523: }
524: if (j != cnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %D != %D",j,cnt);
525: PetscFree(vals);
526: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
527: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
528: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
530: /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
531: for (i = 0, j = 0; i < nl; i++)
532: if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
533: ISCreateGeneral(comm,j,lndmapi,PETSC_USE_POINTER,&is);
534: MatMPIAIJGetLocalMatCondensed(A2,MAT_INITIAL_MATRIX,&is,NULL,&A3);
535: ISDestroy(&is);
536: MatDestroy(&A2);
538: /* extend local to global map to include connected isolated separators */
539: PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is);
540: if (!is) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing column map");
541: ISLocalToGlobalMappingCreateIS(is,&ll2g);
542: MatGetRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
543: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
544: ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is);
545: MatRestoreRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
546: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
547: ISLocalToGlobalMappingApplyIS(ll2g,is,&is2);
548: ISDestroy(&is);
549: ISLocalToGlobalMappingDestroy(&ll2g);
551: /* add new nodes to the local-to-global map */
552: ISLocalToGlobalMappingDestroy(l2g);
553: ISExpand(ndsub,is2,&is);
554: ISDestroy(&is2);
555: ISLocalToGlobalMappingCreateIS(is,l2g);
556: ISDestroy(&is);
558: MatDestroy(&A3);
559: PetscFree(lndmapi);
560: MatPreallocateFinalize(dnz,onz);
561: for (i = 0; i < gcnt; i++) {
562: ISDestroy(&workis[i]);
563: }
564: PetscFree(workis);
565: }
566: ISRestoreIndices(ndmap,&ndmapi);
567: PetscSFDestroy(&sf);
568: PetscFree(ndmapc);
569: ISDestroy(&ndmap);
570: ISDestroy(&ndsub);
571: ISLocalToGlobalMappingSetBlockSize(*l2g,bs);
572: ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view");
573: break;
574: case MAT_IS_DISASSEMBLE_L2G_NATURAL:
575: if (ismpiaij) {
576: MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
577: } else if (ismpibaij) {
578: MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
579: } else SETERRQ1(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
580: if (!garray) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present");
581: if (A->rmap->n) {
582: PetscInt dc,oc,stc,*aux;
584: MatGetLocalSize(A,NULL,&dc);
585: MatGetLocalSize(Ao,NULL,&oc);
586: MatGetOwnershipRangeColumn(A,&stc,NULL);
587: PetscMalloc1((dc+oc)/bs,&aux);
588: for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs;
589: for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
590: ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
591: } else {
592: ISCreateBlock(comm,1,0,NULL,PETSC_OWN_POINTER,&is);
593: }
594: ISLocalToGlobalMappingCreateIS(is,l2g);
595: ISDestroy(&is);
596: break;
597: default:
598: SETERRQ1(comm,PETSC_ERR_ARG_WRONG,"Unsupported l2g disassembling type %D",mode);
599: }
600: return(0);
601: }
603: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
604: {
605: Mat lA,Ad,Ao,B = NULL;
606: ISLocalToGlobalMapping rl2g,cl2g;
607: IS is;
608: MPI_Comm comm;
609: void *ptrs[2];
610: const char *names[2] = {"_convert_csr_aux","_convert_csr_data"};
611: const PetscInt *garray;
612: PetscScalar *dd,*od,*aa,*data;
613: const PetscInt *di,*dj,*oi,*oj;
614: const PetscInt *odi,*odj,*ooi,*ooj;
615: PetscInt *aux,*ii,*jj;
616: PetscInt bs,lc,dr,dc,oc,str,stc,nnz,i,jd,jo,cum;
617: PetscBool flg,ismpiaij,ismpibaij,was_inplace = PETSC_FALSE;
618: PetscMPIInt size;
619: PetscErrorCode ierr;
622: PetscObjectGetComm((PetscObject)A,&comm);
623: MPI_Comm_size(comm,&size);
624: if (size == 1) {
625: MatConvert_SeqXAIJ_IS(A,type,reuse,newmat);
626: return(0);
627: }
628: if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
629: MatMPIXAIJComputeLocalToGlobalMapping_Private(A,&rl2g);
630: MatCreate(comm,&B);
631: MatSetType(B,MATIS);
632: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
633: MatSetLocalToGlobalMapping(B,rl2g,rl2g);
634: MatGetBlockSize(A,&bs);
635: MatSetBlockSize(B,bs);
636: ISLocalToGlobalMappingDestroy(&rl2g);
637: if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
638: reuse = MAT_REUSE_MATRIX;
639: }
640: if (reuse == MAT_REUSE_MATRIX) {
641: Mat *newlA, lA;
642: IS rows, cols;
643: const PetscInt *ridx, *cidx;
644: PetscInt rbs, cbs, nr, nc;
646: if (!B) B = *newmat;
647: MatGetLocalToGlobalMapping(B,&rl2g,&cl2g);
648: ISLocalToGlobalMappingGetBlockIndices(rl2g,&ridx);
649: ISLocalToGlobalMappingGetBlockIndices(cl2g,&cidx);
650: ISLocalToGlobalMappingGetSize(rl2g,&nr);
651: ISLocalToGlobalMappingGetSize(cl2g,&nc);
652: ISLocalToGlobalMappingGetBlockSize(rl2g,&rbs);
653: ISLocalToGlobalMappingGetBlockSize(cl2g,&cbs);
654: ISCreateBlock(comm,rbs,nr/rbs,ridx,PETSC_USE_POINTER,&rows);
655: if (rl2g != cl2g) {
656: ISCreateBlock(comm,cbs,nc/cbs,cidx,PETSC_USE_POINTER,&cols);
657: } else {
658: PetscObjectReference((PetscObject)rows);
659: cols = rows;
660: }
661: MatISGetLocalMat(B,&lA);
662: MatCreateSubMatrices(A,1,&rows,&cols,MAT_INITIAL_MATRIX,&newlA);
663: MatConvert(newlA[0],MATSEQAIJ,MAT_INPLACE_MATRIX,&newlA[0]);
664: ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&ridx);
665: ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&cidx);
666: ISDestroy(&rows);
667: ISDestroy(&cols);
668: if (!lA->preallocated) { /* first time */
669: MatDuplicate(newlA[0],MAT_COPY_VALUES,&lA);
670: MatISSetLocalMat(B,lA);
671: PetscObjectDereference((PetscObject)lA);
672: }
673: MatCopy(newlA[0],lA,SAME_NONZERO_PATTERN);
674: MatDestroySubMatrices(1,&newlA);
675: MatISScaleDisassembling_Private(B);
676: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
677: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
678: if (was_inplace) { MatHeaderReplace(A,&B); }
679: else *newmat = B;
680: return(0);
681: }
682: /* rectangular case, just compress out the column space */
683: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij);
684: PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij);
685: if (ismpiaij) {
686: bs = 1;
687: MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
688: } else if (ismpibaij) {
689: MatGetBlockSize(A,&bs);
690: MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
691: MatConvert(Ad,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ad);
692: MatConvert(Ao,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ao);
693: } else SETERRQ1(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
694: MatSeqAIJGetArray(Ad,&dd);
695: MatSeqAIJGetArray(Ao,&od);
696: if (!garray) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present");
698: /* access relevant information from MPIAIJ */
699: MatGetOwnershipRange(A,&str,NULL);
700: MatGetOwnershipRangeColumn(A,&stc,NULL);
701: MatGetLocalSize(A,&dr,&dc);
702: MatGetLocalSize(Ao,NULL,&oc);
703: MatGetRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&di,&dj,&flg);
704: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
705: MatGetRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&oi,&oj,&flg);
706: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
707: nnz = di[dr] + oi[dr];
708: /* store original pointers to be restored later */
709: odi = di; odj = dj; ooi = oi; ooj = oj;
711: /* generate l2g maps for rows and cols */
712: ISCreateStride(comm,dr/bs,str/bs,1,&is);
713: if (bs > 1) {
714: IS is2;
716: ISGetLocalSize(is,&i);
717: ISGetIndices(is,(const PetscInt**)&aux);
718: ISCreateBlock(comm,bs,i,aux,PETSC_COPY_VALUES,&is2);
719: ISRestoreIndices(is,(const PetscInt**)&aux);
720: ISDestroy(&is);
721: is = is2;
722: }
723: ISLocalToGlobalMappingCreateIS(is,&rl2g);
724: ISDestroy(&is);
725: if (dr) {
726: PetscMalloc1((dc+oc)/bs,&aux);
727: for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs;
728: for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
729: ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
730: lc = dc+oc;
731: } else {
732: ISCreateBlock(comm,bs,0,NULL,PETSC_OWN_POINTER,&is);
733: lc = 0;
734: }
735: ISLocalToGlobalMappingCreateIS(is,&cl2g);
736: ISDestroy(&is);
738: /* create MATIS object */
739: MatCreate(comm,&B);
740: MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
741: MatSetType(B,MATIS);
742: MatSetBlockSize(B,bs);
743: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
744: ISLocalToGlobalMappingDestroy(&rl2g);
745: ISLocalToGlobalMappingDestroy(&cl2g);
747: /* merge local matrices */
748: PetscMalloc1(nnz+dr+1,&aux);
749: PetscMalloc1(nnz,&data);
750: ii = aux;
751: jj = aux+dr+1;
752: aa = data;
753: *ii = *(di++) + *(oi++);
754: for (jd=0,jo=0,cum=0;*ii<nnz;cum++)
755: {
756: for (;jd<*di;jd++) { *jj++ = *dj++; *aa++ = *dd++; }
757: for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; }
758: *(++ii) = *(di++) + *(oi++);
759: }
760: for (;cum<dr;cum++) *(++ii) = nnz;
762: MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg);
763: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
764: MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg);
765: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
766: MatSeqAIJRestoreArray(Ad,&dd);
767: MatSeqAIJRestoreArray(Ao,&od);
769: ii = aux;
770: jj = aux+dr+1;
771: aa = data;
772: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);
774: /* create containers to destroy the data */
775: ptrs[0] = aux;
776: ptrs[1] = data;
777: for (i=0; i<2; i++) {
778: PetscContainer c;
780: PetscContainerCreate(PETSC_COMM_SELF,&c);
781: PetscContainerSetPointer(c,ptrs[i]);
782: PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);
783: PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
784: PetscContainerDestroy(&c);
785: }
786: if (ismpibaij) { /* destroy converted local matrices */
787: MatDestroy(&Ad);
788: MatDestroy(&Ao);
789: }
791: /* finalize matrix */
792: MatISSetLocalMat(B,lA);
793: MatDestroy(&lA);
794: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
795: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
796: if (reuse == MAT_INPLACE_MATRIX) {
797: MatHeaderReplace(A,&B);
798: } else *newmat = B;
799: return(0);
800: }
802: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
803: {
804: Mat **nest,*snest,**rnest,lA,B;
805: IS *iscol,*isrow,*islrow,*islcol;
806: ISLocalToGlobalMapping rl2g,cl2g;
807: MPI_Comm comm;
808: PetscInt *lr,*lc,*l2gidxs;
809: PetscInt i,j,nr,nc,rbs,cbs;
810: PetscBool convert,lreuse,*istrans;
811: PetscErrorCode ierr;
814: MatNestGetSubMats(A,&nr,&nc,&nest);
815: lreuse = PETSC_FALSE;
816: rnest = NULL;
817: if (reuse == MAT_REUSE_MATRIX) {
818: PetscBool ismatis,isnest;
820: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
821: if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
822: MatISGetLocalMat(*newmat,&lA);
823: PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
824: if (isnest) {
825: MatNestGetSubMats(lA,&i,&j,&rnest);
826: lreuse = (PetscBool)(i == nr && j == nc);
827: if (!lreuse) rnest = NULL;
828: }
829: }
830: PetscObjectGetComm((PetscObject)A,&comm);
831: PetscCalloc2(nr,&lr,nc,&lc);
832: PetscCalloc6(nr,&isrow,nc,&iscol,
833: nr,&islrow,nc,&islcol,
834: nr*nc,&snest,nr*nc,&istrans);
835: MatNestGetISs(A,isrow,iscol);
836: for (i=0;i<nr;i++) {
837: for (j=0;j<nc;j++) {
838: PetscBool ismatis;
839: PetscInt l1,l2,lb1,lb2,ij=i*nc+j;
841: /* Null matrix pointers are allowed in MATNEST */
842: if (!nest[i][j]) continue;
844: /* Nested matrices should be of type MATIS */
845: PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
846: if (istrans[ij]) {
847: Mat T,lT;
848: MatTransposeGetMat(nest[i][j],&T);
849: PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
850: if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) (transposed) is not of type MATIS",i,j);
851: MatISGetLocalMat(T,&lT);
852: MatCreateTranspose(lT,&snest[ij]);
853: } else {
854: PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
855: if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
856: MatISGetLocalMat(nest[i][j],&snest[ij]);
857: }
859: /* Check compatibility of local sizes */
860: MatGetSize(snest[ij],&l1,&l2);
861: MatGetBlockSizes(snest[ij],&lb1,&lb2);
862: if (!l1 || !l2) continue;
863: if (lr[i] && l1 != lr[i]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lr[i],l1);
864: if (lc[j] && l2 != lc[j]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lc[j],l2);
865: lr[i] = l1;
866: lc[j] = l2;
868: /* check compatibilty for local matrix reusage */
869: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
870: }
871: }
873: #if defined (PETSC_USE_DEBUG)
874: /* Check compatibility of l2g maps for rows */
875: for (i=0;i<nr;i++) {
876: rl2g = NULL;
877: for (j=0;j<nc;j++) {
878: PetscInt n1,n2;
880: if (!nest[i][j]) continue;
881: if (istrans[i*nc+j]) {
882: Mat T;
884: MatTransposeGetMat(nest[i][j],&T);
885: MatGetLocalToGlobalMapping(T,NULL,&cl2g);
886: } else {
887: MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
888: }
889: ISLocalToGlobalMappingGetSize(cl2g,&n1);
890: if (!n1) continue;
891: if (!rl2g) {
892: rl2g = cl2g;
893: } else {
894: const PetscInt *idxs1,*idxs2;
895: PetscBool same;
897: ISLocalToGlobalMappingGetSize(rl2g,&n2);
898: if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap size %D != %D",i,j,n1,n2);
899: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
900: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
901: PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
902: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
903: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
904: if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap",i,j);
905: }
906: }
907: }
908: /* Check compatibility of l2g maps for columns */
909: for (i=0;i<nc;i++) {
910: rl2g = NULL;
911: for (j=0;j<nr;j++) {
912: PetscInt n1,n2;
914: if (!nest[j][i]) continue;
915: if (istrans[j*nc+i]) {
916: Mat T;
918: MatTransposeGetMat(nest[j][i],&T);
919: MatGetLocalToGlobalMapping(T,&cl2g,NULL);
920: } else {
921: MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
922: }
923: ISLocalToGlobalMappingGetSize(cl2g,&n1);
924: if (!n1) continue;
925: if (!rl2g) {
926: rl2g = cl2g;
927: } else {
928: const PetscInt *idxs1,*idxs2;
929: PetscBool same;
931: ISLocalToGlobalMappingGetSize(rl2g,&n2);
932: if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap size %D != %D",j,i,n1,n2);
933: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
934: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
935: PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
936: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
937: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
938: if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap",j,i);
939: }
940: }
941: }
942: #endif
944: B = NULL;
945: if (reuse != MAT_REUSE_MATRIX) {
946: PetscInt stl;
948: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
949: for (i=0,stl=0;i<nr;i++) stl += lr[i];
950: PetscMalloc1(stl,&l2gidxs);
951: for (i=0,stl=0;i<nr;i++) {
952: Mat usedmat;
953: Mat_IS *matis;
954: const PetscInt *idxs;
956: /* local IS for local NEST */
957: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
959: /* l2gmap */
960: j = 0;
961: usedmat = nest[i][j];
962: while (!usedmat && j < nc-1) usedmat = nest[i][++j];
963: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");
965: if (istrans[i*nc+j]) {
966: Mat T;
967: MatTransposeGetMat(usedmat,&T);
968: usedmat = T;
969: }
970: matis = (Mat_IS*)(usedmat->data);
971: ISGetIndices(isrow[i],&idxs);
972: if (istrans[i*nc+j]) {
973: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
974: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
975: } else {
976: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
977: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
978: }
979: ISRestoreIndices(isrow[i],&idxs);
980: stl += lr[i];
981: }
982: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);
984: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
985: for (i=0,stl=0;i<nc;i++) stl += lc[i];
986: PetscMalloc1(stl,&l2gidxs);
987: for (i=0,stl=0;i<nc;i++) {
988: Mat usedmat;
989: Mat_IS *matis;
990: const PetscInt *idxs;
992: /* local IS for local NEST */
993: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
995: /* l2gmap */
996: j = 0;
997: usedmat = nest[j][i];
998: while (!usedmat && j < nr-1) usedmat = nest[++j][i];
999: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat");
1000: if (istrans[j*nc+i]) {
1001: Mat T;
1002: MatTransposeGetMat(usedmat,&T);
1003: usedmat = T;
1004: }
1005: matis = (Mat_IS*)(usedmat->data);
1006: ISGetIndices(iscol[i],&idxs);
1007: if (istrans[j*nc+i]) {
1008: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
1009: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
1010: } else {
1011: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
1012: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
1013: }
1014: ISRestoreIndices(iscol[i],&idxs);
1015: stl += lc[i];
1016: }
1017: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);
1019: /* Create MATIS */
1020: MatCreate(comm,&B);
1021: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1022: MatGetBlockSizes(A,&rbs,&cbs);
1023: MatSetBlockSizes(B,rbs,cbs);
1024: MatSetType(B,MATIS);
1025: MatISSetLocalMatType(B,MATNEST);
1026: { /* hack : avoid setup of scatters */
1027: Mat_IS *matis = (Mat_IS*)(B->data);
1028: matis->islocalref = PETSC_TRUE;
1029: }
1030: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
1031: ISLocalToGlobalMappingDestroy(&rl2g);
1032: ISLocalToGlobalMappingDestroy(&cl2g);
1033: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1034: MatNestSetVecType(lA,VECNEST);
1035: for (i=0;i<nr*nc;i++) {
1036: if (istrans[i]) {
1037: MatDestroy(&snest[i]);
1038: }
1039: }
1040: MatISSetLocalMat(B,lA);
1041: MatDestroy(&lA);
1042: { /* hack : setup of scatters done here */
1043: Mat_IS *matis = (Mat_IS*)(B->data);
1045: matis->islocalref = PETSC_FALSE;
1046: MatISSetUpScatters_Private(B);
1047: }
1048: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1049: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1050: if (reuse == MAT_INPLACE_MATRIX) {
1051: MatHeaderReplace(A,&B);
1052: } else {
1053: *newmat = B;
1054: }
1055: } else {
1056: if (lreuse) {
1057: MatISGetLocalMat(*newmat,&lA);
1058: for (i=0;i<nr;i++) {
1059: for (j=0;j<nc;j++) {
1060: if (snest[i*nc+j]) {
1061: MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
1062: if (istrans[i*nc+j]) {
1063: MatDestroy(&snest[i*nc+j]);
1064: }
1065: }
1066: }
1067: }
1068: } else {
1069: PetscInt stl;
1070: for (i=0,stl=0;i<nr;i++) {
1071: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
1072: stl += lr[i];
1073: }
1074: for (i=0,stl=0;i<nc;i++) {
1075: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
1076: stl += lc[i];
1077: }
1078: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1079: for (i=0;i<nr*nc;i++) {
1080: if (istrans[i]) {
1081: MatDestroy(&snest[i]);
1082: }
1083: }
1084: MatISSetLocalMat(*newmat,lA);
1085: MatDestroy(&lA);
1086: }
1087: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1088: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1089: }
1091: /* Create local matrix in MATNEST format */
1092: convert = PETSC_FALSE;
1093: PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
1094: if (convert) {
1095: Mat M;
1096: MatISLocalFields lf;
1097: PetscContainer c;
1099: MatISGetLocalMat(*newmat,&lA);
1100: MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
1101: MatISSetLocalMat(*newmat,M);
1102: MatDestroy(&M);
1104: /* attach local fields to the matrix */
1105: PetscNew(&lf);
1106: PetscCalloc2(nr,&lf->rf,nc,&lf->cf);
1107: for (i=0;i<nr;i++) {
1108: PetscInt n,st;
1110: ISGetLocalSize(islrow[i],&n);
1111: ISStrideGetInfo(islrow[i],&st,NULL);
1112: ISCreateStride(comm,n,st,1,&lf->rf[i]);
1113: }
1114: for (i=0;i<nc;i++) {
1115: PetscInt n,st;
1117: ISGetLocalSize(islcol[i],&n);
1118: ISStrideGetInfo(islcol[i],&st,NULL);
1119: ISCreateStride(comm,n,st,1,&lf->cf[i]);
1120: }
1121: lf->nr = nr;
1122: lf->nc = nc;
1123: PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
1124: PetscContainerSetPointer(c,lf);
1125: PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
1126: PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
1127: PetscContainerDestroy(&c);
1128: }
1130: /* Free workspace */
1131: for (i=0;i<nr;i++) {
1132: ISDestroy(&islrow[i]);
1133: }
1134: for (i=0;i<nc;i++) {
1135: ISDestroy(&islcol[i]);
1136: }
1137: PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
1138: PetscFree2(lr,lc);
1139: return(0);
1140: }
1142: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1143: {
1144: Mat_IS *matis = (Mat_IS*)A->data;
1145: Vec ll,rr;
1146: const PetscScalar *Y,*X;
1147: PetscScalar *x,*y;
1148: PetscErrorCode ierr;
1151: if (l) {
1152: ll = matis->y;
1153: VecGetArrayRead(l,&Y);
1154: VecGetArray(ll,&y);
1155: PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y);
1156: } else {
1157: ll = NULL;
1158: }
1159: if (r) {
1160: rr = matis->x;
1161: VecGetArrayRead(r,&X);
1162: VecGetArray(rr,&x);
1163: PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x);
1164: } else {
1165: rr = NULL;
1166: }
1167: if (ll) {
1168: PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y);
1169: VecRestoreArrayRead(l,&Y);
1170: VecRestoreArray(ll,&y);
1171: }
1172: if (rr) {
1173: PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x);
1174: VecRestoreArrayRead(r,&X);
1175: VecRestoreArray(rr,&x);
1176: }
1177: MatDiagonalScale(matis->A,ll,rr);
1178: return(0);
1179: }
1181: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
1182: {
1183: Mat_IS *matis = (Mat_IS*)A->data;
1184: MatInfo info;
1185: PetscReal isend[6],irecv[6];
1186: PetscInt bs;
1190: MatGetBlockSize(A,&bs);
1191: if (matis->A->ops->getinfo) {
1192: MatGetInfo(matis->A,MAT_LOCAL,&info);
1193: isend[0] = info.nz_used;
1194: isend[1] = info.nz_allocated;
1195: isend[2] = info.nz_unneeded;
1196: isend[3] = info.memory;
1197: isend[4] = info.mallocs;
1198: } else {
1199: isend[0] = 0.;
1200: isend[1] = 0.;
1201: isend[2] = 0.;
1202: isend[3] = 0.;
1203: isend[4] = 0.;
1204: }
1205: isend[5] = matis->A->num_ass;
1206: if (flag == MAT_LOCAL) {
1207: ginfo->nz_used = isend[0];
1208: ginfo->nz_allocated = isend[1];
1209: ginfo->nz_unneeded = isend[2];
1210: ginfo->memory = isend[3];
1211: ginfo->mallocs = isend[4];
1212: ginfo->assemblies = isend[5];
1213: } else if (flag == MAT_GLOBAL_MAX) {
1214: MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));
1216: ginfo->nz_used = irecv[0];
1217: ginfo->nz_allocated = irecv[1];
1218: ginfo->nz_unneeded = irecv[2];
1219: ginfo->memory = irecv[3];
1220: ginfo->mallocs = irecv[4];
1221: ginfo->assemblies = irecv[5];
1222: } else if (flag == MAT_GLOBAL_SUM) {
1223: MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));
1225: ginfo->nz_used = irecv[0];
1226: ginfo->nz_allocated = irecv[1];
1227: ginfo->nz_unneeded = irecv[2];
1228: ginfo->memory = irecv[3];
1229: ginfo->mallocs = irecv[4];
1230: ginfo->assemblies = A->num_ass;
1231: }
1232: ginfo->block_size = bs;
1233: ginfo->fill_ratio_given = 0;
1234: ginfo->fill_ratio_needed = 0;
1235: ginfo->factor_mallocs = 0;
1236: return(0);
1237: }
1239: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
1240: {
1241: Mat C,lC,lA;
1242: PetscErrorCode ierr;
1245: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1246: ISLocalToGlobalMapping rl2g,cl2g;
1247: MatCreate(PetscObjectComm((PetscObject)A),&C);
1248: MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
1249: MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
1250: MatSetType(C,MATIS);
1251: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
1252: MatSetLocalToGlobalMapping(C,cl2g,rl2g);
1253: } else {
1254: C = *B;
1255: }
1257: /* perform local transposition */
1258: MatISGetLocalMat(A,&lA);
1259: MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
1260: MatISSetLocalMat(C,lC);
1261: MatDestroy(&lC);
1263: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1264: *B = C;
1265: } else {
1266: MatHeaderMerge(A,&C);
1267: }
1268: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1269: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1270: return(0);
1271: }
1273: PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
1274: {
1275: Mat_IS *is = (Mat_IS*)A->data;
1279: if (D) { /* MatShift_IS pass D = NULL */
1280: VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1281: VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1282: }
1283: VecPointwiseDivide(is->y,is->y,is->counter);
1284: MatDiagonalSet(is->A,is->y,insmode);
1285: return(0);
1286: }
1288: PetscErrorCode MatShift_IS(Mat A,PetscScalar a)
1289: {
1290: Mat_IS *is = (Mat_IS*)A->data;
1294: VecSet(is->y,a);
1295: MatDiagonalSet_IS(A,NULL,ADD_VALUES);
1296: return(0);
1297: }
1299: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1300: {
1302: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1305: #if defined(PETSC_USE_DEBUG)
1306: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
1307: #endif
1308: ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
1309: ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
1310: MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1311: return(0);
1312: }
1314: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1315: {
1317: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1320: #if defined(PETSC_USE_DEBUG)
1321: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
1322: #endif
1323: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
1324: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
1325: MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1326: return(0);
1327: }
1329: static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
1330: {
1331: PetscInt *owners = map->range;
1332: PetscInt n = map->n;
1333: PetscSF sf;
1334: PetscInt *lidxs,*work = NULL;
1335: PetscSFNode *ridxs;
1336: PetscMPIInt rank;
1337: PetscInt r, p = 0, len = 0;
1341: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
1342: /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */
1343: MPI_Comm_rank(map->comm,&rank);
1344: PetscMalloc1(n,&lidxs);
1345: for (r = 0; r < n; ++r) lidxs[r] = -1;
1346: PetscMalloc1(N,&ridxs);
1347: for (r = 0; r < N; ++r) {
1348: const PetscInt idx = idxs[r];
1349: if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
1350: if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
1351: PetscLayoutFindOwner(map,idx,&p);
1352: }
1353: ridxs[r].rank = p;
1354: ridxs[r].index = idxs[r] - owners[p];
1355: }
1356: PetscSFCreate(map->comm,&sf);
1357: PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
1358: PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
1359: PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
1360: if (ogidxs) { /* communicate global idxs */
1361: PetscInt cum = 0,start,*work2;
1363: PetscMalloc1(n,&work);
1364: PetscCalloc1(N,&work2);
1365: for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
1366: MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
1367: start -= cum;
1368: cum = 0;
1369: for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
1370: PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);
1371: PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);
1372: PetscFree(work2);
1373: }
1374: PetscSFDestroy(&sf);
1375: /* Compress and put in indices */
1376: for (r = 0; r < n; ++r)
1377: if (lidxs[r] >= 0) {
1378: if (work) work[len] = work[r];
1379: lidxs[len++] = r;
1380: }
1381: if (on) *on = len;
1382: if (oidxs) *oidxs = lidxs;
1383: if (ogidxs) *ogidxs = work;
1384: return(0);
1385: }
1387: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
1388: {
1389: Mat locmat,newlocmat;
1390: Mat_IS *newmatis;
1391: #if defined(PETSC_USE_DEBUG)
1392: Vec rtest,ltest;
1393: const PetscScalar *array;
1394: #endif
1395: const PetscInt *idxs;
1396: PetscInt i,m,n;
1397: PetscErrorCode ierr;
1400: if (scall == MAT_REUSE_MATRIX) {
1401: PetscBool ismatis;
1403: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
1404: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
1405: newmatis = (Mat_IS*)(*newmat)->data;
1406: if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
1407: if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
1408: }
1409: /* irow and icol may not have duplicate entries */
1410: #if defined(PETSC_USE_DEBUG)
1411: MatCreateVecs(mat,<est,&rtest);
1412: ISGetLocalSize(irow,&n);
1413: ISGetIndices(irow,&idxs);
1414: for (i=0;i<n;i++) {
1415: VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
1416: }
1417: VecAssemblyBegin(rtest);
1418: VecAssemblyEnd(rtest);
1419: VecGetLocalSize(rtest,&n);
1420: VecGetOwnershipRange(rtest,&m,NULL);
1421: VecGetArrayRead(rtest,&array);
1422: for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Irow may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
1423: VecRestoreArrayRead(rtest,&array);
1424: ISRestoreIndices(irow,&idxs);
1425: ISGetLocalSize(icol,&n);
1426: ISGetIndices(icol,&idxs);
1427: for (i=0;i<n;i++) {
1428: VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
1429: }
1430: VecAssemblyBegin(ltest);
1431: VecAssemblyEnd(ltest);
1432: VecGetLocalSize(ltest,&n);
1433: VecGetOwnershipRange(ltest,&m,NULL);
1434: VecGetArrayRead(ltest,&array);
1435: for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Icol may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
1436: VecRestoreArrayRead(ltest,&array);
1437: ISRestoreIndices(icol,&idxs);
1438: VecDestroy(&rtest);
1439: VecDestroy(<est);
1440: #endif
1441: if (scall == MAT_INITIAL_MATRIX) {
1442: Mat_IS *matis = (Mat_IS*)mat->data;
1443: ISLocalToGlobalMapping rl2g;
1444: IS is;
1445: PetscInt *lidxs,*lgidxs,*newgidxs;
1446: PetscInt ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
1447: PetscBool cong;
1448: MPI_Comm comm;
1450: PetscObjectGetComm((PetscObject)mat,&comm);
1451: MatGetBlockSizes(mat,&arbs,&acbs);
1452: ISGetBlockSize(irow,&irbs);
1453: ISGetBlockSize(icol,&icbs);
1454: rbs = arbs == irbs ? irbs : 1;
1455: cbs = acbs == icbs ? icbs : 1;
1456: ISGetLocalSize(irow,&m);
1457: ISGetLocalSize(icol,&n);
1458: MatCreate(comm,newmat);
1459: MatSetType(*newmat,MATIS);
1460: MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
1461: MatSetBlockSizes(*newmat,rbs,cbs);
1462: /* communicate irow to their owners in the layout */
1463: ISGetIndices(irow,&idxs);
1464: PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
1465: ISRestoreIndices(irow,&idxs);
1466: PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));
1467: for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
1468: PetscFree(lidxs);
1469: PetscFree(lgidxs);
1470: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1471: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1472: for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
1473: PetscMalloc1(newloc,&newgidxs);
1474: PetscMalloc1(newloc,&lidxs);
1475: for (i=0,newloc=0;i<matis->sf->nleaves;i++)
1476: if (matis->sf_leafdata[i]) {
1477: lidxs[newloc] = i;
1478: newgidxs[newloc++] = matis->sf_leafdata[i]-1;
1479: }
1480: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1481: ISLocalToGlobalMappingCreateIS(is,&rl2g);
1482: ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);
1483: ISDestroy(&is);
1484: /* local is to extract local submatrix */
1485: newmatis = (Mat_IS*)(*newmat)->data;
1486: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
1487: MatHasCongruentLayouts(mat,&cong);
1488: if (cong && irow == icol && matis->csf == matis->sf) {
1489: MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
1490: PetscObjectReference((PetscObject)newmatis->getsub_ris);
1491: newmatis->getsub_cis = newmatis->getsub_ris;
1492: } else {
1493: ISLocalToGlobalMapping cl2g;
1495: /* communicate icol to their owners in the layout */
1496: ISGetIndices(icol,&idxs);
1497: PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
1498: ISRestoreIndices(icol,&idxs);
1499: PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));
1500: for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
1501: PetscFree(lidxs);
1502: PetscFree(lgidxs);
1503: PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
1504: PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
1505: for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
1506: PetscMalloc1(newloc,&newgidxs);
1507: PetscMalloc1(newloc,&lidxs);
1508: for (i=0,newloc=0;i<matis->csf->nleaves;i++)
1509: if (matis->csf_leafdata[i]) {
1510: lidxs[newloc] = i;
1511: newgidxs[newloc++] = matis->csf_leafdata[i]-1;
1512: }
1513: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1514: ISLocalToGlobalMappingCreateIS(is,&cl2g);
1515: ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);
1516: ISDestroy(&is);
1517: /* local is to extract local submatrix */
1518: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
1519: MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
1520: ISLocalToGlobalMappingDestroy(&cl2g);
1521: }
1522: ISLocalToGlobalMappingDestroy(&rl2g);
1523: } else {
1524: MatISGetLocalMat(*newmat,&newlocmat);
1525: }
1526: MatISGetLocalMat(mat,&locmat);
1527: newmatis = (Mat_IS*)(*newmat)->data;
1528: MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
1529: if (scall == MAT_INITIAL_MATRIX) {
1530: MatISSetLocalMat(*newmat,newlocmat);
1531: MatDestroy(&newlocmat);
1532: }
1533: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1534: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1535: return(0);
1536: }
1538: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
1539: {
1540: Mat_IS *a = (Mat_IS*)A->data,*b;
1541: PetscBool ismatis;
1545: PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
1546: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
1547: b = (Mat_IS*)B->data;
1548: MatCopy(a->A,b->A,str);
1549: PetscObjectStateIncrease((PetscObject)B);
1550: return(0);
1551: }
1553: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d)
1554: {
1555: Vec v;
1556: const PetscScalar *array;
1557: PetscInt i,n;
1558: PetscErrorCode ierr;
1561: *missing = PETSC_FALSE;
1562: MatCreateVecs(A,NULL,&v);
1563: MatGetDiagonal(A,v);
1564: VecGetLocalSize(v,&n);
1565: VecGetArrayRead(v,&array);
1566: for (i=0;i<n;i++) if (array[i] == 0.) break;
1567: VecRestoreArrayRead(v,&array);
1568: VecDestroy(&v);
1569: if (i != n) *missing = PETSC_TRUE;
1570: if (d) {
1571: *d = -1;
1572: if (*missing) {
1573: PetscInt rstart;
1574: MatGetOwnershipRange(A,&rstart,NULL);
1575: *d = i+rstart;
1576: }
1577: }
1578: return(0);
1579: }
1581: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1582: {
1583: Mat_IS *matis = (Mat_IS*)(B->data);
1584: const PetscInt *gidxs;
1585: PetscInt nleaves;
1589: if (matis->sf) return(0);
1590: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
1591: ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
1592: ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
1593: PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1594: ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
1595: PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
1596: if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
1597: ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
1598: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
1599: ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
1600: PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1601: ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
1602: PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
1603: } else {
1604: matis->csf = matis->sf;
1605: matis->csf_leafdata = matis->sf_leafdata;
1606: matis->csf_rootdata = matis->sf_rootdata;
1607: }
1608: return(0);
1609: }
1611: /*@
1612: MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP.
1614: Collective on MPI_Comm
1616: Input Parameters:
1617: + A - the matrix
1618: - store - the boolean flag
1620: Level: advanced
1622: Notes:
1624: .keywords: matrix
1626: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
1627: @*/
1628: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1629: {
1636: PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
1637: return(0);
1638: }
1640: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1641: {
1642: Mat_IS *matis = (Mat_IS*)(A->data);
1646: matis->storel2l = store;
1647: if (!store) {
1648: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
1649: }
1650: return(0);
1651: }
1653: /*@
1654: MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1656: Collective on MPI_Comm
1658: Input Parameters:
1659: + A - the matrix
1660: - fix - the boolean flag
1662: Level: advanced
1664: Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process.
1666: .keywords: matrix
1668: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
1669: @*/
1670: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1671: {
1678: PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
1679: return(0);
1680: }
1682: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1683: {
1684: Mat_IS *matis = (Mat_IS*)(A->data);
1687: matis->locempty = fix;
1688: return(0);
1689: }
1691: /*@
1692: MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.
1694: Collective on MPI_Comm
1696: Input Parameters:
1697: + B - the matrix
1698: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1699: (same value is used for all local rows)
1700: . d_nnz - array containing the number of nonzeros in the various rows of the
1701: DIAGONAL portion of the local submatrix (possibly different for each row)
1702: or NULL, if d_nz is used to specify the nonzero structure.
1703: The size of this array is equal to the number of local rows, i.e 'm'.
1704: For matrices that will be factored, you must leave room for (and set)
1705: the diagonal entry even if it is zero.
1706: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1707: submatrix (same value is used for all local rows).
1708: - o_nnz - array containing the number of nonzeros in the various rows of the
1709: OFF-DIAGONAL portion of the local submatrix (possibly different for
1710: each row) or NULL, if o_nz is used to specify the nonzero
1711: structure. The size of this array is equal to the number
1712: of local rows, i.e 'm'.
1714: If the *_nnz parameter is given then the *_nz parameter is ignored
1716: Level: intermediate
1718: Notes:
1719: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
1720: from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
1721: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1723: .keywords: matrix
1725: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1726: @*/
1727: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1728: {
1734: PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1735: return(0);
1736: }
1738: /* this is used by DMDA */
1739: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1740: {
1741: Mat_IS *matis = (Mat_IS*)(B->data);
1742: PetscInt bs,i,nlocalcols;
1746: if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
1748: if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
1749: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
1751: if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz;
1752: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
1754: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1755: MatGetSize(matis->A,NULL,&nlocalcols);
1756: MatGetBlockSize(matis->A,&bs);
1757: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1759: for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1760: MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1761: #if defined(PETSC_HAVE_HYPRE)
1762: MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1763: #endif
1765: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
1766: MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1768: nlocalcols /= bs;
1769: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols - i);
1770: MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1772: /* for other matrix types */
1773: MatSetUp(matis->A);
1774: return(0);
1775: }
1777: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1778: {
1779: Mat_IS *matis = (Mat_IS*)(A->data);
1780: PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1781: const PetscInt *global_indices_r,*global_indices_c;
1782: PetscInt i,j,bs,rows,cols;
1783: PetscInt lrows,lcols;
1784: PetscInt local_rows,local_cols;
1785: PetscMPIInt size;
1786: PetscBool isdense,issbaij;
1787: PetscErrorCode ierr;
1790: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1791: MatGetSize(A,&rows,&cols);
1792: MatGetBlockSize(A,&bs);
1793: MatGetSize(matis->A,&local_rows,&local_cols);
1794: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1795: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1796: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
1797: if (A->rmap->mapping != A->cmap->mapping) {
1798: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
1799: } else {
1800: global_indices_c = global_indices_r;
1801: }
1803: if (issbaij) {
1804: MatGetRowUpperTriangular(matis->A);
1805: }
1806: /*
1807: An SF reduce is needed to sum up properly on shared rows.
1808: Note that generally preallocation is not exact, since it overestimates nonzeros
1809: */
1810: MatGetLocalSize(A,&lrows,&lcols);
1811: MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1812: /* All processes need to compute entire row ownership */
1813: PetscMalloc1(rows,&row_ownership);
1814: MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1815: for (i=0;i<size;i++) {
1816: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
1817: row_ownership[j] = i;
1818: }
1819: }
1820: MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);
1822: /*
1823: my_dnz and my_onz contains exact contribution to preallocation from each local mat
1824: then, they will be summed up properly. This way, preallocation is always sufficient
1825: */
1826: PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1827: /* preallocation as a MATAIJ */
1828: if (isdense) { /* special case for dense local matrices */
1829: for (i=0;i<local_rows;i++) {
1830: PetscInt owner = row_ownership[global_indices_r[i]];
1831: for (j=0;j<local_cols;j++) {
1832: PetscInt index_col = global_indices_c[j];
1833: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1834: my_dnz[i] += 1;
1835: } else { /* offdiag block */
1836: my_onz[i] += 1;
1837: }
1838: }
1839: }
1840: } else if (matis->A->ops->getrowij) {
1841: const PetscInt *ii,*jj,*jptr;
1842: PetscBool done;
1843: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1844: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1845: jptr = jj;
1846: for (i=0;i<local_rows;i++) {
1847: PetscInt index_row = global_indices_r[i];
1848: for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1849: PetscInt owner = row_ownership[index_row];
1850: PetscInt index_col = global_indices_c[*jptr];
1851: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1852: my_dnz[i] += 1;
1853: } else { /* offdiag block */
1854: my_onz[i] += 1;
1855: }
1856: /* same as before, interchanging rows and cols */
1857: if (issbaij && index_col != index_row) {
1858: owner = row_ownership[index_col];
1859: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1860: my_dnz[*jptr] += 1;
1861: } else {
1862: my_onz[*jptr] += 1;
1863: }
1864: }
1865: }
1866: }
1867: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1868: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1869: } else { /* loop over rows and use MatGetRow */
1870: for (i=0;i<local_rows;i++) {
1871: const PetscInt *cols;
1872: PetscInt ncols,index_row = global_indices_r[i];
1873: MatGetRow(matis->A,i,&ncols,&cols,NULL);
1874: for (j=0;j<ncols;j++) {
1875: PetscInt owner = row_ownership[index_row];
1876: PetscInt index_col = global_indices_c[cols[j]];
1877: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1878: my_dnz[i] += 1;
1879: } else { /* offdiag block */
1880: my_onz[i] += 1;
1881: }
1882: /* same as before, interchanging rows and cols */
1883: if (issbaij && index_col != index_row) {
1884: owner = row_ownership[index_col];
1885: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1886: my_dnz[cols[j]] += 1;
1887: } else {
1888: my_onz[cols[j]] += 1;
1889: }
1890: }
1891: }
1892: MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1893: }
1894: }
1895: if (global_indices_c != global_indices_r) {
1896: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1897: }
1898: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1899: PetscFree(row_ownership);
1901: /* Reduce my_dnz and my_onz */
1902: if (maxreduce) {
1903: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1904: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1905: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1906: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1907: } else {
1908: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1909: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1910: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1911: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1912: }
1913: PetscFree2(my_dnz,my_onz);
1915: /* Resize preallocation if overestimated */
1916: for (i=0;i<lrows;i++) {
1917: dnz[i] = PetscMin(dnz[i],lcols);
1918: onz[i] = PetscMin(onz[i],cols-lcols);
1919: }
1921: /* Set preallocation */
1922: MatSeqAIJSetPreallocation(B,0,dnz);
1923: MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1924: for (i=0;i<lrows;i+=bs) {
1925: PetscInt b, d = dnz[i],o = onz[i];
1927: for (b=1;b<bs;b++) {
1928: d = PetscMax(d,dnz[i+b]);
1929: o = PetscMax(o,onz[i+b]);
1930: }
1931: dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
1932: onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
1933: }
1934: MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1935: MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1936: MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1937: MatPreallocateFinalize(dnz,onz);
1938: if (issbaij) {
1939: MatRestoreRowUpperTriangular(matis->A);
1940: }
1941: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1942: return(0);
1943: }
1945: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1946: {
1947: Mat_IS *matis = (Mat_IS*)(mat->data);
1948: Mat local_mat,MT;
1949: PetscInt rbs,cbs,rows,cols,lrows,lcols;
1950: PetscInt local_rows,local_cols;
1951: PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij;
1952: #if defined (PETSC_USE_DEBUG)
1953: PetscBool lb[4],bb[4];
1954: #endif
1955: PetscMPIInt size;
1956: PetscScalar *array;
1960: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
1961: if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1962: Mat B;
1963: IS irows = NULL,icols = NULL;
1964: PetscInt rbs,cbs;
1966: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
1967: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
1968: if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1969: IS rows,cols;
1970: const PetscInt *ridxs,*cidxs;
1971: PetscInt i,nw,*work;
1973: ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);
1974: ISLocalToGlobalMappingGetSize(mat->rmap->mapping,&nw);
1975: nw = nw/rbs;
1976: PetscCalloc1(nw,&work);
1977: for (i=0;i<nw;i++) work[ridxs[i]] += 1;
1978: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
1979: if (i == nw) {
1980: ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);
1981: ISSetPermutation(rows);
1982: ISInvertPermutation(rows,PETSC_DECIDE,&irows);
1983: ISDestroy(&rows);
1984: }
1985: ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);
1986: PetscFree(work);
1987: if (irows && mat->rmap->mapping != mat->cmap->mapping) {
1988: ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);
1989: ISLocalToGlobalMappingGetSize(mat->cmap->mapping,&nw);
1990: nw = nw/cbs;
1991: PetscCalloc1(nw,&work);
1992: for (i=0;i<nw;i++) work[cidxs[i]] += 1;
1993: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
1994: if (i == nw) {
1995: ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);
1996: ISSetPermutation(cols);
1997: ISInvertPermutation(cols,PETSC_DECIDE,&icols);
1998: ISDestroy(&cols);
1999: }
2000: ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);
2001: PetscFree(work);
2002: } else if (irows) {
2003: PetscObjectReference((PetscObject)irows);
2004: icols = irows;
2005: }
2006: } else {
2007: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);
2008: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);
2009: if (irows) { PetscObjectReference((PetscObject)irows); }
2010: if (icols) { PetscObjectReference((PetscObject)icols); }
2011: }
2012: if (!irows || !icols) {
2013: ISDestroy(&icols);
2014: ISDestroy(&irows);
2015: goto general_assembly;
2016: }
2017: MatConvert(matis->A,(rbs == cbs && rbs > 1) ? MATSEQBAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
2018: if (reuse != MAT_INPLACE_MATRIX) {
2019: MatCreateSubMatrix(B,irows,icols,reuse,M);
2020: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);
2021: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);
2022: } else {
2023: Mat C;
2025: MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
2026: MatHeaderReplace(mat,&C);
2027: }
2028: MatDestroy(&B);
2029: ISDestroy(&icols);
2030: ISDestroy(&irows);
2031: return(0);
2032: }
2033: general_assembly:
2034: MatGetSize(mat,&rows,&cols);
2035: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2036: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2037: MatGetLocalSize(mat,&lrows,&lcols);
2038: MatGetSize(matis->A,&local_rows,&local_cols);
2039: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
2040: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
2041: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
2042: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
2043: if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
2044: #if defined (PETSC_USE_DEBUG)
2045: lb[0] = isseqdense;
2046: lb[1] = isseqaij;
2047: lb[2] = isseqbaij;
2048: lb[3] = isseqsbaij;
2049: MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
2050: if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
2051: #endif
2053: if (reuse != MAT_REUSE_MATRIX) {
2054: MatCreate(PetscObjectComm((PetscObject)mat),&MT);
2055: MatSetSizes(MT,lrows,lcols,rows,cols);
2056: MatSetType(MT,mtype);
2057: MatSetBlockSizes(MT,rbs,cbs);
2058: MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
2059: } else {
2060: PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;
2062: /* some checks */
2063: MT = *M;
2064: MatGetBlockSizes(MT,&mrbs,&mcbs);
2065: MatGetSize(MT,&mrows,&mcols);
2066: MatGetLocalSize(MT,&mlrows,&mlcols);
2067: if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
2068: if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
2069: if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
2070: if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
2071: if (mrbs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%d != %d)",rbs,mrbs);
2072: if (mcbs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%d != %d)",cbs,mcbs);
2073: MatZeroEntries(MT);
2074: }
2076: if (isseqsbaij || isseqbaij) {
2077: MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
2078: isseqaij = PETSC_TRUE;
2079: } else {
2080: PetscObjectReference((PetscObject)matis->A);
2081: local_mat = matis->A;
2082: }
2084: /* Set values */
2085: MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);
2086: if (isseqdense) { /* special case for dense local matrices */
2087: PetscInt i,*dummy;
2089: PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
2090: for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
2091: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
2092: MatDenseGetArray(local_mat,&array);
2093: MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
2094: MatDenseRestoreArray(local_mat,&array);
2095: PetscFree(dummy);
2096: } else if (isseqaij) {
2097: const PetscInt *blocks;
2098: PetscInt i,nvtxs,*xadj,*adjncy, nb;
2099: PetscBool done;
2101: MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2102: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
2103: MatSeqAIJGetArray(local_mat,&array);
2104: MatGetVariableBlockSizes(local_mat,&nb,&blocks);
2105: if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2106: PetscInt sum;
2108: for (i=0,sum=0;i<nb;i++) sum += blocks[i];
2109: if (sum == nvtxs) {
2110: PetscInt r;
2112: for (i=0,r=0;i<nb;i++) {
2113: if (blocks[i] != xadj[r+1] - xadj[r]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %D: expected %D, got %D",i,blocks[i],xadj[r+1] - xadj[r]);
2114: MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],array+xadj[r],ADD_VALUES);
2115: r += blocks[i];
2116: }
2117: } else {
2118: for (i=0;i<nvtxs;i++) {
2119: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
2120: }
2121: }
2122: } else {
2123: for (i=0;i<nvtxs;i++) {
2124: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
2125: }
2126: }
2127: MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2128: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
2129: MatSeqAIJRestoreArray(local_mat,&array);
2130: } else { /* very basic values insertion for all other matrix types */
2131: PetscInt i;
2133: for (i=0;i<local_rows;i++) {
2134: PetscInt j;
2135: const PetscInt *local_indices_cols;
2137: MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
2138: MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
2139: MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
2140: }
2141: }
2142: MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
2143: MatDestroy(&local_mat);
2144: MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
2145: if (isseqdense) {
2146: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
2147: }
2148: if (reuse == MAT_INPLACE_MATRIX) {
2149: MatHeaderReplace(mat,&MT);
2150: } else if (reuse == MAT_INITIAL_MATRIX) {
2151: *M = MT;
2152: }
2153: return(0);
2154: }
2156: /*@
2157: MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format
2159: Input Parameter:
2160: . mat - the matrix (should be of type MATIS)
2161: . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
2163: Output Parameter:
2164: . newmat - the matrix in AIJ format
2166: Level: developer
2168: Notes:
2169: This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface.
2171: .seealso: MATIS, MatConvert()
2172: @*/
2173: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2174: {
2181: if (reuse == MAT_REUSE_MATRIX) {
2184: if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
2185: }
2186: PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
2187: return(0);
2188: }
2190: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
2191: {
2193: Mat_IS *matis = (Mat_IS*)(mat->data);
2194: PetscInt rbs,cbs,m,n,M,N;
2195: Mat B,localmat;
2198: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2199: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2200: MatGetSize(mat,&M,&N);
2201: MatGetLocalSize(mat,&m,&n);
2202: MatCreate(PetscObjectComm((PetscObject)mat),&B);
2203: MatSetSizes(B,m,n,M,N);
2204: MatSetBlockSize(B,rbs == cbs ? rbs : 1);
2205: MatSetType(B,MATIS);
2206: MatISSetLocalMatType(B,matis->lmattype);
2207: MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
2208: MatDuplicate(matis->A,op,&localmat);
2209: MatISSetLocalMat(B,localmat);
2210: MatDestroy(&localmat);
2211: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2212: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2213: *newmat = B;
2214: return(0);
2215: }
2217: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
2218: {
2220: Mat_IS *matis = (Mat_IS*)A->data;
2221: PetscBool local_sym;
2224: MatIsHermitian(matis->A,tol,&local_sym);
2225: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2226: return(0);
2227: }
2229: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
2230: {
2232: Mat_IS *matis = (Mat_IS*)A->data;
2233: PetscBool local_sym;
2236: MatIsSymmetric(matis->A,tol,&local_sym);
2237: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2238: return(0);
2239: }
2241: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg)
2242: {
2244: Mat_IS *matis = (Mat_IS*)A->data;
2245: PetscBool local_sym;
2248: if (A->rmap->mapping != A->cmap->mapping) {
2249: *flg = PETSC_FALSE;
2250: return(0);
2251: }
2252: MatIsStructurallySymmetric(matis->A,&local_sym);
2253: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2254: return(0);
2255: }
2257: static PetscErrorCode MatDestroy_IS(Mat A)
2258: {
2260: Mat_IS *b = (Mat_IS*)A->data;
2263: PetscFree(b->bdiag);
2264: PetscFree(b->lmattype);
2265: MatDestroy(&b->A);
2266: VecScatterDestroy(&b->cctx);
2267: VecScatterDestroy(&b->rctx);
2268: VecDestroy(&b->x);
2269: VecDestroy(&b->y);
2270: VecDestroy(&b->counter);
2271: ISDestroy(&b->getsub_ris);
2272: ISDestroy(&b->getsub_cis);
2273: if (b->sf != b->csf) {
2274: PetscSFDestroy(&b->csf);
2275: PetscFree2(b->csf_rootdata,b->csf_leafdata);
2276: } else b->csf = NULL;
2277: PetscSFDestroy(&b->sf);
2278: PetscFree2(b->sf_rootdata,b->sf_leafdata);
2279: PetscFree(A->data);
2280: PetscObjectChangeTypeName((PetscObject)A,0);
2281: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
2282: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
2283: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
2284: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
2285: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
2286: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
2287: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
2288: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
2289: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
2290: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
2291: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
2292: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
2293: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
2294: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
2295: return(0);
2296: }
2298: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
2299: {
2301: Mat_IS *is = (Mat_IS*)A->data;
2302: PetscScalar zero = 0.0;
2305: /* scatter the global vector x into the local work vector */
2306: VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2307: VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2309: /* multiply the local matrix */
2310: MatMult(is->A,is->x,is->y);
2312: /* scatter product back into global memory */
2313: VecSet(y,zero);
2314: VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2315: VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2316: return(0);
2317: }
2319: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2320: {
2321: Vec temp_vec;
2325: if (v3 != v2) {
2326: MatMult(A,v1,v3);
2327: VecAXPY(v3,1.0,v2);
2328: } else {
2329: VecDuplicate(v2,&temp_vec);
2330: MatMult(A,v1,temp_vec);
2331: VecAXPY(temp_vec,1.0,v2);
2332: VecCopy(temp_vec,v3);
2333: VecDestroy(&temp_vec);
2334: }
2335: return(0);
2336: }
2338: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
2339: {
2340: Mat_IS *is = (Mat_IS*)A->data;
2344: /* scatter the global vector x into the local work vector */
2345: VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2346: VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2348: /* multiply the local matrix */
2349: MatMultTranspose(is->A,is->y,is->x);
2351: /* scatter product back into global vector */
2352: VecSet(x,0);
2353: VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2354: VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2355: return(0);
2356: }
2358: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2359: {
2360: Vec temp_vec;
2364: if (v3 != v2) {
2365: MatMultTranspose(A,v1,v3);
2366: VecAXPY(v3,1.0,v2);
2367: } else {
2368: VecDuplicate(v2,&temp_vec);
2369: MatMultTranspose(A,v1,temp_vec);
2370: VecAXPY(temp_vec,1.0,v2);
2371: VecCopy(temp_vec,v3);
2372: VecDestroy(&temp_vec);
2373: }
2374: return(0);
2375: }
2377: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
2378: {
2379: Mat_IS *a = (Mat_IS*)A->data;
2381: PetscViewer sviewer;
2382: PetscBool isascii,view = PETSC_TRUE;
2385: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2386: if (isascii) {
2387: PetscViewerFormat format;
2389: PetscViewerGetFormat(viewer,&format);
2390: if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2391: }
2392: if (!view) return(0);
2393: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2394: MatView(a->A,sviewer);
2395: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2396: PetscViewerFlush(viewer);
2397: return(0);
2398: }
2400: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values)
2401: {
2402: Mat_IS *is = (Mat_IS*)mat->data;
2403: MPI_Datatype nodeType;
2404: const PetscScalar *lv;
2405: PetscInt bs;
2406: PetscErrorCode ierr;
2409: MatGetBlockSize(mat,&bs);
2410: MatSetBlockSize(is->A,bs);
2411: MatInvertBlockDiagonal(is->A,&lv);
2412: if (!is->bdiag) {
2413: PetscMalloc1(bs*mat->rmap->n,&is->bdiag);
2414: }
2415: MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType);
2416: MPI_Type_commit(&nodeType);
2417: PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPIU_REPLACE);
2418: PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPIU_REPLACE);
2419: MPI_Type_free(&nodeType);
2420: if (values) *values = is->bdiag;
2421: return(0);
2422: }
2424: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2425: {
2426: Vec cglobal,rglobal;
2427: IS from;
2428: Mat_IS *is = (Mat_IS*)A->data;
2429: PetscScalar sum;
2430: const PetscInt *garray;
2431: PetscInt nr,rbs,nc,cbs;
2432: PetscBool iscuda;
2436: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);
2437: ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);
2438: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);
2439: ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);
2440: VecDestroy(&is->x);
2441: VecDestroy(&is->y);
2442: VecDestroy(&is->counter);
2443: VecScatterDestroy(&is->rctx);
2444: VecScatterDestroy(&is->cctx);
2445: MatCreateVecs(is->A,&is->x,&is->y);
2446: PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);
2447: if (iscuda) {
2448: PetscFree(A->defaultvectype);
2449: PetscStrallocpy(VECCUDA,&A->defaultvectype);
2450: }
2451: MatCreateVecs(A,&cglobal,&rglobal);
2452: ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);
2453: ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
2454: VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
2455: ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);
2456: ISDestroy(&from);
2457: if (A->rmap->mapping != A->cmap->mapping) {
2458: ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);
2459: ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
2460: VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
2461: ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);
2462: ISDestroy(&from);
2463: } else {
2464: PetscObjectReference((PetscObject)is->rctx);
2465: is->cctx = is->rctx;
2466: }
2467: VecDestroy(&cglobal);
2469: /* interface counter vector (local) */
2470: VecDuplicate(is->y,&is->counter);
2471: VecSet(is->y,1.);
2472: VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2473: VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2474: VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2475: VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2477: /* special functions for block-diagonal matrices */
2478: VecSum(rglobal,&sum);
2479: if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && A->rmap->mapping == A->cmap->mapping) {
2480: A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2481: } else {
2482: A->ops->invertblockdiagonal = NULL;
2483: }
2484: VecDestroy(&rglobal);
2486: /* setup SF for general purpose shared indices based communications */
2487: MatISSetUpSF_IS(A);
2488: return(0);
2489: }
2491: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
2492: {
2494: PetscInt nr,rbs,nc,cbs;
2495: Mat_IS *is = (Mat_IS*)A->data;
2500: MatDestroy(&is->A);
2501: if (is->csf != is->sf) {
2502: PetscSFDestroy(&is->csf);
2503: PetscFree2(is->csf_rootdata,is->csf_leafdata);
2504: } else is->csf = NULL;
2505: PetscSFDestroy(&is->sf);
2506: PetscFree2(is->sf_rootdata,is->sf_leafdata);
2507: PetscFree(is->bdiag);
2509: /* Setup Layout and set local to global maps */
2510: PetscLayoutSetUp(A->rmap);
2511: PetscLayoutSetUp(A->cmap);
2512: ISLocalToGlobalMappingGetSize(rmapping,&nr);
2513: ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
2514: ISLocalToGlobalMappingGetSize(cmapping,&nc);
2515: ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
2516: /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */
2517: if (rmapping != cmapping && A->rmap->N == A->cmap->N) {
2518: PetscBool same,gsame;
2520: same = PETSC_FALSE;
2521: if (nr == nc && cbs == rbs) {
2522: const PetscInt *idxs1,*idxs2;
2524: ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);
2525: ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);
2526: PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);
2527: ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);
2528: ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);
2529: }
2530: MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2531: if (gsame) cmapping = rmapping;
2532: }
2533: PetscLayoutSetBlockSize(A->rmap,rbs);
2534: PetscLayoutSetBlockSize(A->cmap,cbs);
2535: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
2536: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
2538: /* Create the local matrix A */
2539: MatCreate(PETSC_COMM_SELF,&is->A);
2540: MatSetType(is->A,is->lmattype);
2541: MatSetSizes(is->A,nr,nc,nr,nc);
2542: MatSetBlockSizes(is->A,rbs,cbs);
2543: MatSetOptionsPrefix(is->A,"is_");
2544: MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);
2545: PetscLayoutSetUp(is->A->rmap);
2546: PetscLayoutSetUp(is->A->cmap);
2548: if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
2549: MatISSetUpScatters_Private(A);
2550: }
2551: MatSetUp(A);
2552: return(0);
2553: }
2555: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2556: {
2557: Mat_IS *is = (Mat_IS*)mat->data;
2559: #if defined(PETSC_USE_DEBUG)
2560: PetscInt i,zm,zn;
2561: #endif
2562: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2565: #if defined(PETSC_USE_DEBUG)
2566: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
2567: /* count negative indices */
2568: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2569: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2570: #endif
2571: ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2572: ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2573: #if defined(PETSC_USE_DEBUG)
2574: /* count negative indices (should be the same as before) */
2575: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2576: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2577: if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
2578: if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
2579: #endif
2580: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
2581: return(0);
2582: }
2584: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2585: {
2586: Mat_IS *is = (Mat_IS*)mat->data;
2588: #if defined(PETSC_USE_DEBUG)
2589: PetscInt i,zm,zn;
2590: #endif
2591: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2594: #if defined(PETSC_USE_DEBUG)
2595: if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
2596: /* count negative indices */
2597: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2598: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2599: #endif
2600: ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2601: ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2602: #if defined(PETSC_USE_DEBUG)
2603: /* count negative indices (should be the same as before) */
2604: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2605: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2606: if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
2607: if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
2608: #endif
2609: MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
2610: return(0);
2611: }
2613: static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2614: {
2616: Mat_IS *is = (Mat_IS*)A->data;
2619: if (is->A->rmap->mapping) {
2620: MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
2621: } else {
2622: MatSetValues(is->A,m,rows,n,cols,values,addv);
2623: }
2624: return(0);
2625: }
2627: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2628: {
2630: Mat_IS *is = (Mat_IS*)A->data;
2633: if (is->A->rmap->mapping) {
2634: #if defined(PETSC_USE_DEBUG)
2635: PetscInt ibs,bs;
2637: ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);
2638: MatGetBlockSize(is->A,&bs);
2639: if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs);
2640: #endif
2641: MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
2642: } else {
2643: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
2644: }
2645: return(0);
2646: }
2648: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns)
2649: {
2650: Mat_IS *is = (Mat_IS*)A->data;
2654: if (!n) {
2655: is->pure_neumann = PETSC_TRUE;
2656: } else {
2657: PetscInt i;
2658: is->pure_neumann = PETSC_FALSE;
2660: if (columns) {
2661: MatZeroRowsColumns(is->A,n,rows,diag,0,0);
2662: } else {
2663: MatZeroRows(is->A,n,rows,diag,0,0);
2664: }
2665: if (diag != 0.) {
2666: const PetscScalar *array;
2667: VecGetArrayRead(is->counter,&array);
2668: for (i=0; i<n; i++) {
2669: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
2670: }
2671: VecRestoreArrayRead(is->counter,&array);
2672: }
2673: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
2674: MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
2675: }
2676: return(0);
2677: }
2679: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
2680: {
2681: Mat_IS *matis = (Mat_IS*)A->data;
2682: PetscInt nr,nl,len,i;
2683: PetscInt *lrows;
2687: #if defined(PETSC_USE_DEBUG)
2688: if (columns || diag != 0. || (x && b)) {
2689: PetscBool cong;
2691: PetscLayoutCompare(A->rmap,A->cmap,&cong);
2692: cong = (PetscBool)(cong && matis->sf == matis->csf);
2693: if (!cong && columns) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2694: if (!cong && diag != 0.) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2695: if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
2696: }
2697: #endif
2698: /* get locally owned rows */
2699: PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);
2700: /* fix right hand side if needed */
2701: if (x && b) {
2702: const PetscScalar *xx;
2703: PetscScalar *bb;
2705: VecGetArrayRead(x, &xx);
2706: VecGetArray(b, &bb);
2707: for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
2708: VecRestoreArrayRead(x, &xx);
2709: VecRestoreArray(b, &bb);
2710: }
2711: /* get rows associated to the local matrices */
2712: MatGetSize(matis->A,&nl,NULL);
2713: PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));
2714: PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));
2715: for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
2716: PetscFree(lrows);
2717: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
2718: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
2719: PetscMalloc1(nl,&lrows);
2720: for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
2721: MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
2722: PetscFree(lrows);
2723: return(0);
2724: }
2726: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2727: {
2731: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
2732: return(0);
2733: }
2735: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2736: {
2740: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
2741: return(0);
2742: }
2744: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
2745: {
2746: Mat_IS *is = (Mat_IS*)A->data;
2750: MatAssemblyBegin(is->A,type);
2751: return(0);
2752: }
2754: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
2755: {
2756: Mat_IS *is = (Mat_IS*)A->data;
2760: MatAssemblyEnd(is->A,type);
2761: /* fix for local empty rows/cols */
2762: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2763: Mat newlA;
2764: ISLocalToGlobalMapping rl2g,cl2g;
2765: IS nzr,nzc;
2766: PetscInt nr,nc,nnzr,nnzc;
2767: PetscBool lnewl2g,newl2g;
2769: MatGetSize(is->A,&nr,&nc);
2770: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
2771: if (!nzr) {
2772: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
2773: }
2774: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
2775: if (!nzc) {
2776: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
2777: }
2778: ISGetSize(nzr,&nnzr);
2779: ISGetSize(nzc,&nnzc);
2780: if (nnzr != nr || nnzc != nc) {
2781: ISLocalToGlobalMapping l2g;
2782: IS is1,is2;
2784: /* need new global l2g map */
2785: lnewl2g = PETSC_TRUE;
2786: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2788: /* extract valid submatrix */
2789: MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);
2791: /* attach local l2g maps for successive calls of MatSetValues on the local matrix */
2792: ISLocalToGlobalMappingCreateIS(nzr,&l2g);
2793: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);
2794: ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2795: ISLocalToGlobalMappingDestroy(&l2g);
2796: if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2797: const PetscInt *idxs1,*idxs2;
2798: PetscInt j,i,nl,*tidxs;
2800: ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);
2801: ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);
2802: PetscMalloc1(nl,&tidxs);
2803: ISGetIndices(is2,&idxs2);
2804: for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2805: if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr);
2806: ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);
2807: ISRestoreIndices(is2,&idxs2);
2808: ISDestroy(&is2);
2809: ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2810: }
2811: ISLocalToGlobalMappingCreateIS(is2,&rl2g);
2812: ISDestroy(&is1);
2813: ISDestroy(&is2);
2815: ISLocalToGlobalMappingCreateIS(nzc,&l2g);
2816: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);
2817: ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2818: ISLocalToGlobalMappingDestroy(&l2g);
2819: if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2820: const PetscInt *idxs1,*idxs2;
2821: PetscInt j,i,nl,*tidxs;
2823: ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);
2824: ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);
2825: PetscMalloc1(nl,&tidxs);
2826: ISGetIndices(is2,&idxs2);
2827: for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2828: if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc);
2829: ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);
2830: ISRestoreIndices(is2,&idxs2);
2831: ISDestroy(&is2);
2832: ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2833: }
2834: ISLocalToGlobalMappingCreateIS(is2,&cl2g);
2835: ISDestroy(&is1);
2836: ISDestroy(&is2);
2838: MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);
2840: ISLocalToGlobalMappingDestroy(&rl2g);
2841: ISLocalToGlobalMappingDestroy(&cl2g);
2842: } else { /* local matrix fully populated */
2843: lnewl2g = PETSC_FALSE;
2844: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2845: PetscObjectReference((PetscObject)is->A);
2846: newlA = is->A;
2847: }
2848: /* attach new global l2g map if needed */
2849: if (newl2g) {
2850: IS gnzr,gnzc;
2851: const PetscInt *grid,*gcid;
2853: ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);
2854: ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);
2855: ISGetIndices(gnzr,&grid);
2856: ISGetIndices(gnzc,&gcid);
2857: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);
2858: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);
2859: ISRestoreIndices(gnzr,&grid);
2860: ISRestoreIndices(gnzc,&gcid);
2861: ISDestroy(&gnzr);
2862: ISDestroy(&gnzc);
2863: MatSetLocalToGlobalMapping(A,rl2g,cl2g);
2864: ISLocalToGlobalMappingDestroy(&rl2g);
2865: ISLocalToGlobalMappingDestroy(&cl2g);
2866: }
2867: MatISSetLocalMat(A,newlA);
2868: MatDestroy(&newlA);
2869: ISDestroy(&nzr);
2870: ISDestroy(&nzc);
2871: is->locempty = PETSC_FALSE;
2872: }
2873: return(0);
2874: }
2876: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
2877: {
2878: Mat_IS *is = (Mat_IS*)mat->data;
2881: *local = is->A;
2882: return(0);
2883: }
2885: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
2886: {
2888: *local = NULL;
2889: return(0);
2890: }
2892: /*@
2893: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
2895: Input Parameter:
2896: . mat - the matrix
2898: Output Parameter:
2899: . local - the local matrix
2901: Level: advanced
2903: Notes:
2904: This can be called if you have precomputed the nonzero structure of the
2905: matrix and want to provide it to the inner matrix object to improve the performance
2906: of the MatSetValues() operation.
2908: Call MatISRestoreLocalMat() when finished with the local matrix.
2910: .seealso: MATIS
2911: @*/
2912: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
2913: {
2919: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
2920: return(0);
2921: }
2923: /*@
2924: MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()
2926: Input Parameter:
2927: . mat - the matrix
2929: Output Parameter:
2930: . local - the local matrix
2932: Level: advanced
2934: .seealso: MATIS
2935: @*/
2936: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
2937: {
2943: PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
2944: return(0);
2945: }
2947: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
2948: {
2949: Mat_IS *is = (Mat_IS*)mat->data;
2953: if (is->A) {
2954: MatSetType(is->A,mtype);
2955: }
2956: PetscFree(is->lmattype);
2957: PetscStrallocpy(mtype,&is->lmattype);
2958: return(0);
2959: }
2961: /*@
2962: MatISSetLocalMatType - Specifies the type of local matrix
2964: Input Parameter:
2965: . mat - the matrix
2966: . mtype - the local matrix type
2968: Output Parameter:
2970: Level: advanced
2972: .seealso: MATIS, MatSetType(), MatType
2973: @*/
2974: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
2975: {
2980: PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
2981: return(0);
2982: }
2984: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
2985: {
2986: Mat_IS *is = (Mat_IS*)mat->data;
2987: PetscInt nrows,ncols,orows,ocols;
2989: MatType mtype,otype;
2990: PetscBool sametype = PETSC_TRUE;
2993: if (is->A) {
2994: MatGetSize(is->A,&orows,&ocols);
2995: MatGetSize(local,&nrows,&ncols);
2996: if (orows != nrows || ocols != ncols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local MATIS matrix should be of size %Dx%D (you passed a %Dx%D matrix)",orows,ocols,nrows,ncols);
2997: MatGetType(local,&mtype);
2998: MatGetType(is->A,&otype);
2999: PetscStrcmp(mtype,otype,&sametype);
3000: }
3001: PetscObjectReference((PetscObject)local);
3002: MatDestroy(&is->A);
3003: is->A = local;
3004: MatGetType(is->A,&mtype);
3005: MatISSetLocalMatType(mat,mtype);
3006: if (!sametype && !is->islocalref) {
3007: MatISSetUpScatters_Private(mat);
3008: }
3009: return(0);
3010: }
3012: /*@
3013: MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.
3015: Collective on Mat
3017: Input Parameter:
3018: . mat - the matrix
3019: . local - the local matrix
3021: Output Parameter:
3023: Level: advanced
3025: Notes:
3026: This can be called if you have precomputed the local matrix and
3027: want to provide it to the matrix object MATIS.
3029: .seealso: MATIS
3030: @*/
3031: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
3032: {
3038: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
3039: return(0);
3040: }
3042: static PetscErrorCode MatZeroEntries_IS(Mat A)
3043: {
3044: Mat_IS *a = (Mat_IS*)A->data;
3048: MatZeroEntries(a->A);
3049: return(0);
3050: }
3052: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
3053: {
3054: Mat_IS *is = (Mat_IS*)A->data;
3058: MatScale(is->A,a);
3059: return(0);
3060: }
3062: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3063: {
3064: Mat_IS *is = (Mat_IS*)A->data;
3068: /* get diagonal of the local matrix */
3069: MatGetDiagonal(is->A,is->y);
3071: /* scatter diagonal back into global vector */
3072: VecSet(v,0);
3073: VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3074: VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3075: return(0);
3076: }
3078: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
3079: {
3080: Mat_IS *a = (Mat_IS*)A->data;
3084: MatSetOption(a->A,op,flg);
3085: return(0);
3086: }
3088: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
3089: {
3090: Mat_IS *y = (Mat_IS*)Y->data;
3091: Mat_IS *x;
3092: #if defined(PETSC_USE_DEBUG)
3093: PetscBool ismatis;
3094: #endif
3098: #if defined(PETSC_USE_DEBUG)
3099: PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
3100: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3101: #endif
3102: x = (Mat_IS*)X->data;
3103: MatAXPY(y->A,a,x->A,str);
3104: return(0);
3105: }
3107: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
3108: {
3109: Mat lA;
3110: Mat_IS *matis;
3111: ISLocalToGlobalMapping rl2g,cl2g;
3112: IS is;
3113: const PetscInt *rg,*rl;
3114: PetscInt nrg;
3115: PetscInt N,M,nrl,i,*idxs;
3116: PetscErrorCode ierr;
3119: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
3120: ISGetLocalSize(row,&nrl);
3121: ISGetIndices(row,&rl);
3122: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
3123: #if defined(PETSC_USE_DEBUG)
3124: for (i=0;i<nrl;i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater then maximum possible %D",i,rl[i],nrg);
3125: #endif
3126: PetscMalloc1(nrg,&idxs);
3127: /* map from [0,nrl) to row */
3128: for (i=0;i<nrl;i++) idxs[i] = rl[i];
3129: for (i=nrl;i<nrg;i++) idxs[i] = -1;
3130: ISRestoreIndices(row,&rl);
3131: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
3132: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
3133: ISLocalToGlobalMappingCreateIS(is,&rl2g);
3134: ISDestroy(&is);
3135: /* compute new l2g map for columns */
3136: if (col != row || A->rmap->mapping != A->cmap->mapping) {
3137: const PetscInt *cg,*cl;
3138: PetscInt ncg;
3139: PetscInt ncl;
3141: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
3142: ISGetLocalSize(col,&ncl);
3143: ISGetIndices(col,&cl);
3144: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
3145: #if defined(PETSC_USE_DEBUG)
3146: for (i=0;i<ncl;i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater then maximum possible %D",i,cl[i],ncg);
3147: #endif
3148: PetscMalloc1(ncg,&idxs);
3149: /* map from [0,ncl) to col */
3150: for (i=0;i<ncl;i++) idxs[i] = cl[i];
3151: for (i=ncl;i<ncg;i++) idxs[i] = -1;
3152: ISRestoreIndices(col,&cl);
3153: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
3154: ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
3155: ISLocalToGlobalMappingCreateIS(is,&cl2g);
3156: ISDestroy(&is);
3157: } else {
3158: PetscObjectReference((PetscObject)rl2g);
3159: cl2g = rl2g;
3160: }
3161: /* create the MATIS submatrix */
3162: MatGetSize(A,&M,&N);
3163: MatCreate(PetscObjectComm((PetscObject)A),submat);
3164: MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
3165: MatSetType(*submat,MATIS);
3166: matis = (Mat_IS*)((*submat)->data);
3167: matis->islocalref = PETSC_TRUE;
3168: MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
3169: MatISGetLocalMat(A,&lA);
3170: MatISSetLocalMat(*submat,lA);
3171: MatSetUp(*submat);
3172: MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
3173: MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
3174: ISLocalToGlobalMappingDestroy(&rl2g);
3175: ISLocalToGlobalMappingDestroy(&cl2g);
3176: /* remove unsupported ops */
3177: PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
3178: (*submat)->ops->destroy = MatDestroy_IS;
3179: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3180: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3181: (*submat)->ops->assemblybegin = MatAssemblyBegin_IS;
3182: (*submat)->ops->assemblyend = MatAssemblyEnd_IS;
3183: return(0);
3184: }
3186: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
3187: {
3188: Mat_IS *a = (Mat_IS*)A->data;
3189: char type[256];
3190: PetscBool flg;
3194: PetscOptionsHead(PetscOptionsObject,"MATIS options");
3195: PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
3196: PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
3197: PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
3198: if (flg) {
3199: MatISSetLocalMatType(A,type);
3200: }
3201: if (a->A) {
3202: MatSetFromOptions(a->A);
3203: }
3204: PetscOptionsTail();
3205: return(0);
3206: }
3208: /*@
3209: MatCreateIS - Creates a "process" unassembled matrix, assembled on each
3210: process but not across processes.
3212: Input Parameters:
3213: + comm - MPI communicator that will share the matrix
3214: . bs - block size of the matrix
3215: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3216: . rmap - local to global map for rows
3217: - cmap - local to global map for cols
3219: Output Parameter:
3220: . A - the resulting matrix
3222: Level: advanced
3224: Notes:
3225: See MATIS for more details.
3226: m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors
3227: used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
3228: If either rmap or cmap are NULL, then the matrix is assumed to be square.
3230: .seealso: MATIS, MatSetLocalToGlobalMapping()
3231: @*/
3232: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
3233: {
3237: if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
3238: MatCreate(comm,A);
3239: MatSetSizes(*A,m,n,M,N);
3240: if (bs > 0) {
3241: MatSetBlockSize(*A,bs);
3242: }
3243: MatSetType(*A,MATIS);
3244: if (rmap && cmap) {
3245: MatSetLocalToGlobalMapping(*A,rmap,cmap);
3246: } else if (!rmap) {
3247: MatSetLocalToGlobalMapping(*A,cmap,cmap);
3248: } else {
3249: MatSetLocalToGlobalMapping(*A,rmap,rmap);
3250: }
3251: return(0);
3252: }
3254: /*MC
3255: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP).
3256: This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector
3257: product is handled "implicitly".
3259: Options Database Keys:
3260: + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
3261: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3262: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP().
3264: Notes:
3265: Options prefix for the inner matrix are given by -is_mat_xxx
3267: You must call MatSetLocalToGlobalMapping() before using this matrix type.
3269: You can do matrix preallocation on the local matrix after you obtain it with
3270: MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation()
3272: Level: advanced
3274: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP
3276: M*/
3278: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3279: {
3281: Mat_IS *b;
3284: PetscNewLog(A,&b);
3285: PetscStrallocpy(MATAIJ,&b->lmattype);
3286: A->data = (void*)b;
3288: /* matrix ops */
3289: PetscMemzero(A->ops,sizeof(struct _MatOps));
3290: A->ops->mult = MatMult_IS;
3291: A->ops->multadd = MatMultAdd_IS;
3292: A->ops->multtranspose = MatMultTranspose_IS;
3293: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3294: A->ops->destroy = MatDestroy_IS;
3295: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3296: A->ops->setvalues = MatSetValues_IS;
3297: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3298: A->ops->setvalueslocal = MatSetValuesLocal_IS;
3299: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3300: A->ops->zerorows = MatZeroRows_IS;
3301: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3302: A->ops->assemblybegin = MatAssemblyBegin_IS;
3303: A->ops->assemblyend = MatAssemblyEnd_IS;
3304: A->ops->view = MatView_IS;
3305: A->ops->zeroentries = MatZeroEntries_IS;
3306: A->ops->scale = MatScale_IS;
3307: A->ops->getdiagonal = MatGetDiagonal_IS;
3308: A->ops->setoption = MatSetOption_IS;
3309: A->ops->ishermitian = MatIsHermitian_IS;
3310: A->ops->issymmetric = MatIsSymmetric_IS;
3311: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3312: A->ops->duplicate = MatDuplicate_IS;
3313: A->ops->missingdiagonal = MatMissingDiagonal_IS;
3314: A->ops->copy = MatCopy_IS;
3315: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3316: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3317: A->ops->axpy = MatAXPY_IS;
3318: A->ops->diagonalset = MatDiagonalSet_IS;
3319: A->ops->shift = MatShift_IS;
3320: A->ops->transpose = MatTranspose_IS;
3321: A->ops->getinfo = MatGetInfo_IS;
3322: A->ops->diagonalscale = MatDiagonalScale_IS;
3323: A->ops->setfromoptions = MatSetFromOptions_IS;
3325: /* special MATIS functions */
3326: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
3327: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
3328: PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
3329: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
3330: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
3331: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
3332: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
3333: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
3334: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
3335: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
3336: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
3337: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
3338: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
3339: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
3340: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
3341: PetscObjectChangeTypeName((PetscObject)A,MATIS);
3342: return(0);
3343: }