Actual source code: matis.c
petsc-3.10.5 2019-03-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: /*@
803: MatISSetUpSF - Setup star forest objects used by MatIS.
805: Collective on Mat
807: Input Parameters:
808: + A - the matrix
810: Level: advanced
812: Notes:
813: This function does not need to be called by the user.
815: .keywords: matrix
817: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatISGetLocalMat()
818: @*/
819: PetscErrorCode MatISSetUpSF(Mat A)
820: {
826: PetscTryMethod(A,"MatISSetUpSF_C",(Mat),(A));
827: return(0);
828: }
830: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
831: {
832: Mat **nest,*snest,**rnest,lA,B;
833: IS *iscol,*isrow,*islrow,*islcol;
834: ISLocalToGlobalMapping rl2g,cl2g;
835: MPI_Comm comm;
836: PetscInt *lr,*lc,*l2gidxs;
837: PetscInt i,j,nr,nc,rbs,cbs;
838: PetscBool convert,lreuse,*istrans;
839: PetscErrorCode ierr;
842: MatNestGetSubMats(A,&nr,&nc,&nest);
843: lreuse = PETSC_FALSE;
844: rnest = NULL;
845: if (reuse == MAT_REUSE_MATRIX) {
846: PetscBool ismatis,isnest;
848: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
849: if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
850: MatISGetLocalMat(*newmat,&lA);
851: PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
852: if (isnest) {
853: MatNestGetSubMats(lA,&i,&j,&rnest);
854: lreuse = (PetscBool)(i == nr && j == nc);
855: if (!lreuse) rnest = NULL;
856: }
857: }
858: PetscObjectGetComm((PetscObject)A,&comm);
859: PetscCalloc2(nr,&lr,nc,&lc);
860: PetscCalloc6(nr,&isrow,nc,&iscol,
861: nr,&islrow,nc,&islcol,
862: nr*nc,&snest,nr*nc,&istrans);
863: MatNestGetISs(A,isrow,iscol);
864: for (i=0;i<nr;i++) {
865: for (j=0;j<nc;j++) {
866: PetscBool ismatis;
867: PetscInt l1,l2,lb1,lb2,ij=i*nc+j;
869: /* Null matrix pointers are allowed in MATNEST */
870: if (!nest[i][j]) continue;
872: /* Nested matrices should be of type MATIS */
873: PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
874: if (istrans[ij]) {
875: Mat T,lT;
876: MatTransposeGetMat(nest[i][j],&T);
877: PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
878: 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);
879: MatISGetLocalMat(T,&lT);
880: MatCreateTranspose(lT,&snest[ij]);
881: } else {
882: PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
883: if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
884: MatISGetLocalMat(nest[i][j],&snest[ij]);
885: }
887: /* Check compatibility of local sizes */
888: MatGetSize(snest[ij],&l1,&l2);
889: MatGetBlockSizes(snest[ij],&lb1,&lb2);
890: if (!l1 || !l2) continue;
891: 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);
892: 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);
893: lr[i] = l1;
894: lc[j] = l2;
896: /* check compatibilty for local matrix reusage */
897: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
898: }
899: }
901: #if defined (PETSC_USE_DEBUG)
902: /* Check compatibility of l2g maps for rows */
903: for (i=0;i<nr;i++) {
904: rl2g = NULL;
905: for (j=0;j<nc;j++) {
906: PetscInt n1,n2;
908: if (!nest[i][j]) continue;
909: if (istrans[i*nc+j]) {
910: Mat T;
912: MatTransposeGetMat(nest[i][j],&T);
913: MatGetLocalToGlobalMapping(T,NULL,&cl2g);
914: } else {
915: MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
916: }
917: ISLocalToGlobalMappingGetSize(cl2g,&n1);
918: if (!n1) continue;
919: if (!rl2g) {
920: rl2g = cl2g;
921: } else {
922: const PetscInt *idxs1,*idxs2;
923: PetscBool same;
925: ISLocalToGlobalMappingGetSize(rl2g,&n2);
926: 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);
927: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
928: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
929: PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
930: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
931: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
932: 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);
933: }
934: }
935: }
936: /* Check compatibility of l2g maps for columns */
937: for (i=0;i<nc;i++) {
938: rl2g = NULL;
939: for (j=0;j<nr;j++) {
940: PetscInt n1,n2;
942: if (!nest[j][i]) continue;
943: if (istrans[j*nc+i]) {
944: Mat T;
946: MatTransposeGetMat(nest[j][i],&T);
947: MatGetLocalToGlobalMapping(T,&cl2g,NULL);
948: } else {
949: MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
950: }
951: ISLocalToGlobalMappingGetSize(cl2g,&n1);
952: if (!n1) continue;
953: if (!rl2g) {
954: rl2g = cl2g;
955: } else {
956: const PetscInt *idxs1,*idxs2;
957: PetscBool same;
959: ISLocalToGlobalMappingGetSize(rl2g,&n2);
960: 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);
961: ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
962: ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
963: PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
964: ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
965: ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
966: 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);
967: }
968: }
969: }
970: #endif
972: B = NULL;
973: if (reuse != MAT_REUSE_MATRIX) {
974: PetscInt stl;
976: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
977: for (i=0,stl=0;i<nr;i++) stl += lr[i];
978: PetscMalloc1(stl,&l2gidxs);
979: for (i=0,stl=0;i<nr;i++) {
980: Mat usedmat;
981: Mat_IS *matis;
982: const PetscInt *idxs;
984: /* local IS for local NEST */
985: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
987: /* l2gmap */
988: j = 0;
989: usedmat = nest[i][j];
990: while (!usedmat && j < nc-1) usedmat = nest[i][++j];
991: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");
993: if (istrans[i*nc+j]) {
994: Mat T;
995: MatTransposeGetMat(usedmat,&T);
996: usedmat = T;
997: }
998: MatISSetUpSF(usedmat);
999: matis = (Mat_IS*)(usedmat->data);
1000: ISGetIndices(isrow[i],&idxs);
1001: if (istrans[i*nc+j]) {
1002: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
1003: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
1004: } else {
1005: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
1006: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
1007: }
1008: ISRestoreIndices(isrow[i],&idxs);
1009: stl += lr[i];
1010: }
1011: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);
1013: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
1014: for (i=0,stl=0;i<nc;i++) stl += lc[i];
1015: PetscMalloc1(stl,&l2gidxs);
1016: for (i=0,stl=0;i<nc;i++) {
1017: Mat usedmat;
1018: Mat_IS *matis;
1019: const PetscInt *idxs;
1021: /* local IS for local NEST */
1022: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
1024: /* l2gmap */
1025: j = 0;
1026: usedmat = nest[j][i];
1027: while (!usedmat && j < nr-1) usedmat = nest[++j][i];
1028: if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat");
1029: if (istrans[j*nc+i]) {
1030: Mat T;
1031: MatTransposeGetMat(usedmat,&T);
1032: usedmat = T;
1033: }
1034: MatISSetUpSF(usedmat);
1035: matis = (Mat_IS*)(usedmat->data);
1036: ISGetIndices(iscol[i],&idxs);
1037: if (istrans[j*nc+i]) {
1038: PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
1039: PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
1040: } else {
1041: PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
1042: PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
1043: }
1044: ISRestoreIndices(iscol[i],&idxs);
1045: stl += lc[i];
1046: }
1047: ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);
1049: /* Create MATIS */
1050: MatCreate(comm,&B);
1051: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1052: MatGetBlockSizes(A,&rbs,&cbs);
1053: MatSetBlockSizes(B,rbs,cbs);
1054: MatSetType(B,MATIS);
1055: MatISSetLocalMatType(B,MATNEST);
1056: { /* hack : avoid setup of scatters */
1057: Mat_IS *matis = (Mat_IS*)(B->data);
1058: matis->islocalref = PETSC_TRUE;
1059: }
1060: MatSetLocalToGlobalMapping(B,rl2g,cl2g);
1061: ISLocalToGlobalMappingDestroy(&rl2g);
1062: ISLocalToGlobalMappingDestroy(&cl2g);
1063: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1064: MatNestSetVecType(lA,VECNEST);
1065: for (i=0;i<nr*nc;i++) {
1066: if (istrans[i]) {
1067: MatDestroy(&snest[i]);
1068: }
1069: }
1070: MatISSetLocalMat(B,lA);
1071: MatDestroy(&lA);
1072: { /* hack : setup of scatters done here */
1073: Mat_IS *matis = (Mat_IS*)(B->data);
1075: matis->islocalref = PETSC_FALSE;
1076: MatISSetUpScatters_Private(B);
1077: }
1078: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1079: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1080: if (reuse == MAT_INPLACE_MATRIX) {
1081: MatHeaderReplace(A,&B);
1082: } else {
1083: *newmat = B;
1084: }
1085: } else {
1086: if (lreuse) {
1087: MatISGetLocalMat(*newmat,&lA);
1088: for (i=0;i<nr;i++) {
1089: for (j=0;j<nc;j++) {
1090: if (snest[i*nc+j]) {
1091: MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
1092: if (istrans[i*nc+j]) {
1093: MatDestroy(&snest[i*nc+j]);
1094: }
1095: }
1096: }
1097: }
1098: } else {
1099: PetscInt stl;
1100: for (i=0,stl=0;i<nr;i++) {
1101: ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
1102: stl += lr[i];
1103: }
1104: for (i=0,stl=0;i<nc;i++) {
1105: ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
1106: stl += lc[i];
1107: }
1108: MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
1109: for (i=0;i<nr*nc;i++) {
1110: if (istrans[i]) {
1111: MatDestroy(&snest[i]);
1112: }
1113: }
1114: MatISSetLocalMat(*newmat,lA);
1115: MatDestroy(&lA);
1116: }
1117: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1118: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1119: }
1121: /* Create local matrix in MATNEST format */
1122: convert = PETSC_FALSE;
1123: PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
1124: if (convert) {
1125: Mat M;
1126: MatISLocalFields lf;
1127: PetscContainer c;
1129: MatISGetLocalMat(*newmat,&lA);
1130: MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
1131: MatISSetLocalMat(*newmat,M);
1132: MatDestroy(&M);
1134: /* attach local fields to the matrix */
1135: PetscNew(&lf);
1136: PetscCalloc2(nr,&lf->rf,nc,&lf->cf);
1137: for (i=0;i<nr;i++) {
1138: PetscInt n,st;
1140: ISGetLocalSize(islrow[i],&n);
1141: ISStrideGetInfo(islrow[i],&st,NULL);
1142: ISCreateStride(comm,n,st,1,&lf->rf[i]);
1143: }
1144: for (i=0;i<nc;i++) {
1145: PetscInt n,st;
1147: ISGetLocalSize(islcol[i],&n);
1148: ISStrideGetInfo(islcol[i],&st,NULL);
1149: ISCreateStride(comm,n,st,1,&lf->cf[i]);
1150: }
1151: lf->nr = nr;
1152: lf->nc = nc;
1153: PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
1154: PetscContainerSetPointer(c,lf);
1155: PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
1156: PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
1157: PetscContainerDestroy(&c);
1158: }
1160: /* Free workspace */
1161: for (i=0;i<nr;i++) {
1162: ISDestroy(&islrow[i]);
1163: }
1164: for (i=0;i<nc;i++) {
1165: ISDestroy(&islcol[i]);
1166: }
1167: PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
1168: PetscFree2(lr,lc);
1169: return(0);
1170: }
1172: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1173: {
1174: Mat_IS *matis = (Mat_IS*)A->data;
1175: Vec ll,rr;
1176: const PetscScalar *Y,*X;
1177: PetscScalar *x,*y;
1178: PetscErrorCode ierr;
1181: MatISSetUpSF(A);
1182: if (l) {
1183: ll = matis->y;
1184: VecGetArrayRead(l,&Y);
1185: VecGetArray(ll,&y);
1186: PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y);
1187: } else {
1188: ll = NULL;
1189: }
1190: if (r) {
1191: rr = matis->x;
1192: VecGetArrayRead(r,&X);
1193: VecGetArray(rr,&x);
1194: PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x);
1195: } else {
1196: rr = NULL;
1197: }
1198: if (ll) {
1199: PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y);
1200: VecRestoreArrayRead(l,&Y);
1201: VecRestoreArray(ll,&y);
1202: }
1203: if (rr) {
1204: PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x);
1205: VecRestoreArrayRead(r,&X);
1206: VecRestoreArray(rr,&x);
1207: }
1208: MatDiagonalScale(matis->A,ll,rr);
1209: return(0);
1210: }
1212: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
1213: {
1214: Mat_IS *matis = (Mat_IS*)A->data;
1215: MatInfo info;
1216: PetscReal isend[6],irecv[6];
1217: PetscInt bs;
1221: MatGetBlockSize(A,&bs);
1222: if (matis->A->ops->getinfo) {
1223: MatGetInfo(matis->A,MAT_LOCAL,&info);
1224: isend[0] = info.nz_used;
1225: isend[1] = info.nz_allocated;
1226: isend[2] = info.nz_unneeded;
1227: isend[3] = info.memory;
1228: isend[4] = info.mallocs;
1229: } else {
1230: isend[0] = 0.;
1231: isend[1] = 0.;
1232: isend[2] = 0.;
1233: isend[3] = 0.;
1234: isend[4] = 0.;
1235: }
1236: isend[5] = matis->A->num_ass;
1237: if (flag == MAT_LOCAL) {
1238: ginfo->nz_used = isend[0];
1239: ginfo->nz_allocated = isend[1];
1240: ginfo->nz_unneeded = isend[2];
1241: ginfo->memory = isend[3];
1242: ginfo->mallocs = isend[4];
1243: ginfo->assemblies = isend[5];
1244: } else if (flag == MAT_GLOBAL_MAX) {
1245: MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));
1247: ginfo->nz_used = irecv[0];
1248: ginfo->nz_allocated = irecv[1];
1249: ginfo->nz_unneeded = irecv[2];
1250: ginfo->memory = irecv[3];
1251: ginfo->mallocs = irecv[4];
1252: ginfo->assemblies = irecv[5];
1253: } else if (flag == MAT_GLOBAL_SUM) {
1254: MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));
1256: ginfo->nz_used = irecv[0];
1257: ginfo->nz_allocated = irecv[1];
1258: ginfo->nz_unneeded = irecv[2];
1259: ginfo->memory = irecv[3];
1260: ginfo->mallocs = irecv[4];
1261: ginfo->assemblies = A->num_ass;
1262: }
1263: ginfo->block_size = bs;
1264: ginfo->fill_ratio_given = 0;
1265: ginfo->fill_ratio_needed = 0;
1266: ginfo->factor_mallocs = 0;
1267: return(0);
1268: }
1270: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
1271: {
1272: Mat C,lC,lA;
1273: PetscErrorCode ierr;
1276: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1277: ISLocalToGlobalMapping rl2g,cl2g;
1278: MatCreate(PetscObjectComm((PetscObject)A),&C);
1279: MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
1280: MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
1281: MatSetType(C,MATIS);
1282: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
1283: MatSetLocalToGlobalMapping(C,cl2g,rl2g);
1284: } else {
1285: C = *B;
1286: }
1288: /* perform local transposition */
1289: MatISGetLocalMat(A,&lA);
1290: MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
1291: MatISSetLocalMat(C,lC);
1292: MatDestroy(&lC);
1294: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1295: *B = C;
1296: } else {
1297: MatHeaderMerge(A,&C);
1298: }
1299: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1300: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1301: return(0);
1302: }
1304: PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
1305: {
1306: Mat_IS *is = (Mat_IS*)A->data;
1310: if (D) { /* MatShift_IS pass D = NULL */
1311: VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1312: VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1313: }
1314: VecPointwiseDivide(is->y,is->y,is->counter);
1315: MatDiagonalSet(is->A,is->y,insmode);
1316: return(0);
1317: }
1319: PetscErrorCode MatShift_IS(Mat A,PetscScalar a)
1320: {
1321: Mat_IS *is = (Mat_IS*)A->data;
1325: VecSet(is->y,a);
1326: MatDiagonalSet_IS(A,NULL,ADD_VALUES);
1327: return(0);
1328: }
1330: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1331: {
1333: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1336: #if defined(PETSC_USE_DEBUG)
1337: 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);
1338: #endif
1339: ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
1340: ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
1341: MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1342: return(0);
1343: }
1345: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
1346: {
1348: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
1351: #if defined(PETSC_USE_DEBUG)
1352: 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);
1353: #endif
1354: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
1355: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
1356: MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1357: return(0);
1358: }
1360: static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
1361: {
1362: PetscInt *owners = map->range;
1363: PetscInt n = map->n;
1364: PetscSF sf;
1365: PetscInt *lidxs,*work = NULL;
1366: PetscSFNode *ridxs;
1367: PetscMPIInt rank;
1368: PetscInt r, p = 0, len = 0;
1372: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
1373: /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */
1374: MPI_Comm_rank(map->comm,&rank);
1375: PetscMalloc1(n,&lidxs);
1376: for (r = 0; r < n; ++r) lidxs[r] = -1;
1377: PetscMalloc1(N,&ridxs);
1378: for (r = 0; r < N; ++r) {
1379: const PetscInt idx = idxs[r];
1380: if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
1381: if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
1382: PetscLayoutFindOwner(map,idx,&p);
1383: }
1384: ridxs[r].rank = p;
1385: ridxs[r].index = idxs[r] - owners[p];
1386: }
1387: PetscSFCreate(map->comm,&sf);
1388: PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
1389: PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
1390: PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
1391: if (ogidxs) { /* communicate global idxs */
1392: PetscInt cum = 0,start,*work2;
1394: PetscMalloc1(n,&work);
1395: PetscCalloc1(N,&work2);
1396: for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
1397: MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
1398: start -= cum;
1399: cum = 0;
1400: for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
1401: PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);
1402: PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);
1403: PetscFree(work2);
1404: }
1405: PetscSFDestroy(&sf);
1406: /* Compress and put in indices */
1407: for (r = 0; r < n; ++r)
1408: if (lidxs[r] >= 0) {
1409: if (work) work[len] = work[r];
1410: lidxs[len++] = r;
1411: }
1412: if (on) *on = len;
1413: if (oidxs) *oidxs = lidxs;
1414: if (ogidxs) *ogidxs = work;
1415: return(0);
1416: }
1418: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
1419: {
1420: Mat locmat,newlocmat;
1421: Mat_IS *newmatis;
1422: #if defined(PETSC_USE_DEBUG)
1423: Vec rtest,ltest;
1424: const PetscScalar *array;
1425: #endif
1426: const PetscInt *idxs;
1427: PetscInt i,m,n;
1428: PetscErrorCode ierr;
1431: if (scall == MAT_REUSE_MATRIX) {
1432: PetscBool ismatis;
1434: PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
1435: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
1436: newmatis = (Mat_IS*)(*newmat)->data;
1437: if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
1438: if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
1439: }
1440: /* irow and icol may not have duplicate entries */
1441: #if defined(PETSC_USE_DEBUG)
1442: MatCreateVecs(mat,<est,&rtest);
1443: ISGetLocalSize(irow,&n);
1444: ISGetIndices(irow,&idxs);
1445: for (i=0;i<n;i++) {
1446: VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
1447: }
1448: VecAssemblyBegin(rtest);
1449: VecAssemblyEnd(rtest);
1450: VecGetLocalSize(rtest,&n);
1451: VecGetOwnershipRange(rtest,&m,NULL);
1452: VecGetArrayRead(rtest,&array);
1453: 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]));
1454: VecRestoreArrayRead(rtest,&array);
1455: ISRestoreIndices(irow,&idxs);
1456: ISGetLocalSize(icol,&n);
1457: ISGetIndices(icol,&idxs);
1458: for (i=0;i<n;i++) {
1459: VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
1460: }
1461: VecAssemblyBegin(ltest);
1462: VecAssemblyEnd(ltest);
1463: VecGetLocalSize(ltest,&n);
1464: VecGetOwnershipRange(ltest,&m,NULL);
1465: VecGetArrayRead(ltest,&array);
1466: 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]));
1467: VecRestoreArrayRead(ltest,&array);
1468: ISRestoreIndices(icol,&idxs);
1469: VecDestroy(&rtest);
1470: VecDestroy(<est);
1471: #endif
1472: if (scall == MAT_INITIAL_MATRIX) {
1473: Mat_IS *matis = (Mat_IS*)mat->data;
1474: ISLocalToGlobalMapping rl2g;
1475: IS is;
1476: PetscInt *lidxs,*lgidxs,*newgidxs;
1477: PetscInt ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
1478: PetscBool cong;
1479: MPI_Comm comm;
1481: PetscObjectGetComm((PetscObject)mat,&comm);
1482: MatGetBlockSizes(mat,&arbs,&acbs);
1483: ISGetBlockSize(irow,&irbs);
1484: ISGetBlockSize(icol,&icbs);
1485: rbs = arbs == irbs ? irbs : 1;
1486: cbs = acbs == icbs ? icbs : 1;
1487: ISGetLocalSize(irow,&m);
1488: ISGetLocalSize(icol,&n);
1489: MatCreate(comm,newmat);
1490: MatSetType(*newmat,MATIS);
1491: MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
1492: MatSetBlockSizes(*newmat,rbs,cbs);
1493: /* communicate irow to their owners in the layout */
1494: ISGetIndices(irow,&idxs);
1495: PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
1496: ISRestoreIndices(irow,&idxs);
1497: MatISSetUpSF(mat);
1498: PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));
1499: for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
1500: PetscFree(lidxs);
1501: PetscFree(lgidxs);
1502: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1503: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1504: for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
1505: PetscMalloc1(newloc,&newgidxs);
1506: PetscMalloc1(newloc,&lidxs);
1507: for (i=0,newloc=0;i<matis->sf->nleaves;i++)
1508: if (matis->sf_leafdata[i]) {
1509: lidxs[newloc] = i;
1510: newgidxs[newloc++] = matis->sf_leafdata[i]-1;
1511: }
1512: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1513: ISLocalToGlobalMappingCreateIS(is,&rl2g);
1514: ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);
1515: ISDestroy(&is);
1516: /* local is to extract local submatrix */
1517: newmatis = (Mat_IS*)(*newmat)->data;
1518: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
1519: MatHasCongruentLayouts(mat,&cong);
1520: if (cong && irow == icol && matis->csf == matis->sf) {
1521: MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
1522: PetscObjectReference((PetscObject)newmatis->getsub_ris);
1523: newmatis->getsub_cis = newmatis->getsub_ris;
1524: } else {
1525: ISLocalToGlobalMapping cl2g;
1527: /* communicate icol to their owners in the layout */
1528: ISGetIndices(icol,&idxs);
1529: PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
1530: ISRestoreIndices(icol,&idxs);
1531: PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));
1532: for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
1533: PetscFree(lidxs);
1534: PetscFree(lgidxs);
1535: PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
1536: PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
1537: for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
1538: PetscMalloc1(newloc,&newgidxs);
1539: PetscMalloc1(newloc,&lidxs);
1540: for (i=0,newloc=0;i<matis->csf->nleaves;i++)
1541: if (matis->csf_leafdata[i]) {
1542: lidxs[newloc] = i;
1543: newgidxs[newloc++] = matis->csf_leafdata[i]-1;
1544: }
1545: ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1546: ISLocalToGlobalMappingCreateIS(is,&cl2g);
1547: ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);
1548: ISDestroy(&is);
1549: /* local is to extract local submatrix */
1550: ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
1551: MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
1552: ISLocalToGlobalMappingDestroy(&cl2g);
1553: }
1554: ISLocalToGlobalMappingDestroy(&rl2g);
1555: } else {
1556: MatISGetLocalMat(*newmat,&newlocmat);
1557: }
1558: MatISGetLocalMat(mat,&locmat);
1559: newmatis = (Mat_IS*)(*newmat)->data;
1560: MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
1561: if (scall == MAT_INITIAL_MATRIX) {
1562: MatISSetLocalMat(*newmat,newlocmat);
1563: MatDestroy(&newlocmat);
1564: }
1565: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1566: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1567: return(0);
1568: }
1570: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
1571: {
1572: Mat_IS *a = (Mat_IS*)A->data,*b;
1573: PetscBool ismatis;
1577: PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
1578: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
1579: b = (Mat_IS*)B->data;
1580: MatCopy(a->A,b->A,str);
1581: PetscObjectStateIncrease((PetscObject)B);
1582: return(0);
1583: }
1585: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d)
1586: {
1587: Vec v;
1588: const PetscScalar *array;
1589: PetscInt i,n;
1590: PetscErrorCode ierr;
1593: *missing = PETSC_FALSE;
1594: MatCreateVecs(A,NULL,&v);
1595: MatGetDiagonal(A,v);
1596: VecGetLocalSize(v,&n);
1597: VecGetArrayRead(v,&array);
1598: for (i=0;i<n;i++) if (array[i] == 0.) break;
1599: VecRestoreArrayRead(v,&array);
1600: VecDestroy(&v);
1601: if (i != n) *missing = PETSC_TRUE;
1602: if (d) {
1603: *d = -1;
1604: if (*missing) {
1605: PetscInt rstart;
1606: MatGetOwnershipRange(A,&rstart,NULL);
1607: *d = i+rstart;
1608: }
1609: }
1610: return(0);
1611: }
1613: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1614: {
1615: Mat_IS *matis = (Mat_IS*)(B->data);
1616: const PetscInt *gidxs;
1617: PetscInt nleaves;
1621: if (matis->sf) return(0);
1622: PetscLayoutSetUp(B->rmap);
1623: PetscLayoutSetUp(B->cmap);
1624: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
1625: ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
1626: ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
1627: PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1628: ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
1629: PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
1630: if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
1631: ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
1632: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
1633: ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
1634: PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1635: ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
1636: PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
1637: } else {
1638: matis->csf = matis->sf;
1639: matis->csf_leafdata = matis->sf_leafdata;
1640: matis->csf_rootdata = matis->sf_rootdata;
1641: }
1642: return(0);
1643: }
1645: /*@
1646: MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP.
1648: Collective on MPI_Comm
1650: Input Parameters:
1651: + A - the matrix
1652: - store - the boolean flag
1654: Level: advanced
1656: Notes:
1658: .keywords: matrix
1660: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
1661: @*/
1662: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1663: {
1670: PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
1671: return(0);
1672: }
1674: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1675: {
1676: Mat_IS *matis = (Mat_IS*)(A->data);
1680: matis->storel2l = store;
1681: if (!store) {
1682: PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
1683: }
1684: return(0);
1685: }
1687: /*@
1688: MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1690: Collective on MPI_Comm
1692: Input Parameters:
1693: + A - the matrix
1694: - fix - the boolean flag
1696: Level: advanced
1698: Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process.
1700: .keywords: matrix
1702: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
1703: @*/
1704: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1705: {
1712: PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
1713: return(0);
1714: }
1716: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1717: {
1718: Mat_IS *matis = (Mat_IS*)(A->data);
1721: matis->locempty = fix;
1722: return(0);
1723: }
1725: /*@
1726: MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.
1728: Collective on MPI_Comm
1730: Input Parameters:
1731: + B - the matrix
1732: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1733: (same value is used for all local rows)
1734: . d_nnz - array containing the number of nonzeros in the various rows of the
1735: DIAGONAL portion of the local submatrix (possibly different for each row)
1736: or NULL, if d_nz is used to specify the nonzero structure.
1737: The size of this array is equal to the number of local rows, i.e 'm'.
1738: For matrices that will be factored, you must leave room for (and set)
1739: the diagonal entry even if it is zero.
1740: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1741: submatrix (same value is used for all local rows).
1742: - o_nnz - array containing the number of nonzeros in the various rows of the
1743: OFF-DIAGONAL portion of the local submatrix (possibly different for
1744: each row) or NULL, if o_nz is used to specify the nonzero
1745: structure. The size of this array is equal to the number
1746: of local rows, i.e 'm'.
1748: If the *_nnz parameter is given then the *_nz parameter is ignored
1750: Level: intermediate
1752: Notes:
1753: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
1754: from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
1755: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1757: .keywords: matrix
1759: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1760: @*/
1761: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1762: {
1768: PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1769: return(0);
1770: }
1772: /* this is used by DMDA */
1773: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1774: {
1775: Mat_IS *matis = (Mat_IS*)(B->data);
1776: PetscInt bs,i,nlocalcols;
1780: if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
1781: MatISSetUpSF(B);
1783: if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
1784: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
1786: if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz;
1787: else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
1789: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1790: MatGetSize(matis->A,NULL,&nlocalcols);
1791: MatGetBlockSize(matis->A,&bs);
1792: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1794: for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1795: MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1796: #if defined(PETSC_HAVE_HYPRE)
1797: MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1798: #endif
1800: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
1801: MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1803: nlocalcols /= bs;
1804: for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols - i);
1805: MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
1807: /* for other matrix types */
1808: MatSetUp(matis->A);
1809: return(0);
1810: }
1812: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1813: {
1814: Mat_IS *matis = (Mat_IS*)(A->data);
1815: PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1816: const PetscInt *global_indices_r,*global_indices_c;
1817: PetscInt i,j,bs,rows,cols;
1818: PetscInt lrows,lcols;
1819: PetscInt local_rows,local_cols;
1820: PetscMPIInt size;
1821: PetscBool isdense,issbaij;
1822: PetscErrorCode ierr;
1825: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1826: MatGetSize(A,&rows,&cols);
1827: MatGetBlockSize(A,&bs);
1828: MatGetSize(matis->A,&local_rows,&local_cols);
1829: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1830: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1831: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
1832: if (A->rmap->mapping != A->cmap->mapping) {
1833: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
1834: } else {
1835: global_indices_c = global_indices_r;
1836: }
1838: if (issbaij) {
1839: MatGetRowUpperTriangular(matis->A);
1840: }
1841: /*
1842: An SF reduce is needed to sum up properly on shared rows.
1843: Note that generally preallocation is not exact, since it overestimates nonzeros
1844: */
1845: MatISSetUpSF(A);
1846: MatGetLocalSize(A,&lrows,&lcols);
1847: MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1848: /* All processes need to compute entire row ownership */
1849: PetscMalloc1(rows,&row_ownership);
1850: MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1851: for (i=0;i<size;i++) {
1852: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
1853: row_ownership[j] = i;
1854: }
1855: }
1856: MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);
1858: /*
1859: my_dnz and my_onz contains exact contribution to preallocation from each local mat
1860: then, they will be summed up properly. This way, preallocation is always sufficient
1861: */
1862: PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1863: /* preallocation as a MATAIJ */
1864: if (isdense) { /* special case for dense local matrices */
1865: for (i=0;i<local_rows;i++) {
1866: PetscInt owner = row_ownership[global_indices_r[i]];
1867: for (j=0;j<local_cols;j++) {
1868: PetscInt index_col = global_indices_c[j];
1869: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1870: my_dnz[i] += 1;
1871: } else { /* offdiag block */
1872: my_onz[i] += 1;
1873: }
1874: }
1875: }
1876: } else if (matis->A->ops->getrowij) {
1877: const PetscInt *ii,*jj,*jptr;
1878: PetscBool done;
1879: MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1880: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1881: jptr = jj;
1882: for (i=0;i<local_rows;i++) {
1883: PetscInt index_row = global_indices_r[i];
1884: for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1885: PetscInt owner = row_ownership[index_row];
1886: PetscInt index_col = global_indices_c[*jptr];
1887: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1888: my_dnz[i] += 1;
1889: } else { /* offdiag block */
1890: my_onz[i] += 1;
1891: }
1892: /* same as before, interchanging rows and cols */
1893: if (issbaij && index_col != index_row) {
1894: owner = row_ownership[index_col];
1895: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1896: my_dnz[*jptr] += 1;
1897: } else {
1898: my_onz[*jptr] += 1;
1899: }
1900: }
1901: }
1902: }
1903: MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1904: if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1905: } else { /* loop over rows and use MatGetRow */
1906: for (i=0;i<local_rows;i++) {
1907: const PetscInt *cols;
1908: PetscInt ncols,index_row = global_indices_r[i];
1909: MatGetRow(matis->A,i,&ncols,&cols,NULL);
1910: for (j=0;j<ncols;j++) {
1911: PetscInt owner = row_ownership[index_row];
1912: PetscInt index_col = global_indices_c[cols[j]];
1913: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1914: my_dnz[i] += 1;
1915: } else { /* offdiag block */
1916: my_onz[i] += 1;
1917: }
1918: /* same as before, interchanging rows and cols */
1919: if (issbaij && index_col != index_row) {
1920: owner = row_ownership[index_col];
1921: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1922: my_dnz[cols[j]] += 1;
1923: } else {
1924: my_onz[cols[j]] += 1;
1925: }
1926: }
1927: }
1928: MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1929: }
1930: }
1931: if (global_indices_c != global_indices_r) {
1932: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1933: }
1934: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1935: PetscFree(row_ownership);
1937: /* Reduce my_dnz and my_onz */
1938: if (maxreduce) {
1939: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1940: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1941: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1942: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1943: } else {
1944: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1945: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1946: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1947: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1948: }
1949: PetscFree2(my_dnz,my_onz);
1951: /* Resize preallocation if overestimated */
1952: for (i=0;i<lrows;i++) {
1953: dnz[i] = PetscMin(dnz[i],lcols);
1954: onz[i] = PetscMin(onz[i],cols-lcols);
1955: }
1957: /* Set preallocation */
1958: MatSeqAIJSetPreallocation(B,0,dnz);
1959: MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1960: for (i=0;i<lrows;i+=bs) {
1961: PetscInt b, d = dnz[i],o = onz[i];
1963: for (b=1;b<bs;b++) {
1964: d = PetscMax(d,dnz[i+b]);
1965: o = PetscMax(o,onz[i+b]);
1966: }
1967: dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
1968: onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
1969: }
1970: MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1971: MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1972: MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1973: MatPreallocateFinalize(dnz,onz);
1974: if (issbaij) {
1975: MatRestoreRowUpperTriangular(matis->A);
1976: }
1977: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1978: return(0);
1979: }
1981: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1982: {
1983: Mat_IS *matis = (Mat_IS*)(mat->data);
1984: Mat local_mat,MT;
1985: PetscInt rbs,cbs,rows,cols,lrows,lcols;
1986: PetscInt local_rows,local_cols;
1987: PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij;
1988: #if defined (PETSC_USE_DEBUG)
1989: PetscBool lb[4],bb[4];
1990: #endif
1991: PetscMPIInt size;
1992: PetscScalar *array;
1996: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
1997: if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1998: Mat B;
1999: IS irows = NULL,icols = NULL;
2000: PetscInt rbs,cbs;
2002: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2003: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2004: if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
2005: IS rows,cols;
2006: const PetscInt *ridxs,*cidxs;
2007: PetscInt i,nw,*work;
2009: ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);
2010: ISLocalToGlobalMappingGetSize(mat->rmap->mapping,&nw);
2011: nw = nw/rbs;
2012: PetscCalloc1(nw,&work);
2013: for (i=0;i<nw;i++) work[ridxs[i]] += 1;
2014: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
2015: if (i == nw) {
2016: ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);
2017: ISSetPermutation(rows);
2018: ISInvertPermutation(rows,PETSC_DECIDE,&irows);
2019: ISDestroy(&rows);
2020: }
2021: ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);
2022: PetscFree(work);
2023: if (irows && mat->rmap->mapping != mat->cmap->mapping) {
2024: ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);
2025: ISLocalToGlobalMappingGetSize(mat->cmap->mapping,&nw);
2026: nw = nw/cbs;
2027: PetscCalloc1(nw,&work);
2028: for (i=0;i<nw;i++) work[cidxs[i]] += 1;
2029: for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
2030: if (i == nw) {
2031: ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);
2032: ISSetPermutation(cols);
2033: ISInvertPermutation(cols,PETSC_DECIDE,&icols);
2034: ISDestroy(&cols);
2035: }
2036: ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);
2037: PetscFree(work);
2038: } else if (irows) {
2039: PetscObjectReference((PetscObject)irows);
2040: icols = irows;
2041: }
2042: } else {
2043: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);
2044: PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);
2045: if (irows) { PetscObjectReference((PetscObject)irows); }
2046: if (icols) { PetscObjectReference((PetscObject)icols); }
2047: }
2048: if (!irows || !icols) {
2049: ISDestroy(&icols);
2050: ISDestroy(&irows);
2051: goto general_assembly;
2052: }
2053: MatConvert(matis->A,(rbs == cbs && rbs > 1) ? MATSEQBAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
2054: if (reuse != MAT_INPLACE_MATRIX) {
2055: MatCreateSubMatrix(B,irows,icols,reuse,M);
2056: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);
2057: PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);
2058: } else {
2059: Mat C;
2061: MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
2062: MatHeaderReplace(mat,&C);
2063: }
2064: MatDestroy(&B);
2065: ISDestroy(&icols);
2066: ISDestroy(&irows);
2067: return(0);
2068: }
2069: general_assembly:
2070: MatGetSize(mat,&rows,&cols);
2071: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2072: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2073: MatGetLocalSize(mat,&lrows,&lcols);
2074: MatGetSize(matis->A,&local_rows,&local_cols);
2075: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
2076: PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
2077: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
2078: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
2079: if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
2080: #if defined (PETSC_USE_DEBUG)
2081: lb[0] = isseqdense;
2082: lb[1] = isseqaij;
2083: lb[2] = isseqbaij;
2084: lb[3] = isseqsbaij;
2085: MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
2086: if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
2087: #endif
2089: if (reuse != MAT_REUSE_MATRIX) {
2090: MatCreate(PetscObjectComm((PetscObject)mat),&MT);
2091: MatSetSizes(MT,lrows,lcols,rows,cols);
2092: MatSetType(MT,mtype);
2093: MatSetBlockSizes(MT,rbs,cbs);
2094: MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
2095: } else {
2096: PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;
2098: /* some checks */
2099: MT = *M;
2100: MatGetBlockSizes(MT,&mrbs,&mcbs);
2101: MatGetSize(MT,&mrows,&mcols);
2102: MatGetLocalSize(MT,&mlrows,&mlcols);
2103: if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
2104: if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
2105: if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
2106: if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
2107: if (mrbs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%d != %d)",rbs,mrbs);
2108: if (mcbs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%d != %d)",cbs,mcbs);
2109: MatZeroEntries(MT);
2110: }
2112: if (isseqsbaij || isseqbaij) {
2113: MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
2114: isseqaij = PETSC_TRUE;
2115: } else {
2116: PetscObjectReference((PetscObject)matis->A);
2117: local_mat = matis->A;
2118: }
2120: /* Set values */
2121: MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);
2122: if (isseqdense) { /* special case for dense local matrices */
2123: PetscInt i,*dummy;
2125: PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
2126: for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
2127: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
2128: MatDenseGetArray(local_mat,&array);
2129: MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
2130: MatDenseRestoreArray(local_mat,&array);
2131: PetscFree(dummy);
2132: } else if (isseqaij) {
2133: const PetscInt *blocks;
2134: PetscInt i,nvtxs,*xadj,*adjncy, nb;
2135: PetscBool done;
2137: MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2138: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
2139: MatSeqAIJGetArray(local_mat,&array);
2140: MatGetVariableBlockSizes(local_mat,&nb,&blocks);
2141: if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2142: PetscInt sum;
2144: for (i=0,sum=0;i<nb;i++) sum += blocks[i];
2145: if (sum == nvtxs) {
2146: PetscInt r;
2148: for (i=0,r=0;i<nb;i++) {
2149: 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]);
2150: MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],array+xadj[r],ADD_VALUES);
2151: r += blocks[i];
2152: }
2153: } else {
2154: for (i=0;i<nvtxs;i++) {
2155: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
2156: }
2157: }
2158: } else {
2159: for (i=0;i<nvtxs;i++) {
2160: MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
2161: }
2162: }
2163: MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2164: if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
2165: MatSeqAIJRestoreArray(local_mat,&array);
2166: } else { /* very basic values insertion for all other matrix types */
2167: PetscInt i;
2169: for (i=0;i<local_rows;i++) {
2170: PetscInt j;
2171: const PetscInt *local_indices_cols;
2173: MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
2174: MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
2175: MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
2176: }
2177: }
2178: MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
2179: MatDestroy(&local_mat);
2180: MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
2181: if (isseqdense) {
2182: MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
2183: }
2184: if (reuse == MAT_INPLACE_MATRIX) {
2185: MatHeaderReplace(mat,&MT);
2186: } else if (reuse == MAT_INITIAL_MATRIX) {
2187: *M = MT;
2188: }
2189: return(0);
2190: }
2192: /*@
2193: MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format
2195: Input Parameter:
2196: . mat - the matrix (should be of type MATIS)
2197: . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
2199: Output Parameter:
2200: . newmat - the matrix in AIJ format
2202: Level: developer
2204: Notes:
2205: This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface.
2207: .seealso: MATIS, MatConvert()
2208: @*/
2209: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2210: {
2217: if (reuse == MAT_REUSE_MATRIX) {
2220: if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
2221: }
2222: PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
2223: return(0);
2224: }
2226: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
2227: {
2229: Mat_IS *matis = (Mat_IS*)(mat->data);
2230: PetscInt rbs,cbs,m,n,M,N;
2231: Mat B,localmat;
2234: ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2235: ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2236: MatGetSize(mat,&M,&N);
2237: MatGetLocalSize(mat,&m,&n);
2238: MatCreate(PetscObjectComm((PetscObject)mat),&B);
2239: MatSetSizes(B,m,n,M,N);
2240: MatSetBlockSize(B,rbs == cbs ? rbs : 1);
2241: MatSetType(B,MATIS);
2242: MatISSetLocalMatType(B,matis->lmattype);
2243: MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
2244: MatDuplicate(matis->A,op,&localmat);
2245: MatISSetLocalMat(B,localmat);
2246: MatDestroy(&localmat);
2247: if (matis->sf) {
2248: Mat_IS *bmatis = (Mat_IS*)(B->data);
2250: PetscObjectReference((PetscObject)matis->sf);
2251: bmatis->sf = matis->sf;
2252: PetscMalloc2(matis->sf->nroots,&bmatis->sf_rootdata,matis->sf->nleaves,&bmatis->sf_leafdata);
2253: if (matis->sf != matis->csf) {
2254: PetscObjectReference((PetscObject)matis->csf);
2255: bmatis->csf = matis->csf;
2256: PetscMalloc2(matis->csf->nroots,&bmatis->csf_rootdata,matis->csf->nleaves,&bmatis->csf_leafdata);
2257: } else {
2258: bmatis->csf = bmatis->sf;
2259: bmatis->csf_leafdata = bmatis->sf_leafdata;
2260: bmatis->csf_rootdata = bmatis->sf_rootdata;
2261: }
2262: }
2263: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2264: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2265: *newmat = B;
2266: return(0);
2267: }
2269: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
2270: {
2272: Mat_IS *matis = (Mat_IS*)A->data;
2273: PetscBool local_sym;
2276: MatIsHermitian(matis->A,tol,&local_sym);
2277: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2278: return(0);
2279: }
2281: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
2282: {
2284: Mat_IS *matis = (Mat_IS*)A->data;
2285: PetscBool local_sym;
2288: MatIsSymmetric(matis->A,tol,&local_sym);
2289: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2290: return(0);
2291: }
2293: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg)
2294: {
2296: Mat_IS *matis = (Mat_IS*)A->data;
2297: PetscBool local_sym;
2300: if (A->rmap->mapping != A->cmap->mapping) {
2301: *flg = PETSC_FALSE;
2302: return(0);
2303: }
2304: MatIsStructurallySymmetric(matis->A,&local_sym);
2305: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2306: return(0);
2307: }
2309: static PetscErrorCode MatDestroy_IS(Mat A)
2310: {
2312: Mat_IS *b = (Mat_IS*)A->data;
2315: PetscFree(b->lmattype);
2316: MatDestroy(&b->A);
2317: VecScatterDestroy(&b->cctx);
2318: VecScatterDestroy(&b->rctx);
2319: VecDestroy(&b->x);
2320: VecDestroy(&b->y);
2321: VecDestroy(&b->counter);
2322: ISDestroy(&b->getsub_ris);
2323: ISDestroy(&b->getsub_cis);
2324: if (b->sf != b->csf) {
2325: PetscSFDestroy(&b->csf);
2326: PetscFree2(b->csf_rootdata,b->csf_leafdata);
2327: } else b->csf = NULL;
2328: PetscSFDestroy(&b->sf);
2329: PetscFree2(b->sf_rootdata,b->sf_leafdata);
2330: PetscFree(A->data);
2331: PetscObjectChangeTypeName((PetscObject)A,0);
2332: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
2333: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
2334: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
2335: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
2336: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
2337: PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);
2338: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
2339: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
2340: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
2341: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
2342: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
2343: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
2344: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
2345: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
2346: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
2347: return(0);
2348: }
2350: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
2351: {
2353: Mat_IS *is = (Mat_IS*)A->data;
2354: PetscScalar zero = 0.0;
2357: /* scatter the global vector x into the local work vector */
2358: VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2359: VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2361: /* multiply the local matrix */
2362: MatMult(is->A,is->x,is->y);
2364: /* scatter product back into global memory */
2365: VecSet(y,zero);
2366: VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2367: VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2368: return(0);
2369: }
2371: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2372: {
2373: Vec temp_vec;
2377: if (v3 != v2) {
2378: MatMult(A,v1,v3);
2379: VecAXPY(v3,1.0,v2);
2380: } else {
2381: VecDuplicate(v2,&temp_vec);
2382: MatMult(A,v1,temp_vec);
2383: VecAXPY(temp_vec,1.0,v2);
2384: VecCopy(temp_vec,v3);
2385: VecDestroy(&temp_vec);
2386: }
2387: return(0);
2388: }
2390: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
2391: {
2392: Mat_IS *is = (Mat_IS*)A->data;
2396: /* scatter the global vector x into the local work vector */
2397: VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2398: VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
2400: /* multiply the local matrix */
2401: MatMultTranspose(is->A,is->y,is->x);
2403: /* scatter product back into global vector */
2404: VecSet(x,0);
2405: VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2406: VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2407: return(0);
2408: }
2410: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2411: {
2412: Vec temp_vec;
2416: if (v3 != v2) {
2417: MatMultTranspose(A,v1,v3);
2418: VecAXPY(v3,1.0,v2);
2419: } else {
2420: VecDuplicate(v2,&temp_vec);
2421: MatMultTranspose(A,v1,temp_vec);
2422: VecAXPY(temp_vec,1.0,v2);
2423: VecCopy(temp_vec,v3);
2424: VecDestroy(&temp_vec);
2425: }
2426: return(0);
2427: }
2429: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
2430: {
2431: Mat_IS *a = (Mat_IS*)A->data;
2433: PetscViewer sviewer;
2434: PetscBool isascii,view = PETSC_TRUE;
2437: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2438: if (isascii) {
2439: PetscViewerFormat format;
2441: PetscViewerGetFormat(viewer,&format);
2442: if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2443: }
2444: if (!view) return(0);
2445: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2446: MatView(a->A,sviewer);
2447: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2448: PetscViewerFlush(viewer);
2449: return(0);
2450: }
2452: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2453: {
2454: Vec cglobal,rglobal;
2455: IS from;
2456: Mat_IS *is = (Mat_IS*)A->data;
2457: const PetscInt *garray;
2458: PetscInt nr,rbs,nc,cbs;
2459: PetscBool iscuda;
2463: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);
2464: ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);
2465: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);
2466: ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);
2467: VecDestroy(&is->x);
2468: VecDestroy(&is->y);
2469: VecDestroy(&is->counter);
2470: VecScatterDestroy(&is->rctx);
2471: VecScatterDestroy(&is->cctx);
2472: MatCreateVecs(is->A,&is->x,&is->y);
2473: PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);
2474: if (iscuda) {
2475: PetscFree(A->defaultvectype);
2476: PetscStrallocpy(VECCUDA,&A->defaultvectype);
2477: }
2478: MatCreateVecs(A,&cglobal,&rglobal);
2479: ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);
2480: ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
2481: VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
2482: ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);
2483: ISDestroy(&from);
2484: if (A->rmap->mapping != A->cmap->mapping) {
2485: ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);
2486: ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
2487: VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
2488: ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);
2489: ISDestroy(&from);
2490: } else {
2491: PetscObjectReference((PetscObject)is->rctx);
2492: is->cctx = is->rctx;
2493: }
2494: /* interface counter vector (local) */
2495: VecDuplicate(is->y,&is->counter);
2496: VecSet(is->y,1.);
2497: VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2498: VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2499: VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2500: VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2501: VecDestroy(&rglobal);
2502: VecDestroy(&cglobal);
2503: ISDestroy(&from);
2504: return(0);
2505: }
2507: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
2508: {
2510: PetscInt nr,rbs,nc,cbs;
2511: Mat_IS *is = (Mat_IS*)A->data;
2516: MatDestroy(&is->A);
2517: if (is->csf != is->sf) {
2518: PetscSFDestroy(&is->csf);
2519: PetscFree2(is->csf_rootdata,is->csf_leafdata);
2520: } else is->csf = NULL;
2521: PetscSFDestroy(&is->sf);
2522: PetscFree2(is->sf_rootdata,is->sf_leafdata);
2524: /* Setup Layout and set local to global maps */
2525: PetscLayoutSetUp(A->rmap);
2526: PetscLayoutSetUp(A->cmap);
2527: ISLocalToGlobalMappingGetSize(rmapping,&nr);
2528: ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
2529: ISLocalToGlobalMappingGetSize(cmapping,&nc);
2530: ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
2531: /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */
2532: if (rmapping != cmapping && A->rmap->N == A->cmap->N) {
2533: PetscBool same,gsame;
2535: same = PETSC_FALSE;
2536: if (nr == nc && cbs == rbs) {
2537: const PetscInt *idxs1,*idxs2;
2539: ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);
2540: ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);
2541: PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);
2542: ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);
2543: ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);
2544: }
2545: MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2546: if (gsame) cmapping = rmapping;
2547: }
2548: PetscLayoutSetBlockSize(A->rmap,rbs);
2549: PetscLayoutSetBlockSize(A->cmap,cbs);
2550: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
2551: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
2553: /* Create the local matrix A */
2554: MatCreate(PETSC_COMM_SELF,&is->A);
2555: MatSetType(is->A,is->lmattype);
2556: MatSetSizes(is->A,nr,nc,nr,nc);
2557: MatSetBlockSizes(is->A,rbs,cbs);
2558: MatSetOptionsPrefix(is->A,"is_");
2559: MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);
2560: PetscLayoutSetUp(is->A->rmap);
2561: PetscLayoutSetUp(is->A->cmap);
2563: if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
2564: MatISSetUpScatters_Private(A);
2565: }
2566: MatSetUp(A);
2567: return(0);
2568: }
2570: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2571: {
2572: Mat_IS *is = (Mat_IS*)mat->data;
2574: #if defined(PETSC_USE_DEBUG)
2575: PetscInt i,zm,zn;
2576: #endif
2577: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2580: #if defined(PETSC_USE_DEBUG)
2581: 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);
2582: /* count negative indices */
2583: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2584: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2585: #endif
2586: ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2587: ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2588: #if defined(PETSC_USE_DEBUG)
2589: /* count negative indices (should be the same as before) */
2590: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2591: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2592: 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");
2593: 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");
2594: #endif
2595: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
2596: return(0);
2597: }
2599: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2600: {
2601: Mat_IS *is = (Mat_IS*)mat->data;
2603: #if defined(PETSC_USE_DEBUG)
2604: PetscInt i,zm,zn;
2605: #endif
2606: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];
2609: #if defined(PETSC_USE_DEBUG)
2610: 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);
2611: /* count negative indices */
2612: for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2613: for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2614: #endif
2615: ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2616: ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2617: #if defined(PETSC_USE_DEBUG)
2618: /* count negative indices (should be the same as before) */
2619: for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2620: for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2621: 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");
2622: 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");
2623: #endif
2624: MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
2625: return(0);
2626: }
2628: static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2629: {
2631: Mat_IS *is = (Mat_IS*)A->data;
2634: if (is->A->rmap->mapping) {
2635: MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
2636: } else {
2637: MatSetValues(is->A,m,rows,n,cols,values,addv);
2638: }
2639: return(0);
2640: }
2642: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
2643: {
2645: Mat_IS *is = (Mat_IS*)A->data;
2648: if (is->A->rmap->mapping) {
2649: #if defined(PETSC_USE_DEBUG)
2650: PetscInt ibs,bs;
2652: ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);
2653: MatGetBlockSize(is->A,&bs);
2654: if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs);
2655: #endif
2656: MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
2657: } else {
2658: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
2659: }
2660: return(0);
2661: }
2663: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns)
2664: {
2665: Mat_IS *is = (Mat_IS*)A->data;
2669: if (!n) {
2670: is->pure_neumann = PETSC_TRUE;
2671: } else {
2672: PetscInt i;
2673: is->pure_neumann = PETSC_FALSE;
2675: if (columns) {
2676: MatZeroRowsColumns(is->A,n,rows,diag,0,0);
2677: } else {
2678: MatZeroRows(is->A,n,rows,diag,0,0);
2679: }
2680: if (diag != 0.) {
2681: const PetscScalar *array;
2682: VecGetArrayRead(is->counter,&array);
2683: for (i=0; i<n; i++) {
2684: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
2685: }
2686: VecRestoreArrayRead(is->counter,&array);
2687: }
2688: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
2689: MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
2690: }
2691: return(0);
2692: }
2694: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
2695: {
2696: Mat_IS *matis = (Mat_IS*)A->data;
2697: PetscInt nr,nl,len,i;
2698: PetscInt *lrows;
2702: #if defined(PETSC_USE_DEBUG)
2703: if (columns || diag != 0. || (x && b)) {
2704: PetscBool cong;
2706: PetscLayoutCompare(A->rmap,A->cmap,&cong);
2707: cong = (PetscBool)(cong && matis->sf == matis->csf);
2708: 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");
2709: 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");
2710: 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");
2711: }
2712: #endif
2713: /* get locally owned rows */
2714: PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);
2715: /* fix right hand side if needed */
2716: if (x && b) {
2717: const PetscScalar *xx;
2718: PetscScalar *bb;
2720: VecGetArrayRead(x, &xx);
2721: VecGetArray(b, &bb);
2722: for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
2723: VecRestoreArrayRead(x, &xx);
2724: VecRestoreArray(b, &bb);
2725: }
2726: /* get rows associated to the local matrices */
2727: MatISSetUpSF(A);
2728: MatGetSize(matis->A,&nl,NULL);
2729: PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));
2730: PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));
2731: for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
2732: PetscFree(lrows);
2733: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
2734: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
2735: PetscMalloc1(nl,&lrows);
2736: for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
2737: MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
2738: PetscFree(lrows);
2739: return(0);
2740: }
2742: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2743: {
2747: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
2748: return(0);
2749: }
2751: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2752: {
2756: MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
2757: return(0);
2758: }
2760: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
2761: {
2762: Mat_IS *is = (Mat_IS*)A->data;
2766: MatAssemblyBegin(is->A,type);
2767: return(0);
2768: }
2770: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
2771: {
2772: Mat_IS *is = (Mat_IS*)A->data;
2776: MatAssemblyEnd(is->A,type);
2777: /* fix for local empty rows/cols */
2778: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2779: Mat newlA;
2780: ISLocalToGlobalMapping rl2g,cl2g;
2781: IS nzr,nzc;
2782: PetscInt nr,nc,nnzr,nnzc;
2783: PetscBool lnewl2g,newl2g;
2785: MatGetSize(is->A,&nr,&nc);
2786: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
2787: if (!nzr) {
2788: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
2789: }
2790: MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
2791: if (!nzc) {
2792: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
2793: }
2794: ISGetSize(nzr,&nnzr);
2795: ISGetSize(nzc,&nnzc);
2796: if (nnzr != nr || nnzc != nc) {
2797: ISLocalToGlobalMapping l2g;
2798: IS is1,is2;
2800: /* need new global l2g map */
2801: lnewl2g = PETSC_TRUE;
2802: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2804: /* extract valid submatrix */
2805: MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);
2807: /* attach local l2g maps for successive calls of MatSetValues on the local matrix */
2808: ISLocalToGlobalMappingCreateIS(nzr,&l2g);
2809: ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);
2810: ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2811: ISLocalToGlobalMappingDestroy(&l2g);
2812: if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2813: const PetscInt *idxs1,*idxs2;
2814: PetscInt j,i,nl,*tidxs;
2816: ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);
2817: ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);
2818: PetscMalloc1(nl,&tidxs);
2819: ISGetIndices(is2,&idxs2);
2820: for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2821: if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr);
2822: ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);
2823: ISRestoreIndices(is2,&idxs2);
2824: ISDestroy(&is2);
2825: ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2826: }
2827: ISLocalToGlobalMappingCreateIS(is2,&rl2g);
2828: ISDestroy(&is1);
2829: ISDestroy(&is2);
2831: ISLocalToGlobalMappingCreateIS(nzc,&l2g);
2832: ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);
2833: ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2834: ISLocalToGlobalMappingDestroy(&l2g);
2835: if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2836: const PetscInt *idxs1,*idxs2;
2837: PetscInt j,i,nl,*tidxs;
2839: ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);
2840: ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);
2841: PetscMalloc1(nl,&tidxs);
2842: ISGetIndices(is2,&idxs2);
2843: for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2844: if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc);
2845: ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);
2846: ISRestoreIndices(is2,&idxs2);
2847: ISDestroy(&is2);
2848: ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2849: }
2850: ISLocalToGlobalMappingCreateIS(is2,&cl2g);
2851: ISDestroy(&is1);
2852: ISDestroy(&is2);
2854: MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);
2856: ISLocalToGlobalMappingDestroy(&rl2g);
2857: ISLocalToGlobalMappingDestroy(&cl2g);
2858: } else { /* local matrix fully populated */
2859: lnewl2g = PETSC_FALSE;
2860: MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2861: PetscObjectReference((PetscObject)is->A);
2862: newlA = is->A;
2863: }
2864: /* attach new global l2g map if needed */
2865: if (newl2g) {
2866: IS gnzr,gnzc;
2867: const PetscInt *grid,*gcid;
2869: ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);
2870: ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);
2871: ISGetIndices(gnzr,&grid);
2872: ISGetIndices(gnzc,&gcid);
2873: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);
2874: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);
2875: ISRestoreIndices(gnzr,&grid);
2876: ISRestoreIndices(gnzc,&gcid);
2877: ISDestroy(&gnzr);
2878: ISDestroy(&gnzc);
2879: MatSetLocalToGlobalMapping(A,rl2g,cl2g);
2880: ISLocalToGlobalMappingDestroy(&rl2g);
2881: ISLocalToGlobalMappingDestroy(&cl2g);
2882: }
2883: MatISSetLocalMat(A,newlA);
2884: MatDestroy(&newlA);
2885: ISDestroy(&nzr);
2886: ISDestroy(&nzc);
2887: is->locempty = PETSC_FALSE;
2888: }
2889: return(0);
2890: }
2892: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
2893: {
2894: Mat_IS *is = (Mat_IS*)mat->data;
2897: *local = is->A;
2898: return(0);
2899: }
2901: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
2902: {
2904: *local = NULL;
2905: return(0);
2906: }
2908: /*@
2909: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
2911: Input Parameter:
2912: . mat - the matrix
2914: Output Parameter:
2915: . local - the local matrix
2917: Level: advanced
2919: Notes:
2920: This can be called if you have precomputed the nonzero structure of the
2921: matrix and want to provide it to the inner matrix object to improve the performance
2922: of the MatSetValues() operation.
2924: Call MatISRestoreLocalMat() when finished with the local matrix.
2926: .seealso: MATIS
2927: @*/
2928: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
2929: {
2935: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
2936: return(0);
2937: }
2939: /*@
2940: MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()
2942: Input Parameter:
2943: . mat - the matrix
2945: Output Parameter:
2946: . local - the local matrix
2948: Level: advanced
2950: .seealso: MATIS
2951: @*/
2952: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
2953: {
2959: PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
2960: return(0);
2961: }
2963: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
2964: {
2965: Mat_IS *is = (Mat_IS*)mat->data;
2969: if (is->A) {
2970: MatSetType(is->A,mtype);
2971: }
2972: PetscFree(is->lmattype);
2973: PetscStrallocpy(mtype,&is->lmattype);
2974: return(0);
2975: }
2977: /*@
2978: MatISSetLocalMatType - Specifies the type of local matrix
2980: Input Parameter:
2981: . mat - the matrix
2982: . mtype - the local matrix type
2984: Output Parameter:
2986: Level: advanced
2988: .seealso: MATIS, MatSetType(), MatType
2989: @*/
2990: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
2991: {
2996: PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
2997: return(0);
2998: }
3000: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
3001: {
3002: Mat_IS *is = (Mat_IS*)mat->data;
3003: PetscInt nrows,ncols,orows,ocols;
3005: MatType mtype,otype;
3006: PetscBool sametype = PETSC_TRUE;
3009: if (is->A) {
3010: MatGetSize(is->A,&orows,&ocols);
3011: MatGetSize(local,&nrows,&ncols);
3012: 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);
3013: MatGetType(local,&mtype);
3014: MatGetType(is->A,&otype);
3015: PetscStrcmp(mtype,otype,&sametype);
3016: }
3017: PetscObjectReference((PetscObject)local);
3018: MatDestroy(&is->A);
3019: is->A = local;
3020: MatGetType(is->A,&mtype);
3021: MatISSetLocalMatType(mat,mtype);
3022: if (!sametype && !is->islocalref) {
3023: MatISSetUpScatters_Private(mat);
3024: }
3025: return(0);
3026: }
3028: /*@
3029: MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.
3031: Collective on Mat
3033: Input Parameter:
3034: . mat - the matrix
3035: . local - the local matrix
3037: Output Parameter:
3039: Level: advanced
3041: Notes:
3042: This can be called if you have precomputed the local matrix and
3043: want to provide it to the matrix object MATIS.
3045: .seealso: MATIS
3046: @*/
3047: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
3048: {
3054: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
3055: return(0);
3056: }
3058: static PetscErrorCode MatZeroEntries_IS(Mat A)
3059: {
3060: Mat_IS *a = (Mat_IS*)A->data;
3064: MatZeroEntries(a->A);
3065: return(0);
3066: }
3068: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
3069: {
3070: Mat_IS *is = (Mat_IS*)A->data;
3074: MatScale(is->A,a);
3075: return(0);
3076: }
3078: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3079: {
3080: Mat_IS *is = (Mat_IS*)A->data;
3084: /* get diagonal of the local matrix */
3085: MatGetDiagonal(is->A,is->y);
3087: /* scatter diagonal back into global vector */
3088: VecSet(v,0);
3089: VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3090: VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3091: return(0);
3092: }
3094: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
3095: {
3096: Mat_IS *a = (Mat_IS*)A->data;
3100: MatSetOption(a->A,op,flg);
3101: return(0);
3102: }
3104: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
3105: {
3106: Mat_IS *y = (Mat_IS*)Y->data;
3107: Mat_IS *x;
3108: #if defined(PETSC_USE_DEBUG)
3109: PetscBool ismatis;
3110: #endif
3114: #if defined(PETSC_USE_DEBUG)
3115: PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
3116: if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3117: #endif
3118: x = (Mat_IS*)X->data;
3119: MatAXPY(y->A,a,x->A,str);
3120: return(0);
3121: }
3123: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
3124: {
3125: Mat lA;
3126: Mat_IS *matis;
3127: ISLocalToGlobalMapping rl2g,cl2g;
3128: IS is;
3129: const PetscInt *rg,*rl;
3130: PetscInt nrg;
3131: PetscInt N,M,nrl,i,*idxs;
3132: PetscErrorCode ierr;
3135: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
3136: ISGetLocalSize(row,&nrl);
3137: ISGetIndices(row,&rl);
3138: ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
3139: #if defined(PETSC_USE_DEBUG)
3140: 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);
3141: #endif
3142: PetscMalloc1(nrg,&idxs);
3143: /* map from [0,nrl) to row */
3144: for (i=0;i<nrl;i++) idxs[i] = rl[i];
3145: #if defined(PETSC_USE_DEBUG)
3146: for (i=nrl;i<nrg;i++) idxs[i] = nrg;
3147: #else
3148: for (i=nrl;i<nrg;i++) idxs[i] = -1;
3149: #endif
3150: ISRestoreIndices(row,&rl);
3151: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
3152: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
3153: ISLocalToGlobalMappingCreateIS(is,&rl2g);
3154: ISDestroy(&is);
3155: /* compute new l2g map for columns */
3156: if (col != row || A->rmap->mapping != A->cmap->mapping) {
3157: const PetscInt *cg,*cl;
3158: PetscInt ncg;
3159: PetscInt ncl;
3161: ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
3162: ISGetLocalSize(col,&ncl);
3163: ISGetIndices(col,&cl);
3164: ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
3165: #if defined(PETSC_USE_DEBUG)
3166: 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);
3167: #endif
3168: PetscMalloc1(ncg,&idxs);
3169: /* map from [0,ncl) to col */
3170: for (i=0;i<ncl;i++) idxs[i] = cl[i];
3171: #if defined(PETSC_USE_DEBUG)
3172: for (i=ncl;i<ncg;i++) idxs[i] = ncg;
3173: #else
3174: for (i=ncl;i<ncg;i++) idxs[i] = -1;
3175: #endif
3176: ISRestoreIndices(col,&cl);
3177: ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
3178: ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
3179: ISLocalToGlobalMappingCreateIS(is,&cl2g);
3180: ISDestroy(&is);
3181: } else {
3182: PetscObjectReference((PetscObject)rl2g);
3183: cl2g = rl2g;
3184: }
3185: /* create the MATIS submatrix */
3186: MatGetSize(A,&M,&N);
3187: MatCreate(PetscObjectComm((PetscObject)A),submat);
3188: MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
3189: MatSetType(*submat,MATIS);
3190: matis = (Mat_IS*)((*submat)->data);
3191: matis->islocalref = PETSC_TRUE;
3192: MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
3193: MatISGetLocalMat(A,&lA);
3194: MatISSetLocalMat(*submat,lA);
3195: MatSetUp(*submat);
3196: MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
3197: MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
3198: ISLocalToGlobalMappingDestroy(&rl2g);
3199: ISLocalToGlobalMappingDestroy(&cl2g);
3200: /* remove unsupported ops */
3201: PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
3202: (*submat)->ops->destroy = MatDestroy_IS;
3203: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3204: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3205: (*submat)->ops->assemblybegin = MatAssemblyBegin_IS;
3206: (*submat)->ops->assemblyend = MatAssemblyEnd_IS;
3207: return(0);
3208: }
3210: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
3211: {
3212: Mat_IS *a = (Mat_IS*)A->data;
3213: char type[256];
3214: PetscBool flg;
3218: PetscOptionsHead(PetscOptionsObject,"MATIS options");
3219: PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
3220: PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
3221: PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
3222: if (flg) {
3223: MatISSetLocalMatType(A,type);
3224: }
3225: if (a->A) {
3226: MatSetFromOptions(a->A);
3227: }
3228: PetscOptionsTail();
3229: return(0);
3230: }
3232: /*@
3233: MatCreateIS - Creates a "process" unassembled matrix, assembled on each
3234: process but not across processes.
3236: Input Parameters:
3237: + comm - MPI communicator that will share the matrix
3238: . bs - block size of the matrix
3239: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3240: . rmap - local to global map for rows
3241: - cmap - local to global map for cols
3243: Output Parameter:
3244: . A - the resulting matrix
3246: Level: advanced
3248: Notes:
3249: See MATIS for more details.
3250: m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors
3251: used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
3252: If either rmap or cmap are NULL, then the matrix is assumed to be square.
3254: .seealso: MATIS, MatSetLocalToGlobalMapping()
3255: @*/
3256: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
3257: {
3261: if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
3262: MatCreate(comm,A);
3263: MatSetSizes(*A,m,n,M,N);
3264: if (bs > 0) {
3265: MatSetBlockSize(*A,bs);
3266: }
3267: MatSetType(*A,MATIS);
3268: if (rmap && cmap) {
3269: MatSetLocalToGlobalMapping(*A,rmap,cmap);
3270: } else if (!rmap) {
3271: MatSetLocalToGlobalMapping(*A,cmap,cmap);
3272: } else {
3273: MatSetLocalToGlobalMapping(*A,rmap,rmap);
3274: }
3275: return(0);
3276: }
3278: /*MC
3279: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP).
3280: This stores the matrices in globally unassembled form. Each processor
3281: assembles only its local Neumann problem and the parallel matrix vector
3282: product is handled "implicitly".
3284: Operations Provided:
3285: + MatMult()
3286: . MatMultAdd()
3287: . MatMultTranspose()
3288: . MatMultTransposeAdd()
3289: . MatZeroEntries()
3290: . MatSetOption()
3291: . MatZeroRows()
3292: . MatSetValues()
3293: . MatSetValuesBlocked()
3294: . MatSetValuesLocal()
3295: . MatSetValuesBlockedLocal()
3296: . MatScale()
3297: . MatGetDiagonal()
3298: . MatMissingDiagonal()
3299: . MatDuplicate()
3300: . MatCopy()
3301: . MatAXPY()
3302: . MatCreateSubMatrix()
3303: . MatGetLocalSubMatrix()
3304: . MatTranspose()
3305: . MatPtAP() (with P of AIJ type)
3306: - MatSetLocalToGlobalMapping()
3308: Options Database Keys:
3309: + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
3310: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3311: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP().
3313: Notes:
3314: Options prefix for the inner matrix are given by -is_mat_xxx
3316: You must call MatSetLocalToGlobalMapping() before using this matrix type.
3318: You can do matrix preallocation on the local matrix after you obtain it with
3319: MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation()
3321: Level: advanced
3323: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP
3325: M*/
3327: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3328: {
3330: Mat_IS *b;
3333: PetscNewLog(A,&b);
3334: PetscStrallocpy(MATAIJ,&b->lmattype);
3335: A->data = (void*)b;
3337: /* matrix ops */
3338: PetscMemzero(A->ops,sizeof(struct _MatOps));
3339: A->ops->mult = MatMult_IS;
3340: A->ops->multadd = MatMultAdd_IS;
3341: A->ops->multtranspose = MatMultTranspose_IS;
3342: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3343: A->ops->destroy = MatDestroy_IS;
3344: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3345: A->ops->setvalues = MatSetValues_IS;
3346: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3347: A->ops->setvalueslocal = MatSetValuesLocal_IS;
3348: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3349: A->ops->zerorows = MatZeroRows_IS;
3350: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3351: A->ops->assemblybegin = MatAssemblyBegin_IS;
3352: A->ops->assemblyend = MatAssemblyEnd_IS;
3353: A->ops->view = MatView_IS;
3354: A->ops->zeroentries = MatZeroEntries_IS;
3355: A->ops->scale = MatScale_IS;
3356: A->ops->getdiagonal = MatGetDiagonal_IS;
3357: A->ops->setoption = MatSetOption_IS;
3358: A->ops->ishermitian = MatIsHermitian_IS;
3359: A->ops->issymmetric = MatIsSymmetric_IS;
3360: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3361: A->ops->duplicate = MatDuplicate_IS;
3362: A->ops->missingdiagonal = MatMissingDiagonal_IS;
3363: A->ops->copy = MatCopy_IS;
3364: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3365: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3366: A->ops->axpy = MatAXPY_IS;
3367: A->ops->diagonalset = MatDiagonalSet_IS;
3368: A->ops->shift = MatShift_IS;
3369: A->ops->transpose = MatTranspose_IS;
3370: A->ops->getinfo = MatGetInfo_IS;
3371: A->ops->diagonalscale = MatDiagonalScale_IS;
3372: A->ops->setfromoptions = MatSetFromOptions_IS;
3374: /* special MATIS functions */
3375: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
3376: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
3377: PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
3378: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
3379: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
3380: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
3381: PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);
3382: PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
3383: PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
3384: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
3385: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
3386: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
3387: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
3388: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
3389: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
3390: PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
3391: PetscObjectChangeTypeName((PetscObject)A,MATIS);
3392: return(0);
3393: }