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