Actual source code: matnest.c
petsc-3.7.5 2017-01-01
2: #include <../src/mat/impls/nest/matnestimpl.h> /*I "petscmat.h" I*/
3: #include <petscsf.h>
5: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
6: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left);
8: /* private functions */
11: static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
12: {
13: Mat_Nest *bA = (Mat_Nest*)A->data;
14: PetscInt i,j;
18: *m = *n = *M = *N = 0;
19: for (i=0; i<bA->nr; i++) { /* rows */
20: PetscInt sm,sM;
21: ISGetLocalSize(bA->isglobal.row[i],&sm);
22: ISGetSize(bA->isglobal.row[i],&sM);
23: *m += sm;
24: *M += sM;
25: }
26: for (j=0; j<bA->nc; j++) { /* cols */
27: PetscInt sn,sN;
28: ISGetLocalSize(bA->isglobal.col[j],&sn);
29: ISGetSize(bA->isglobal.col[j],&sN);
30: *n += sn;
31: *N += sN;
32: }
33: return(0);
34: }
36: /* operations */
39: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
40: {
41: Mat_Nest *bA = (Mat_Nest*)A->data;
42: Vec *bx = bA->right,*by = bA->left;
43: PetscInt i,j,nr = bA->nr,nc = bA->nc;
47: for (i=0; i<nr; i++) {VecGetSubVector(y,bA->isglobal.row[i],&by[i]);}
48: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
49: for (i=0; i<nr; i++) {
50: VecZeroEntries(by[i]);
51: for (j=0; j<nc; j++) {
52: if (!bA->m[i][j]) continue;
53: /* y[i] <- y[i] + A[i][j] * x[j] */
54: MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);
55: }
56: }
57: for (i=0; i<nr; i++) {VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);}
58: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
59: return(0);
60: }
64: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
65: {
66: Mat_Nest *bA = (Mat_Nest*)A->data;
67: Vec *bx = bA->right,*bz = bA->left;
68: PetscInt i,j,nr = bA->nr,nc = bA->nc;
72: for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
73: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
74: for (i=0; i<nr; i++) {
75: if (y != z) {
76: Vec by;
77: VecGetSubVector(y,bA->isglobal.row[i],&by);
78: VecCopy(by,bz[i]);
79: VecRestoreSubVector(y,bA->isglobal.row[i],&by);
80: }
81: for (j=0; j<nc; j++) {
82: if (!bA->m[i][j]) continue;
83: /* y[i] <- y[i] + A[i][j] * x[j] */
84: MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
85: }
86: }
87: for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
88: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
89: return(0);
90: }
94: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
95: {
96: Mat_Nest *bA = (Mat_Nest*)A->data;
97: Vec *bx = bA->left,*by = bA->right;
98: PetscInt i,j,nr = bA->nr,nc = bA->nc;
102: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
103: for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
104: for (j=0; j<nc; j++) {
105: VecZeroEntries(by[j]);
106: for (i=0; i<nr; i++) {
107: if (!bA->m[i][j]) continue;
108: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
109: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
110: }
111: }
112: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
113: for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
114: return(0);
115: }
119: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
120: {
121: Mat_Nest *bA = (Mat_Nest*)A->data;
122: Vec *bx = bA->left,*bz = bA->right;
123: PetscInt i,j,nr = bA->nr,nc = bA->nc;
127: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
128: for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
129: for (j=0; j<nc; j++) {
130: if (y != z) {
131: Vec by;
132: VecGetSubVector(y,bA->isglobal.col[j],&by);
133: VecCopy(by,bz[j]);
134: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
135: }
136: for (i=0; i<nr; i++) {
137: if (!bA->m[i][j]) continue;
138: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
139: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
140: }
141: }
142: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
143: for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
144: return(0);
145: }
149: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
150: {
152: IS *lst = *list;
153: PetscInt i;
156: if (!lst) return(0);
157: for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
158: PetscFree(lst);
159: *list = NULL;
160: return(0);
161: }
165: static PetscErrorCode MatDestroy_Nest(Mat A)
166: {
167: Mat_Nest *vs = (Mat_Nest*)A->data;
168: PetscInt i,j;
172: /* release the matrices and the place holders */
173: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
174: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
175: MatNestDestroyISList(vs->nr,&vs->islocal.row);
176: MatNestDestroyISList(vs->nc,&vs->islocal.col);
178: PetscFree(vs->row_len);
179: PetscFree(vs->col_len);
181: PetscFree2(vs->left,vs->right);
183: /* release the matrices and the place holders */
184: if (vs->m) {
185: for (i=0; i<vs->nr; i++) {
186: for (j=0; j<vs->nc; j++) {
187: MatDestroy(&vs->m[i][j]);
188: }
189: PetscFree(vs->m[i]);
190: }
191: PetscFree(vs->m);
192: }
193: PetscFree(A->data);
195: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
196: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
197: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
198: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
199: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
200: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
201: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
202: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
203: return(0);
204: }
208: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
209: {
210: Mat_Nest *vs = (Mat_Nest*)A->data;
211: PetscInt i,j;
215: for (i=0; i<vs->nr; i++) {
216: for (j=0; j<vs->nc; j++) {
217: if (vs->m[i][j]) {
218: MatAssemblyBegin(vs->m[i][j],type);
219: if (!vs->splitassembly) {
220: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
221: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
222: * already performing an assembly, but the result would by more complicated and appears to offer less
223: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
224: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
225: */
226: MatAssemblyEnd(vs->m[i][j],type);
227: }
228: }
229: }
230: }
231: return(0);
232: }
236: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
237: {
238: Mat_Nest *vs = (Mat_Nest*)A->data;
239: PetscInt i,j;
243: for (i=0; i<vs->nr; i++) {
244: for (j=0; j<vs->nc; j++) {
245: if (vs->m[i][j]) {
246: if (vs->splitassembly) {
247: MatAssemblyEnd(vs->m[i][j],type);
248: }
249: }
250: }
251: }
252: return(0);
253: }
257: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
258: {
260: Mat_Nest *vs = (Mat_Nest*)A->data;
261: PetscInt j;
262: Mat sub;
265: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
266: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
267: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
268: *B = sub;
269: return(0);
270: }
274: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
275: {
277: Mat_Nest *vs = (Mat_Nest*)A->data;
278: PetscInt i;
279: Mat sub;
282: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
283: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
284: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
285: *B = sub;
286: return(0);
287: }
291: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
292: {
294: PetscInt i;
295: PetscBool flg;
301: *found = -1;
302: for (i=0; i<n; i++) {
303: if (!list[i]) continue;
304: ISEqual(list[i],is,&flg);
305: if (flg) {
306: *found = i;
307: return(0);
308: }
309: }
310: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
311: return(0);
312: }
316: /* Get a block row as a new MatNest */
317: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
318: {
319: Mat_Nest *vs = (Mat_Nest*)A->data;
320: char keyname[256];
324: *B = NULL;
325: PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
326: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
327: if (*B) return(0);
329: MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);
331: (*B)->assembled = A->assembled;
333: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
334: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
335: return(0);
336: }
340: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
341: {
342: Mat_Nest *vs = (Mat_Nest*)A->data;
344: PetscInt row,col;
345: PetscBool same,isFullCol,isFullColGlobal;
348: /* Check if full column space. This is a hack */
349: isFullCol = PETSC_FALSE;
350: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
351: if (same) {
352: PetscInt n,first,step,i,an,am,afirst,astep;
353: ISStrideGetInfo(iscol,&first,&step);
354: ISGetLocalSize(iscol,&n);
355: isFullCol = PETSC_TRUE;
356: for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
357: ISStrideGetInfo(is->col[i],&afirst,&astep);
358: ISGetLocalSize(is->col[i],&am);
359: if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
360: an += am;
361: }
362: if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
363: }
364: MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));
366: if (isFullColGlobal) {
367: PetscInt row;
368: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
369: MatNestGetRow(A,row,B);
370: } else {
371: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
372: MatNestFindIS(A,vs->nc,is->col,iscol,&col);
373: if (!vs->m[row][col]) {
374: PetscInt lr,lc;
376: MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
377: ISGetLocalSize(vs->isglobal.row[row],&lr);
378: ISGetLocalSize(vs->isglobal.col[col],&lc);
379: MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
380: MatSetUp(vs->m[row][col]);
381: MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
382: MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
383: }
384: *B = vs->m[row][col];
385: }
386: return(0);
387: }
391: static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
392: {
394: Mat_Nest *vs = (Mat_Nest*)A->data;
395: Mat sub;
398: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
399: switch (reuse) {
400: case MAT_INITIAL_MATRIX:
401: if (sub) { PetscObjectReference((PetscObject)sub); }
402: *B = sub;
403: break;
404: case MAT_REUSE_MATRIX:
405: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
406: break;
407: case MAT_IGNORE_MATRIX: /* Nothing to do */
408: break;
409: case MAT_INPLACE_MATRIX: /* Nothing to do */
410: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
411: break;
412: }
413: return(0);
414: }
418: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
419: {
421: Mat_Nest *vs = (Mat_Nest*)A->data;
422: Mat sub;
425: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
426: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
427: if (sub) {PetscObjectReference((PetscObject)sub);}
428: *B = sub;
429: return(0);
430: }
434: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
435: {
437: Mat_Nest *vs = (Mat_Nest*)A->data;
438: Mat sub;
441: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
442: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
443: if (sub) {
444: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
445: MatDestroy(B);
446: }
447: return(0);
448: }
452: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
453: {
454: Mat_Nest *bA = (Mat_Nest*)A->data;
455: PetscInt i;
459: for (i=0; i<bA->nr; i++) {
460: Vec bv;
461: VecGetSubVector(v,bA->isglobal.row[i],&bv);
462: if (bA->m[i][i]) {
463: MatGetDiagonal(bA->m[i][i],bv);
464: } else {
465: VecSet(bv,0.0);
466: }
467: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
468: }
469: return(0);
470: }
474: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
475: {
476: Mat_Nest *bA = (Mat_Nest*)A->data;
477: Vec bl,*br;
478: PetscInt i,j;
482: PetscCalloc1(bA->nc,&br);
483: if (r) {
484: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
485: }
486: bl = NULL;
487: for (i=0; i<bA->nr; i++) {
488: if (l) {
489: VecGetSubVector(l,bA->isglobal.row[i],&bl);
490: }
491: for (j=0; j<bA->nc; j++) {
492: if (bA->m[i][j]) {
493: MatDiagonalScale(bA->m[i][j],bl,br[j]);
494: }
495: }
496: if (l) {
497: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
498: }
499: }
500: if (r) {
501: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
502: }
503: PetscFree(br);
504: return(0);
505: }
509: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
510: {
511: Mat_Nest *bA = (Mat_Nest*)A->data;
512: PetscInt i,j;
516: for (i=0; i<bA->nr; i++) {
517: for (j=0; j<bA->nc; j++) {
518: if (bA->m[i][j]) {
519: MatScale(bA->m[i][j],a);
520: }
521: }
522: }
523: return(0);
524: }
528: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
529: {
530: Mat_Nest *bA = (Mat_Nest*)A->data;
531: PetscInt i;
535: for (i=0; i<bA->nr; i++) {
536: if (!bA->m[i][i]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for shifting an empty diagonal block, insert a matrix in block (%D,%D)",i,i);
537: MatShift(bA->m[i][i],a);
538: }
539: return(0);
540: }
544: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
545: {
546: Mat_Nest *bA = (Mat_Nest*)A->data;
547: Vec *L,*R;
548: MPI_Comm comm;
549: PetscInt i,j;
553: PetscObjectGetComm((PetscObject)A,&comm);
554: if (right) {
555: /* allocate R */
556: PetscMalloc1(bA->nc, &R);
557: /* Create the right vectors */
558: for (j=0; j<bA->nc; j++) {
559: for (i=0; i<bA->nr; i++) {
560: if (bA->m[i][j]) {
561: MatCreateVecs(bA->m[i][j],&R[j],NULL);
562: break;
563: }
564: }
565: if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
566: }
567: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
568: /* hand back control to the nest vector */
569: for (j=0; j<bA->nc; j++) {
570: VecDestroy(&R[j]);
571: }
572: PetscFree(R);
573: }
575: if (left) {
576: /* allocate L */
577: PetscMalloc1(bA->nr, &L);
578: /* Create the left vectors */
579: for (i=0; i<bA->nr; i++) {
580: for (j=0; j<bA->nc; j++) {
581: if (bA->m[i][j]) {
582: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
583: break;
584: }
585: }
586: if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
587: }
589: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
590: for (i=0; i<bA->nr; i++) {
591: VecDestroy(&L[i]);
592: }
594: PetscFree(L);
595: }
596: return(0);
597: }
601: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
602: {
603: Mat_Nest *bA = (Mat_Nest*)A->data;
604: PetscBool isascii;
605: PetscInt i,j;
609: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
610: if (isascii) {
612: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
613: PetscViewerASCIIPushTab(viewer); /* push0 */
614: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc);
616: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
617: for (i=0; i<bA->nr; i++) {
618: for (j=0; j<bA->nc; j++) {
619: MatType type;
620: char name[256] = "",prefix[256] = "";
621: PetscInt NR,NC;
622: PetscBool isNest = PETSC_FALSE;
624: if (!bA->m[i][j]) {
625: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
626: continue;
627: }
628: MatGetSize(bA->m[i][j],&NR,&NC);
629: MatGetType(bA->m[i][j], &type);
630: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
631: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
632: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
634: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
636: if (isNest) {
637: PetscViewerASCIIPushTab(viewer); /* push1 */
638: MatView(bA->m[i][j],viewer);
639: PetscViewerASCIIPopTab(viewer); /* pop1 */
640: }
641: }
642: }
643: PetscViewerASCIIPopTab(viewer); /* pop0 */
644: }
645: return(0);
646: }
650: static PetscErrorCode MatZeroEntries_Nest(Mat A)
651: {
652: Mat_Nest *bA = (Mat_Nest*)A->data;
653: PetscInt i,j;
657: for (i=0; i<bA->nr; i++) {
658: for (j=0; j<bA->nc; j++) {
659: if (!bA->m[i][j]) continue;
660: MatZeroEntries(bA->m[i][j]);
661: }
662: }
663: return(0);
664: }
668: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
669: {
670: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
671: PetscInt i,j,nr = bA->nr,nc = bA->nc;
675: if (nr != bB->nr || nc != bB->nc) SETERRQ4(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot copy a Mat_Nest of block size (%D,%D) to a Mat_Nest of block size (%D,%D)",bB->nr,bB->nc,nr,nc);
676: for (i=0; i<nr; i++) {
677: for (j=0; j<nc; j++) {
678: if (bA->m[i][j] && bB->m[i][j]) {
679: MatCopy(bA->m[i][j],bB->m[i][j],str);
680: } else if (bA->m[i][j] || bB->m[i][j]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j);
681: }
682: }
683: return(0);
684: }
688: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
689: {
690: Mat_Nest *bA = (Mat_Nest*)A->data;
691: Mat *b;
692: PetscInt i,j,nr = bA->nr,nc = bA->nc;
696: PetscMalloc1(nr*nc,&b);
697: for (i=0; i<nr; i++) {
698: for (j=0; j<nc; j++) {
699: if (bA->m[i][j]) {
700: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
701: } else {
702: b[i*nc+j] = NULL;
703: }
704: }
705: }
706: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
707: /* Give the new MatNest exclusive ownership */
708: for (i=0; i<nr*nc; i++) {
709: MatDestroy(&b[i]);
710: }
711: PetscFree(b);
713: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
714: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
715: return(0);
716: }
718: /* nest api */
721: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
722: {
723: Mat_Nest *bA = (Mat_Nest*)A->data;
726: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
727: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
728: *mat = bA->m[idxm][jdxm];
729: return(0);
730: }
734: /*@
735: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
737: Not collective
739: Input Parameters:
740: + A - nest matrix
741: . idxm - index of the matrix within the nest matrix
742: - jdxm - index of the matrix within the nest matrix
744: Output Parameter:
745: . sub - matrix at index idxm,jdxm within the nest matrix
747: Level: developer
749: .seealso: MatNestGetSize(), MatNestGetSubMats()
750: @*/
751: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
752: {
756: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
757: return(0);
758: }
762: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
763: {
764: Mat_Nest *bA = (Mat_Nest*)A->data;
765: PetscInt m,n,M,N,mi,ni,Mi,Ni;
769: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
770: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
771: MatGetLocalSize(mat,&m,&n);
772: MatGetSize(mat,&M,&N);
773: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
774: ISGetSize(bA->isglobal.row[idxm],&Mi);
775: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
776: ISGetSize(bA->isglobal.col[jdxm],&Ni);
777: if (M != Mi || N != Ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix dimension (%D,%D) incompatible with nest block (%D,%D)",M,N,Mi,Ni);
778: if (m != mi || n != ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix local dimension (%D,%D) incompatible with nest block (%D,%D)",m,n,mi,ni);
780: PetscObjectReference((PetscObject)mat);
781: MatDestroy(&bA->m[idxm][jdxm]);
782: bA->m[idxm][jdxm] = mat;
783: return(0);
784: }
788: /*@
789: MatNestSetSubMat - Set a single submatrix in the nest matrix.
791: Logically collective on the submatrix communicator
793: Input Parameters:
794: + A - nest matrix
795: . idxm - index of the matrix within the nest matrix
796: . jdxm - index of the matrix within the nest matrix
797: - sub - matrix at index idxm,jdxm within the nest matrix
799: Notes:
800: The new submatrix must have the same size and communicator as that block of the nest.
802: This increments the reference count of the submatrix.
804: Level: developer
806: .seealso: MatNestSetSubMats(), MatNestGetSubMat()
807: @*/
808: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
809: {
813: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
814: return(0);
815: }
819: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
820: {
821: Mat_Nest *bA = (Mat_Nest*)A->data;
824: if (M) *M = bA->nr;
825: if (N) *N = bA->nc;
826: if (mat) *mat = bA->m;
827: return(0);
828: }
832: /*@C
833: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
835: Not collective
837: Input Parameters:
838: . A - nest matrix
840: Output Parameter:
841: + M - number of rows in the nest matrix
842: . N - number of cols in the nest matrix
843: - mat - 2d array of matrices
845: Notes:
847: The user should not free the array mat.
849: Level: developer
851: .seealso: MatNestGetSize(), MatNestGetSubMat()
852: @*/
853: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
854: {
858: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
859: return(0);
860: }
864: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
865: {
866: Mat_Nest *bA = (Mat_Nest*)A->data;
869: if (M) *M = bA->nr;
870: if (N) *N = bA->nc;
871: return(0);
872: }
876: /*@
877: MatNestGetSize - Returns the size of the nest matrix.
879: Not collective
881: Input Parameters:
882: . A - nest matrix
884: Output Parameter:
885: + M - number of rows in the nested mat
886: - N - number of cols in the nested mat
888: Notes:
890: Level: developer
892: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
893: @*/
894: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
895: {
899: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
900: return(0);
901: }
905: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
906: {
907: Mat_Nest *vs = (Mat_Nest*)A->data;
908: PetscInt i;
911: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
912: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
913: return(0);
914: }
918: /*@C
919: MatNestGetISs - Returns the index sets partitioning the row and column spaces
921: Not collective
923: Input Parameters:
924: . A - nest matrix
926: Output Parameter:
927: + rows - array of row index sets
928: - cols - array of column index sets
930: Level: advanced
932: Notes:
933: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
935: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
936: @*/
937: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
938: {
943: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
944: return(0);
945: }
949: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
950: {
951: Mat_Nest *vs = (Mat_Nest*)A->data;
952: PetscInt i;
955: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
956: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
957: return(0);
958: }
962: /*@C
963: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
965: Not collective
967: Input Parameters:
968: . A - nest matrix
970: Output Parameter:
971: + rows - array of row index sets (or NULL to ignore)
972: - cols - array of column index sets (or NULL to ignore)
974: Level: advanced
976: Notes:
977: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
979: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
980: @*/
981: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
982: {
987: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
988: return(0);
989: }
993: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
994: {
996: PetscBool flg;
999: PetscStrcmp(vtype,VECNEST,&flg);
1000: /* In reality, this only distinguishes VECNEST and "other" */
1001: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1002: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1003: return(0);
1004: }
1008: /*@C
1009: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1011: Not collective
1013: Input Parameters:
1014: + A - nest matrix
1015: - vtype - type to use for creating vectors
1017: Notes:
1019: Level: developer
1021: .seealso: MatCreateVecs()
1022: @*/
1023: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1024: {
1028: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1029: return(0);
1030: }
1034: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1035: {
1036: Mat_Nest *s = (Mat_Nest*)A->data;
1037: PetscInt i,j,m,n,M,N;
1041: s->nr = nr;
1042: s->nc = nc;
1044: /* Create space for submatrices */
1045: PetscMalloc1(nr,&s->m);
1046: for (i=0; i<nr; i++) {
1047: PetscMalloc1(nc,&s->m[i]);
1048: }
1049: for (i=0; i<nr; i++) {
1050: for (j=0; j<nc; j++) {
1051: s->m[i][j] = a[i*nc+j];
1052: if (a[i*nc+j]) {
1053: PetscObjectReference((PetscObject)a[i*nc+j]);
1054: }
1055: }
1056: }
1058: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1060: PetscMalloc1(nr,&s->row_len);
1061: PetscMalloc1(nc,&s->col_len);
1062: for (i=0; i<nr; i++) s->row_len[i]=-1;
1063: for (j=0; j<nc; j++) s->col_len[j]=-1;
1065: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1067: PetscLayoutSetSize(A->rmap,M);
1068: PetscLayoutSetLocalSize(A->rmap,m);
1069: PetscLayoutSetSize(A->cmap,N);
1070: PetscLayoutSetLocalSize(A->cmap,n);
1072: PetscLayoutSetUp(A->rmap);
1073: PetscLayoutSetUp(A->cmap);
1075: PetscCalloc2(nr,&s->left,nc,&s->right);
1076: return(0);
1077: }
1081: /*@
1082: MatNestSetSubMats - Sets the nested submatrices
1084: Collective on Mat
1086: Input Parameter:
1087: + N - nested matrix
1088: . nr - number of nested row blocks
1089: . is_row - index sets for each nested row block, or NULL to make contiguous
1090: . nc - number of nested column blocks
1091: . is_col - index sets for each nested column block, or NULL to make contiguous
1092: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1094: Level: advanced
1096: .seealso: MatCreateNest(), MATNEST
1097: @*/
1098: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1099: {
1101: PetscInt i;
1105: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1106: if (nr && is_row) {
1109: }
1110: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1111: if (nc && is_col) {
1114: }
1116: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1117: return(0);
1118: }
1122: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1123: {
1125: PetscBool flg;
1126: PetscInt i,j,m,mi,*ix;
1129: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1130: if (islocal[i]) {
1131: ISGetSize(islocal[i],&mi);
1132: flg = PETSC_TRUE; /* We found a non-trivial entry */
1133: } else {
1134: ISGetSize(isglobal[i],&mi);
1135: }
1136: m += mi;
1137: }
1138: if (flg) {
1139: PetscMalloc1(m,&ix);
1140: for (i=0,n=0; i<n; i++) {
1141: ISLocalToGlobalMapping smap = NULL;
1142: VecScatter scat;
1143: IS isreq;
1144: Vec lvec,gvec;
1145: union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x;
1146: Mat sub;
1148: if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers");
1149: if (colflg) {
1150: MatNestFindNonzeroSubMatRow(A,i,&sub);
1151: } else {
1152: MatNestFindNonzeroSubMatCol(A,i,&sub);
1153: }
1154: if (sub) {MatGetLocalToGlobalMapping(sub,&smap,NULL);}
1155: if (islocal[i]) {
1156: ISGetSize(islocal[i],&mi);
1157: } else {
1158: ISGetSize(isglobal[i],&mi);
1159: }
1160: for (j=0; j<mi; j++) ix[m+j] = j;
1161: if (smap) {ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);}
1162: /*
1163: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1164: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1165: The approach here is ugly because it uses VecScatter to move indices.
1166: */
1167: VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);
1168: VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);
1169: ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);
1170: VecScatterCreate(gvec,isreq,lvec,NULL,&scat);
1171: VecGetArray(gvec,(PetscScalar**)&x);
1172: for (j=0; j<mi; j++) x[j].integer = ix[m+j];
1173: VecRestoreArray(gvec,(PetscScalar**)&x);
1174: VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1175: VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1176: VecGetArray(lvec,(PetscScalar**)&x);
1177: for (j=0; j<mi; j++) ix[m+j] = x[j].integer;
1178: VecRestoreArray(lvec,(PetscScalar**)&x);
1179: VecDestroy(&lvec);
1180: VecDestroy(&gvec);
1181: ISDestroy(&isreq);
1182: VecScatterDestroy(&scat);
1183: m += mi;
1184: }
1185: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1186: } else {
1187: *ltog = NULL;
1188: }
1189: return(0);
1190: }
1193: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1194: /*
1195: nprocessors = NP
1196: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1197: proc 0: => (g_0,h_0,)
1198: proc 1: => (g_1,h_1,)
1199: ...
1200: proc nprocs-1: => (g_NP-1,h_NP-1,)
1202: proc 0: proc 1: proc nprocs-1:
1203: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1205: proc 0:
1206: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1207: proc 1:
1208: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1210: proc NP-1:
1211: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1212: */
1215: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1216: {
1217: Mat_Nest *vs = (Mat_Nest*)A->data;
1218: PetscInt i,j,offset,n,nsum,bs;
1220: Mat sub = NULL;
1223: PetscMalloc1(nr,&vs->isglobal.row);
1224: PetscMalloc1(nc,&vs->isglobal.col);
1225: if (is_row) { /* valid IS is passed in */
1226: /* refs on is[] are incremeneted */
1227: for (i=0; i<vs->nr; i++) {
1228: PetscObjectReference((PetscObject)is_row[i]);
1230: vs->isglobal.row[i] = is_row[i];
1231: }
1232: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1233: nsum = 0;
1234: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1235: MatNestFindNonzeroSubMatRow(A,i,&sub);
1236: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1237: MatGetLocalSize(sub,&n,NULL);
1238: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1239: nsum += n;
1240: }
1241: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1242: offset -= nsum;
1243: for (i=0; i<vs->nr; i++) {
1244: MatNestFindNonzeroSubMatRow(A,i,&sub);
1245: MatGetLocalSize(sub,&n,NULL);
1246: MatGetBlockSize(sub,&bs);
1247: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1248: ISSetBlockSize(vs->isglobal.row[i],bs);
1249: offset += n;
1250: }
1251: }
1253: if (is_col) { /* valid IS is passed in */
1254: /* refs on is[] are incremeneted */
1255: for (j=0; j<vs->nc; j++) {
1256: PetscObjectReference((PetscObject)is_col[j]);
1258: vs->isglobal.col[j] = is_col[j];
1259: }
1260: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1261: offset = A->cmap->rstart;
1262: nsum = 0;
1263: for (j=0; j<vs->nc; j++) {
1264: MatNestFindNonzeroSubMatCol(A,j,&sub);
1265: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1266: MatGetLocalSize(sub,NULL,&n);
1267: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1268: nsum += n;
1269: }
1270: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1271: offset -= nsum;
1272: for (j=0; j<vs->nc; j++) {
1273: MatNestFindNonzeroSubMatCol(A,j,&sub);
1274: MatGetLocalSize(sub,NULL,&n);
1275: MatGetBlockSize(sub,&bs);
1276: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1277: ISSetBlockSize(vs->isglobal.col[j],bs);
1278: offset += n;
1279: }
1280: }
1282: /* Set up the local ISs */
1283: PetscMalloc1(vs->nr,&vs->islocal.row);
1284: PetscMalloc1(vs->nc,&vs->islocal.col);
1285: for (i=0,offset=0; i<vs->nr; i++) {
1286: IS isloc;
1287: ISLocalToGlobalMapping rmap = NULL;
1288: PetscInt nlocal,bs;
1289: MatNestFindNonzeroSubMatRow(A,i,&sub);
1290: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1291: if (rmap) {
1292: MatGetBlockSize(sub,&bs);
1293: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1294: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1295: ISSetBlockSize(isloc,bs);
1296: } else {
1297: nlocal = 0;
1298: isloc = NULL;
1299: }
1300: vs->islocal.row[i] = isloc;
1301: offset += nlocal;
1302: }
1303: for (i=0,offset=0; i<vs->nc; i++) {
1304: IS isloc;
1305: ISLocalToGlobalMapping cmap = NULL;
1306: PetscInt nlocal,bs;
1307: MatNestFindNonzeroSubMatCol(A,i,&sub);
1308: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1309: if (cmap) {
1310: MatGetBlockSize(sub,&bs);
1311: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1312: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1313: ISSetBlockSize(isloc,bs);
1314: } else {
1315: nlocal = 0;
1316: isloc = NULL;
1317: }
1318: vs->islocal.col[i] = isloc;
1319: offset += nlocal;
1320: }
1322: /* Set up the aggregate ISLocalToGlobalMapping */
1323: {
1324: ISLocalToGlobalMapping rmap,cmap;
1325: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1326: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1327: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1328: ISLocalToGlobalMappingDestroy(&rmap);
1329: ISLocalToGlobalMappingDestroy(&cmap);
1330: }
1332: #if defined(PETSC_USE_DEBUG)
1333: for (i=0; i<vs->nr; i++) {
1334: for (j=0; j<vs->nc; j++) {
1335: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1336: Mat B = vs->m[i][j];
1337: if (!B) continue;
1338: MatGetSize(B,&M,&N);
1339: MatGetLocalSize(B,&m,&n);
1340: ISGetSize(vs->isglobal.row[i],&Mi);
1341: ISGetSize(vs->isglobal.col[j],&Ni);
1342: ISGetLocalSize(vs->isglobal.row[i],&mi);
1343: ISGetLocalSize(vs->isglobal.col[j],&ni);
1344: if (M != Mi || N != Ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Global sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",M,N,i,j,Mi,Ni);
1345: if (m != mi || n != ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Local sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",m,n,i,j,mi,ni);
1346: }
1347: }
1348: #endif
1350: /* Set A->assembled if all non-null blocks are currently assembled */
1351: for (i=0; i<vs->nr; i++) {
1352: for (j=0; j<vs->nc; j++) {
1353: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1354: }
1355: }
1356: A->assembled = PETSC_TRUE;
1357: return(0);
1358: }
1362: /*@C
1363: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1365: Collective on Mat
1367: Input Parameter:
1368: + comm - Communicator for the new Mat
1369: . nr - number of nested row blocks
1370: . is_row - index sets for each nested row block, or NULL to make contiguous
1371: . nc - number of nested column blocks
1372: . is_col - index sets for each nested column block, or NULL to make contiguous
1373: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1375: Output Parameter:
1376: . B - new matrix
1378: Level: advanced
1380: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1381: @*/
1382: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1383: {
1384: Mat A;
1388: *B = 0;
1389: MatCreate(comm,&A);
1390: MatSetType(A,MATNEST);
1391: MatSetUp(A);
1392: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1393: *B = A;
1394: return(0);
1395: }
1399: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1400: {
1402: Mat_Nest *nest = (Mat_Nest*)A->data;
1403: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1404: PetscInt cstart,cend;
1405: Mat C;
1408: MatGetSize(A,&M,&N);
1409: MatGetLocalSize(A,&m,&n);
1410: MatGetOwnershipRangeColumn(A,&cstart,&cend);
1411: switch (reuse) {
1412: case MAT_INITIAL_MATRIX:
1413: MatCreate(PetscObjectComm((PetscObject)A),&C);
1414: MatSetType(C,newtype);
1415: MatSetSizes(C,m,n,M,N);
1416: *newmat = C;
1417: break;
1418: case MAT_REUSE_MATRIX:
1419: C = *newmat;
1420: break;
1421: default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1422: }
1423: PetscMalloc1(2*m,&dnnz);
1424: onnz = dnnz + m;
1425: for (k=0; k<m; k++) {
1426: dnnz[k] = 0;
1427: onnz[k] = 0;
1428: }
1429: for (j=0; j<nest->nc; ++j) {
1430: IS bNis;
1431: PetscInt bN;
1432: const PetscInt *bNindices;
1433: /* Using global column indices and ISAllGather() is not scalable. */
1434: ISAllGather(nest->isglobal.col[j], &bNis);
1435: ISGetSize(bNis, &bN);
1436: ISGetIndices(bNis,&bNindices);
1437: for (i=0; i<nest->nr; ++i) {
1438: PetscSF bmsf;
1439: PetscSFNode *iremote;
1440: Mat B;
1441: PetscInt bm, *sub_dnnz,*sub_onnz, br;
1442: const PetscInt *bmindices;
1443: B = nest->m[i][j];
1444: if (!B) continue;
1445: ISGetLocalSize(nest->isglobal.row[i],&bm);
1446: ISGetIndices(nest->isglobal.row[i],&bmindices);
1447: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1448: PetscMalloc1(bm,&iremote);
1449: PetscMalloc1(bm,&sub_dnnz);
1450: PetscMalloc1(bm,&sub_onnz);
1451: for (k = 0; k < bm; ++k){
1452: sub_dnnz[k] = 0;
1453: sub_onnz[k] = 0;
1454: }
1455: /*
1456: Locate the owners for all of the locally-owned global row indices for this row block.
1457: These determine the roots of PetscSF used to communicate preallocation data to row owners.
1458: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1459: */
1460: MatGetOwnershipRange(B,&rstart,NULL);
1461: for (br = 0; br < bm; ++br) {
1462: PetscInt row = bmindices[br], rowowner = 0, brncols, col;
1463: const PetscInt *brcols;
1464: PetscInt rowrel = 0; /* row's relative index on its owner rank */
1465: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1466: /* how many roots */
1467: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
1468: /* get nonzero pattern */
1469: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
1470: for (k=0; k<brncols; k++) {
1471: col = bNindices[brcols[k]];
1472: if(col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]){
1473: sub_dnnz[br]++;
1474: }else{
1475: sub_onnz[br]++;
1476: }
1477: }
1478: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
1479: }
1480: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1481: /* bsf will have to take care of disposing of bedges. */
1482: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1483: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1484: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1485: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1486: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1487: PetscFree(sub_dnnz);
1488: PetscFree(sub_onnz);
1489: PetscSFDestroy(&bmsf);
1490: }
1491: ISRestoreIndices(bNis,&bNindices);
1492: ISDestroy(&bNis);
1493: }
1494: MatSeqAIJSetPreallocation(C,0,dnnz);
1495: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1496: PetscFree(dnnz);
1498: /* Fill by row */
1499: for (j=0; j<nest->nc; ++j) {
1500: /* Using global column indices and ISAllGather() is not scalable. */
1501: IS bNis;
1502: PetscInt bN;
1503: const PetscInt *bNindices;
1504: ISAllGather(nest->isglobal.col[j], &bNis);
1505: ISGetSize(bNis,&bN);
1506: ISGetIndices(bNis,&bNindices);
1507: for (i=0; i<nest->nr; ++i) {
1508: Mat B;
1509: PetscInt bm, br;
1510: const PetscInt *bmindices;
1511: B = nest->m[i][j];
1512: if (!B) continue;
1513: ISGetLocalSize(nest->isglobal.row[i],&bm);
1514: ISGetIndices(nest->isglobal.row[i],&bmindices);
1515: MatGetOwnershipRange(B,&rstart,NULL);
1516: for (br = 0; br < bm; ++br) {
1517: PetscInt row = bmindices[br], brncols, *cols;
1518: const PetscInt *brcols;
1519: const PetscScalar *brcoldata;
1520: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1521: PetscMalloc1(brncols,&cols);
1522: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1523: /*
1524: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1525: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1526: */
1527: MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1528: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1529: PetscFree(cols);
1530: }
1531: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1532: }
1533: ISRestoreIndices(bNis,&bNindices);
1534: ISDestroy(&bNis);
1535: }
1536: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1537: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1538: return(0);
1539: }
1541: /*MC
1542: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
1544: Level: intermediate
1546: Notes:
1547: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
1548: It allows the use of symmetric and block formats for parts of multi-physics simulations.
1549: It is usually used with DMComposite and DMCreateMatrix()
1551: .seealso: MatCreate(), MatType, MatCreateNest()
1552: M*/
1555: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1556: {
1557: Mat_Nest *s;
1561: PetscNewLog(A,&s);
1562: A->data = (void*)s;
1564: s->nr = -1;
1565: s->nc = -1;
1566: s->m = NULL;
1567: s->splitassembly = PETSC_FALSE;
1569: PetscMemzero(A->ops,sizeof(*A->ops));
1571: A->ops->mult = MatMult_Nest;
1572: A->ops->multadd = MatMultAdd_Nest;
1573: A->ops->multtranspose = MatMultTranspose_Nest;
1574: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
1575: A->ops->assemblybegin = MatAssemblyBegin_Nest;
1576: A->ops->assemblyend = MatAssemblyEnd_Nest;
1577: A->ops->zeroentries = MatZeroEntries_Nest;
1578: A->ops->copy = MatCopy_Nest;
1579: A->ops->duplicate = MatDuplicate_Nest;
1580: A->ops->getsubmatrix = MatGetSubMatrix_Nest;
1581: A->ops->destroy = MatDestroy_Nest;
1582: A->ops->view = MatView_Nest;
1583: A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1584: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
1585: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1586: A->ops->getdiagonal = MatGetDiagonal_Nest;
1587: A->ops->diagonalscale = MatDiagonalScale_Nest;
1588: A->ops->scale = MatScale_Nest;
1589: A->ops->shift = MatShift_Nest;
1591: A->spptr = 0;
1592: A->assembled = PETSC_FALSE;
1594: /* expose Nest api's */
1595: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
1596: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
1597: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
1598: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
1599: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
1600: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
1601: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
1602: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
1603: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);
1605: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1606: return(0);
1607: }