Actual source code: matnest.c
petsc-3.4.5 2014-06-29
2: #include ../src/mat/impls/nest/matnestimpl.h
3: #include <petscsf.h>
5: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
6: static PetscErrorCode MatGetVecs_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: MPI_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: *B = vs->m[row][col];
374: }
375: return(0);
376: }
380: static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
381: {
383: Mat_Nest *vs = (Mat_Nest*)A->data;
384: Mat sub;
387: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
388: switch (reuse) {
389: case MAT_INITIAL_MATRIX:
390: if (sub) { PetscObjectReference((PetscObject)sub); }
391: *B = sub;
392: break;
393: case MAT_REUSE_MATRIX:
394: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
395: break;
396: case MAT_IGNORE_MATRIX: /* Nothing to do */
397: break;
398: }
399: return(0);
400: }
404: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
405: {
407: Mat_Nest *vs = (Mat_Nest*)A->data;
408: Mat sub;
411: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
412: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
413: if (sub) {PetscObjectReference((PetscObject)sub);}
414: *B = sub;
415: return(0);
416: }
420: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
421: {
423: Mat_Nest *vs = (Mat_Nest*)A->data;
424: Mat sub;
427: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
428: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
429: if (sub) {
430: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
431: MatDestroy(B);
432: }
433: return(0);
434: }
438: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
439: {
440: Mat_Nest *bA = (Mat_Nest*)A->data;
441: PetscInt i;
445: for (i=0; i<bA->nr; i++) {
446: Vec bv;
447: VecGetSubVector(v,bA->isglobal.row[i],&bv);
448: if (bA->m[i][i]) {
449: MatGetDiagonal(bA->m[i][i],bv);
450: } else {
451: VecSet(bv,1.0);
452: }
453: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
454: }
455: return(0);
456: }
460: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
461: {
462: Mat_Nest *bA = (Mat_Nest*)A->data;
463: Vec bl,*br;
464: PetscInt i,j;
468: PetscMalloc(bA->nc*sizeof(Vec),&br);
469: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
470: for (i=0; i<bA->nr; i++) {
471: VecGetSubVector(l,bA->isglobal.row[i],&bl);
472: for (j=0; j<bA->nc; j++) {
473: if (bA->m[i][j]) {
474: MatDiagonalScale(bA->m[i][j],bl,br[j]);
475: }
476: }
477: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
478: }
479: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
480: PetscFree(br);
481: return(0);
482: }
486: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
487: {
488: Mat_Nest *bA = (Mat_Nest*)A->data;
489: PetscInt i,j;
493: for (i=0; i<bA->nr; i++) {
494: for (j=0; j<bA->nc; j++) {
495: if (bA->m[i][j]) {
496: MatScale(bA->m[i][j],a);
497: }
498: }
499: }
500: return(0);
501: }
505: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
506: {
507: Mat_Nest *bA = (Mat_Nest*)A->data;
508: PetscInt i;
512: for (i=0; i<bA->nr; i++) {
513: 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);
514: MatShift(bA->m[i][i],a);
515: }
516: return(0);
517: }
521: static PetscErrorCode MatGetVecs_Nest(Mat A,Vec *right,Vec *left)
522: {
523: Mat_Nest *bA = (Mat_Nest*)A->data;
524: Vec *L,*R;
525: MPI_Comm comm;
526: PetscInt i,j;
530: PetscObjectGetComm((PetscObject)A,&comm);
531: if (right) {
532: /* allocate R */
533: PetscMalloc(sizeof(Vec) * bA->nc, &R);
534: /* Create the right vectors */
535: for (j=0; j<bA->nc; j++) {
536: for (i=0; i<bA->nr; i++) {
537: if (bA->m[i][j]) {
538: MatGetVecs(bA->m[i][j],&R[j],NULL);
539: break;
540: }
541: }
542: if (i==bA->nr) {
543: /* have an empty column */
544: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
545: }
546: }
547: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
548: /* hand back control to the nest vector */
549: for (j=0; j<bA->nc; j++) {
550: VecDestroy(&R[j]);
551: }
552: PetscFree(R);
553: }
555: if (left) {
556: /* allocate L */
557: PetscMalloc(sizeof(Vec) * bA->nr, &L);
558: /* Create the left vectors */
559: for (i=0; i<bA->nr; i++) {
560: for (j=0; j<bA->nc; j++) {
561: if (bA->m[i][j]) {
562: MatGetVecs(bA->m[i][j],NULL,&L[i]);
563: break;
564: }
565: }
566: if (j==bA->nc) {
567: /* have an empty row */
568: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
569: }
570: }
572: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
573: for (i=0; i<bA->nr; i++) {
574: VecDestroy(&L[i]);
575: }
577: PetscFree(L);
578: }
579: return(0);
580: }
584: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
585: {
586: Mat_Nest *bA = (Mat_Nest*)A->data;
587: PetscBool isascii;
588: PetscInt i,j;
592: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
593: if (isascii) {
595: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
596: PetscViewerASCIIPushTab(viewer); /* push0 */
597: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc);
599: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
600: for (i=0; i<bA->nr; i++) {
601: for (j=0; j<bA->nc; j++) {
602: MatType type;
603: char name[256] = "",prefix[256] = "";
604: PetscInt NR,NC;
605: PetscBool isNest = PETSC_FALSE;
607: if (!bA->m[i][j]) {
608: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
609: continue;
610: }
611: MatGetSize(bA->m[i][j],&NR,&NC);
612: MatGetType(bA->m[i][j], &type);
613: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
614: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
615: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
617: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
619: if (isNest) {
620: PetscViewerASCIIPushTab(viewer); /* push1 */
621: MatView(bA->m[i][j],viewer);
622: PetscViewerASCIIPopTab(viewer); /* pop1 */
623: }
624: }
625: }
626: PetscViewerASCIIPopTab(viewer); /* pop0 */
627: }
628: return(0);
629: }
633: static PetscErrorCode MatZeroEntries_Nest(Mat A)
634: {
635: Mat_Nest *bA = (Mat_Nest*)A->data;
636: PetscInt i,j;
640: for (i=0; i<bA->nr; i++) {
641: for (j=0; j<bA->nc; j++) {
642: if (!bA->m[i][j]) continue;
643: MatZeroEntries(bA->m[i][j]);
644: }
645: }
646: return(0);
647: }
651: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
652: {
653: Mat_Nest *bA = (Mat_Nest*)A->data;
654: Mat *b;
655: PetscInt i,j,nr = bA->nr,nc = bA->nc;
659: PetscMalloc(nr*nc*sizeof(Mat),&b);
660: for (i=0; i<nr; i++) {
661: for (j=0; j<nc; j++) {
662: if (bA->m[i][j]) {
663: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
664: } else {
665: b[i*nc+j] = NULL;
666: }
667: }
668: }
669: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
670: /* Give the new MatNest exclusive ownership */
671: for (i=0; i<nr*nc; i++) {
672: MatDestroy(&b[i]);
673: }
674: PetscFree(b);
676: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
677: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
678: return(0);
679: }
681: /* nest api */
684: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
685: {
686: Mat_Nest *bA = (Mat_Nest*)A->data;
689: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
690: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
691: *mat = bA->m[idxm][jdxm];
692: return(0);
693: }
697: /*@
698: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
700: Not collective
702: Input Parameters:
703: + A - nest matrix
704: . idxm - index of the matrix within the nest matrix
705: - jdxm - index of the matrix within the nest matrix
707: Output Parameter:
708: . sub - matrix at index idxm,jdxm within the nest matrix
710: Level: developer
712: .seealso: MatNestGetSize(), MatNestGetSubMats()
713: @*/
714: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
715: {
719: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
720: return(0);
721: }
725: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
726: {
727: Mat_Nest *bA = (Mat_Nest*)A->data;
728: PetscInt m,n,M,N,mi,ni,Mi,Ni;
732: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
733: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
734: MatGetLocalSize(mat,&m,&n);
735: MatGetSize(mat,&M,&N);
736: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
737: ISGetSize(bA->isglobal.row[idxm],&Mi);
738: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
739: ISGetSize(bA->isglobal.col[jdxm],&Ni);
740: 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);
741: 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);
743: PetscObjectReference((PetscObject)mat);
744: MatDestroy(&bA->m[idxm][jdxm]);
745: bA->m[idxm][jdxm] = mat;
746: return(0);
747: }
751: /*@
752: MatNestSetSubMat - Set a single submatrix in the nest matrix.
754: Logically collective on the submatrix communicator
756: Input Parameters:
757: + A - nest matrix
758: . idxm - index of the matrix within the nest matrix
759: . jdxm - index of the matrix within the nest matrix
760: - sub - matrix at index idxm,jdxm within the nest matrix
762: Notes:
763: The new submatrix must have the same size and communicator as that block of the nest.
765: This increments the reference count of the submatrix.
767: Level: developer
769: .seealso: MatNestSetSubMats(), MatNestGetSubMat()
770: @*/
771: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
772: {
776: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
777: return(0);
778: }
782: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
783: {
784: Mat_Nest *bA = (Mat_Nest*)A->data;
787: if (M) *M = bA->nr;
788: if (N) *N = bA->nc;
789: if (mat) *mat = bA->m;
790: return(0);
791: }
795: /*@C
796: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
798: Not collective
800: Input Parameters:
801: . A - nest matrix
803: Output Parameter:
804: + M - number of rows in the nest matrix
805: . N - number of cols in the nest matrix
806: - mat - 2d array of matrices
808: Notes:
810: The user should not free the array mat.
812: Level: developer
814: .seealso: MatNestGetSize(), MatNestGetSubMat()
815: @*/
816: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
817: {
821: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
822: return(0);
823: }
827: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
828: {
829: Mat_Nest *bA = (Mat_Nest*)A->data;
832: if (M) *M = bA->nr;
833: if (N) *N = bA->nc;
834: return(0);
835: }
839: /*@
840: MatNestGetSize - Returns the size of the nest matrix.
842: Not collective
844: Input Parameters:
845: . A - nest matrix
847: Output Parameter:
848: + M - number of rows in the nested mat
849: - N - number of cols in the nested mat
851: Notes:
853: Level: developer
855: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
856: @*/
857: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
858: {
862: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
863: return(0);
864: }
868: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
869: {
870: Mat_Nest *vs = (Mat_Nest*)A->data;
871: PetscInt i;
874: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
875: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
876: return(0);
877: }
881: /*@C
882: MatNestGetISs - Returns the index sets partitioning the row and column spaces
884: Not collective
886: Input Parameters:
887: . A - nest matrix
889: Output Parameter:
890: + rows - array of row index sets
891: - cols - array of column index sets
893: Level: advanced
895: Notes:
896: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
898: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
899: @*/
900: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
901: {
906: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
907: return(0);
908: }
912: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
913: {
914: Mat_Nest *vs = (Mat_Nest*)A->data;
915: PetscInt i;
918: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
919: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
920: return(0);
921: }
925: /*@C
926: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
928: Not collective
930: Input Parameters:
931: . A - nest matrix
933: Output Parameter:
934: + rows - array of row index sets (or NULL to ignore)
935: - cols - array of column index sets (or NULL to ignore)
937: Level: advanced
939: Notes:
940: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
942: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
943: @*/
944: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
945: {
950: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
951: return(0);
952: }
956: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
957: {
959: PetscBool flg;
962: PetscStrcmp(vtype,VECNEST,&flg);
963: /* In reality, this only distinguishes VECNEST and "other" */
964: if (flg) A->ops->getvecs = MatGetVecs_Nest;
965: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
966: return(0);
967: }
971: /*@C
972: MatNestSetVecType - Sets the type of Vec returned by MatGetVecs()
974: Not collective
976: Input Parameters:
977: + A - nest matrix
978: - vtype - type to use for creating vectors
980: Notes:
982: Level: developer
984: .seealso: MatGetVecs()
985: @*/
986: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
987: {
991: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
992: return(0);
993: }
997: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
998: {
999: Mat_Nest *s = (Mat_Nest*)A->data;
1000: PetscInt i,j,m,n,M,N;
1004: s->nr = nr;
1005: s->nc = nc;
1007: /* Create space for submatrices */
1008: PetscMalloc(sizeof(Mat*)*nr,&s->m);
1009: for (i=0; i<nr; i++) {
1010: PetscMalloc(sizeof(Mat)*nc,&s->m[i]);
1011: }
1012: for (i=0; i<nr; i++) {
1013: for (j=0; j<nc; j++) {
1014: s->m[i][j] = a[i*nc+j];
1015: if (a[i*nc+j]) {
1016: PetscObjectReference((PetscObject)a[i*nc+j]);
1017: }
1018: }
1019: }
1021: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1023: PetscMalloc(sizeof(PetscInt)*nr,&s->row_len);
1024: PetscMalloc(sizeof(PetscInt)*nc,&s->col_len);
1025: for (i=0; i<nr; i++) s->row_len[i]=-1;
1026: for (j=0; j<nc; j++) s->col_len[j]=-1;
1028: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1030: PetscLayoutSetSize(A->rmap,M);
1031: PetscLayoutSetLocalSize(A->rmap,m);
1032: PetscLayoutSetSize(A->cmap,N);
1033: PetscLayoutSetLocalSize(A->cmap,n);
1035: PetscLayoutSetUp(A->rmap);
1036: PetscLayoutSetUp(A->cmap);
1038: PetscMalloc2(nr,Vec,&s->left,nc,Vec,&s->right);
1039: PetscMemzero(s->left,nr*sizeof(Vec));
1040: PetscMemzero(s->right,nc*sizeof(Vec));
1041: return(0);
1042: }
1046: /*@
1047: MatNestSetSubMats - Sets the nested submatrices
1049: Collective on Mat
1051: Input Parameter:
1052: + N - nested matrix
1053: . nr - number of nested row blocks
1054: . is_row - index sets for each nested row block, or NULL to make contiguous
1055: . nc - number of nested column blocks
1056: . is_col - index sets for each nested column block, or NULL to make contiguous
1057: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1059: Level: advanced
1061: .seealso: MatCreateNest(), MATNEST
1062: @*/
1063: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1064: {
1066: PetscInt i;
1070: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1071: if (nr && is_row) {
1074: }
1075: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1076: if (nc && is_col) {
1079: }
1081: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1082: return(0);
1083: }
1087: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog,ISLocalToGlobalMapping *ltogb)
1088: {
1090: PetscBool flg;
1091: PetscInt i,j,m,mi,*ix;
1094: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1095: if (islocal[i]) {
1096: ISGetSize(islocal[i],&mi);
1097: flg = PETSC_TRUE; /* We found a non-trivial entry */
1098: } else {
1099: ISGetSize(isglobal[i],&mi);
1100: }
1101: m += mi;
1102: }
1103: if (flg) {
1104: PetscMalloc(m*sizeof(*ix),&ix);
1105: for (i=0,n=0; i<n; i++) {
1106: ISLocalToGlobalMapping smap = NULL;
1107: VecScatter scat;
1108: IS isreq;
1109: Vec lvec,gvec;
1110: union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x;
1111: Mat sub;
1113: if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers");
1114: if (colflg) {
1115: MatNestFindNonzeroSubMatRow(A,i,&sub);
1116: } else {
1117: MatNestFindNonzeroSubMatCol(A,i,&sub);
1118: }
1119: if (sub) {MatGetLocalToGlobalMapping(sub,&smap,NULL);}
1120: if (islocal[i]) {
1121: ISGetSize(islocal[i],&mi);
1122: } else {
1123: ISGetSize(isglobal[i],&mi);
1124: }
1125: for (j=0; j<mi; j++) ix[m+j] = j;
1126: if (smap) {ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);}
1127: /*
1128: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1129: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1130: The approach here is ugly because it uses VecScatter to move indices.
1131: */
1132: VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);
1133: VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);
1134: ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);
1135: VecScatterCreate(gvec,isreq,lvec,NULL,&scat);
1136: VecGetArray(gvec,(PetscScalar**)&x);
1137: for (j=0; j<mi; j++) x[j].integer = ix[m+j];
1138: VecRestoreArray(gvec,(PetscScalar**)&x);
1139: VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1140: VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1141: VecGetArray(lvec,(PetscScalar**)&x);
1142: for (j=0; j<mi; j++) ix[m+j] = x[j].integer;
1143: VecRestoreArray(lvec,(PetscScalar**)&x);
1144: VecDestroy(&lvec);
1145: VecDestroy(&gvec);
1146: ISDestroy(&isreq);
1147: VecScatterDestroy(&scat);
1148: m += mi;
1149: }
1150: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),m,ix,PETSC_OWN_POINTER,ltog);
1151: *ltogb = NULL;
1152: } else {
1153: *ltog = NULL;
1154: *ltogb = NULL;
1155: }
1156: return(0);
1157: }
1160: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1161: /*
1162: nprocessors = NP
1163: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1164: proc 0: => (g_0,h_0,)
1165: proc 1: => (g_1,h_1,)
1166: ...
1167: proc nprocs-1: => (g_NP-1,h_NP-1,)
1169: proc 0: proc 1: proc nprocs-1:
1170: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1172: proc 0:
1173: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1174: proc 1:
1175: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1177: proc NP-1:
1178: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1179: */
1182: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1183: {
1184: Mat_Nest *vs = (Mat_Nest*)A->data;
1185: PetscInt i,j,offset,n,nsum,bs;
1187: Mat sub = NULL;
1190: PetscMalloc(sizeof(IS)*nr,&vs->isglobal.row);
1191: PetscMalloc(sizeof(IS)*nc,&vs->isglobal.col);
1192: if (is_row) { /* valid IS is passed in */
1193: /* refs on is[] are incremeneted */
1194: for (i=0; i<vs->nr; i++) {
1195: PetscObjectReference((PetscObject)is_row[i]);
1197: vs->isglobal.row[i] = is_row[i];
1198: }
1199: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1200: nsum = 0;
1201: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1202: MatNestFindNonzeroSubMatRow(A,i,&sub);
1203: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1204: MatGetLocalSize(sub,&n,NULL);
1205: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1206: nsum += n;
1207: }
1208: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1209: offset -= nsum;
1210: for (i=0; i<vs->nr; i++) {
1211: MatNestFindNonzeroSubMatRow(A,i,&sub);
1212: MatGetLocalSize(sub,&n,NULL);
1213: MatGetBlockSize(sub,&bs);
1214: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1215: ISSetBlockSize(vs->isglobal.row[i],bs);
1216: offset += n;
1217: }
1218: }
1220: if (is_col) { /* valid IS is passed in */
1221: /* refs on is[] are incremeneted */
1222: for (j=0; j<vs->nc; j++) {
1223: PetscObjectReference((PetscObject)is_col[j]);
1225: vs->isglobal.col[j] = is_col[j];
1226: }
1227: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1228: offset = A->cmap->rstart;
1229: nsum = 0;
1230: for (j=0; j<vs->nc; j++) {
1231: MatNestFindNonzeroSubMatCol(A,j,&sub);
1232: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1233: MatGetLocalSize(sub,NULL,&n);
1234: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1235: nsum += n;
1236: }
1237: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1238: offset -= nsum;
1239: for (j=0; j<vs->nc; j++) {
1240: MatNestFindNonzeroSubMatCol(A,j,&sub);
1241: MatGetLocalSize(sub,NULL,&n);
1242: MatGetBlockSize(sub,&bs);
1243: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1244: ISSetBlockSize(vs->isglobal.col[j],bs);
1245: offset += n;
1246: }
1247: }
1249: /* Set up the local ISs */
1250: PetscMalloc(vs->nr*sizeof(IS),&vs->islocal.row);
1251: PetscMalloc(vs->nc*sizeof(IS),&vs->islocal.col);
1252: for (i=0,offset=0; i<vs->nr; i++) {
1253: IS isloc;
1254: ISLocalToGlobalMapping rmap = NULL;
1255: PetscInt nlocal,bs;
1256: MatNestFindNonzeroSubMatRow(A,i,&sub);
1257: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1258: if (rmap) {
1259: MatGetBlockSize(sub,&bs);
1260: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1261: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1262: ISSetBlockSize(isloc,bs);
1263: } else {
1264: nlocal = 0;
1265: isloc = NULL;
1266: }
1267: vs->islocal.row[i] = isloc;
1268: offset += nlocal;
1269: }
1270: for (i=0,offset=0; i<vs->nc; i++) {
1271: IS isloc;
1272: ISLocalToGlobalMapping cmap = NULL;
1273: PetscInt nlocal,bs;
1274: MatNestFindNonzeroSubMatCol(A,i,&sub);
1275: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1276: if (cmap) {
1277: MatGetBlockSize(sub,&bs);
1278: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1279: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1280: ISSetBlockSize(isloc,bs);
1281: } else {
1282: nlocal = 0;
1283: isloc = NULL;
1284: }
1285: vs->islocal.col[i] = isloc;
1286: offset += nlocal;
1287: }
1289: /* Set up the aggregate ISLocalToGlobalMapping */
1290: {
1291: ISLocalToGlobalMapping rmap,rmapb,cmap,cmapb;
1292: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap,&rmapb);
1293: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap,&cmapb);
1294: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1295: if (rmapb && cmapb) {MatSetLocalToGlobalMappingBlock(A,rmapb,cmapb);}
1296: ISLocalToGlobalMappingDestroy(&rmap);
1297: ISLocalToGlobalMappingDestroy(&rmapb);
1298: ISLocalToGlobalMappingDestroy(&cmap);
1299: ISLocalToGlobalMappingDestroy(&cmapb);
1300: }
1302: #if defined(PETSC_USE_DEBUG)
1303: for (i=0; i<vs->nr; i++) {
1304: for (j=0; j<vs->nc; j++) {
1305: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1306: Mat B = vs->m[i][j];
1307: if (!B) continue;
1308: MatGetSize(B,&M,&N);
1309: MatGetLocalSize(B,&m,&n);
1310: ISGetSize(vs->isglobal.row[i],&Mi);
1311: ISGetSize(vs->isglobal.col[j],&Ni);
1312: ISGetLocalSize(vs->isglobal.row[i],&mi);
1313: ISGetLocalSize(vs->isglobal.col[j],&ni);
1314: 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);
1315: 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);
1316: }
1317: }
1318: #endif
1320: /* Set A->assembled if all non-null blocks are currently assembled */
1321: for (i=0; i<vs->nr; i++) {
1322: for (j=0; j<vs->nc; j++) {
1323: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1324: }
1325: }
1326: A->assembled = PETSC_TRUE;
1327: return(0);
1328: }
1332: /*@C
1333: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1335: Collective on Mat
1337: Input Parameter:
1338: + comm - Communicator for the new Mat
1339: . nr - number of nested row blocks
1340: . is_row - index sets for each nested row block, or NULL to make contiguous
1341: . nc - number of nested column blocks
1342: . is_col - index sets for each nested column block, or NULL to make contiguous
1343: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1345: Output Parameter:
1346: . B - new matrix
1348: Level: advanced
1350: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1351: @*/
1352: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1353: {
1354: Mat A;
1358: *B = 0;
1359: MatCreate(comm,&A);
1360: MatSetType(A,MATNEST);
1361: MatSetUp(A);
1362: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1363: *B = A;
1364: return(0);
1365: }
1369: PETSC_EXTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1370: {
1372: Mat_Nest *nest = (Mat_Nest*)A->data;
1373: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz;
1374: Mat C;
1377: MatGetSize(A,&M,&N);
1378: MatGetLocalSize(A,&m,&n);
1379: switch (reuse) {
1380: case MAT_INITIAL_MATRIX:
1381: MatCreate(PetscObjectComm((PetscObject)A),&C);
1382: MatSetType(C,newtype);
1383: MatSetSizes(C,m,n,M,N);
1384: *newmat = C;
1385: break;
1386: case MAT_REUSE_MATRIX:
1387: C = *newmat;
1388: break;
1389: default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1390: }
1392: /* Preallocation */
1393: PetscMalloc(2*m*sizeof(PetscInt),&dnnz);
1394: onnz = dnnz + m;
1395: for (k=0; k<m; k++) {
1396: dnnz[k] = 0;
1397: onnz[k] = 0;
1398: }
1399: for (j=0; j<nest->nc; ++j) {
1400: IS bNis;
1401: PetscInt bN;
1402: const PetscInt *bNindices;
1403: /* Using global column indices and ISAllGather() is not scalable. */
1404: ISAllGather(nest->isglobal.col[j], &bNis);
1405: ISGetSize(bNis, &bN);
1406: ISGetIndices(bNis,&bNindices);
1407: for (i=0; i<nest->nr; ++i) {
1408: PetscSF bmsf;
1409: PetscSFNode *bmedges;
1410: Mat B;
1411: PetscInt bm, *bmdnnz, br;
1412: const PetscInt *bmindices;
1413: B = nest->m[i][j];
1414: if (!B) continue;
1415: ISGetLocalSize(nest->isglobal.row[i],&bm);
1416: ISGetIndices(nest->isglobal.row[i],&bmindices);
1417: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1418: PetscMalloc(2*bm*sizeof(PetscSFNode),&bmedges);
1419: PetscMalloc(2*bm*sizeof(PetscInt),&bmdnnz);
1420: for (k = 0; k < 2*bm; ++k) bmdnnz[k] = 0;
1421: /*
1422: Locate the owners for all of the locally-owned global row indices for this row block.
1423: These determine the roots of PetscSF used to communicate preallocation data to row owners.
1424: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1425: */
1426: for (br = 0; br < bm; ++br) {
1427: PetscInt row = bmindices[br], rowowner = 0, brncols, col, colowner = 0;
1428: const PetscInt *brcols;
1429: PetscInt rowrel = 0; /* row's relative index on its owner rank */
1430: PetscInt rowownerm; /* local row size on row's owning rank. */
1432: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1433: rowownerm = A->rmap->range[rowowner+1]-A->rmap->range[rowowner];
1435: bmedges[br].rank = rowowner; bmedges[br].index = rowrel; /* edge from bmdnnz to dnnz */
1436: bmedges[br].rank = rowowner; bmedges[br].index = rowrel+rowownerm; /* edge from bmonnz to onnz */
1437: /* Now actually compute the data -- bmdnnz and bmonnz by looking at the global columns in the br row of this block. */
1438: /* Note that this is not a pessimistic bound only because we assume the index sets embedding the blocks do not overlap. */
1439: MatGetRow(B,br,&brncols,&brcols,NULL);
1440: for (k=0; k<brncols; k++) {
1441: col = bNindices[brcols[k]];
1442: PetscLayoutFindOwnerIndex(A->cmap,col,&colowner,NULL);
1443: if (colowner == rowowner) bmdnnz[br]++;
1444: else onnz[br]++;
1445: }
1446: MatRestoreRow(B,br,&brncols,&brcols,NULL);
1447: }
1448: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1449: /* bsf will have to take care of disposing of bedges. */
1450: PetscSFSetGraph(bmsf,m,2*bm,NULL,PETSC_COPY_VALUES,bmedges,PETSC_OWN_POINTER);
1451: PetscSFReduceBegin(bmsf,MPIU_INT,bmdnnz,dnnz,MPIU_SUM);
1452: PetscSFReduceEnd(bmsf,MPIU_INT,bmdnnz,dnnz,MPIU_SUM);
1453: PetscFree(bmdnnz);
1454: PetscSFDestroy(&bmsf);
1455: }
1456: ISRestoreIndices(bNis,&bNindices);
1457: ISDestroy(&bNis);
1458: }
1459: MatSeqAIJSetPreallocation(C,0,dnnz);
1460: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1461: PetscFree(dnnz);
1463: /* Fill by row */
1464: for (j=0; j<nest->nc; ++j) {
1465: /* Using global column indices and ISAllGather() is not scalable. */
1466: IS bNis;
1467: PetscInt bN;
1468: const PetscInt *bNindices;
1469: ISAllGather(nest->isglobal.col[j], &bNis);
1470: ISGetSize(bNis,&bN);
1471: ISGetIndices(bNis,&bNindices);
1472: for (i=0; i<nest->nr; ++i) {
1473: Mat B;
1474: PetscInt bm, br;
1475: const PetscInt *bmindices;
1476: B = nest->m[i][j];
1477: if (!B) continue;
1478: ISGetLocalSize(nest->isglobal.row[i],&bm);
1479: ISGetIndices(nest->isglobal.row[i],&bmindices);
1480: for (br = 0; br < bm; ++br) {
1481: PetscInt row = bmindices[br], brncols, *cols;
1482: const PetscInt *brcols;
1483: const PetscScalar *brcoldata;
1484: MatGetRow(B,br,&brncols,&brcols,&brcoldata);
1485: PetscMalloc(brncols*sizeof(PetscInt),&cols);
1486: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1487: /*
1488: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1489: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1490: */
1491: MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1492: MatRestoreRow(B,br,&brncols,&brcols,&brcoldata);
1493: PetscFree(cols);
1494: }
1495: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1496: }
1497: ISRestoreIndices(bNis,&bNindices);
1498: ISDestroy(&bNis);
1499: }
1500: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1501: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1502: return(0);
1503: }
1505: /*MC
1506: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
1508: Level: intermediate
1510: Notes:
1511: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
1512: It allows the use of symmetric and block formats for parts of multi-physics simulations.
1513: It is usually used with DMComposite and DMCreateMatrix()
1515: .seealso: MatCreate(), MatType, MatCreateNest()
1516: M*/
1519: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1520: {
1521: Mat_Nest *s;
1525: PetscNewLog(A,Mat_Nest,&s);
1526: A->data = (void*)s;
1528: s->nr = -1;
1529: s->nc = -1;
1530: s->m = NULL;
1531: s->splitassembly = PETSC_FALSE;
1533: PetscMemzero(A->ops,sizeof(*A->ops));
1535: A->ops->mult = MatMult_Nest;
1536: A->ops->multadd = MatMultAdd_Nest;
1537: A->ops->multtranspose = MatMultTranspose_Nest;
1538: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
1539: A->ops->assemblybegin = MatAssemblyBegin_Nest;
1540: A->ops->assemblyend = MatAssemblyEnd_Nest;
1541: A->ops->zeroentries = MatZeroEntries_Nest;
1542: A->ops->duplicate = MatDuplicate_Nest;
1543: A->ops->getsubmatrix = MatGetSubMatrix_Nest;
1544: A->ops->destroy = MatDestroy_Nest;
1545: A->ops->view = MatView_Nest;
1546: A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1547: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
1548: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1549: A->ops->getdiagonal = MatGetDiagonal_Nest;
1550: A->ops->diagonalscale = MatDiagonalScale_Nest;
1551: A->ops->scale = MatScale_Nest;
1552: A->ops->shift = MatShift_Nest;
1554: A->spptr = 0;
1555: A->same_nonzero = PETSC_FALSE;
1556: A->assembled = PETSC_FALSE;
1558: /* expose Nest api's */
1559: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
1560: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
1561: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
1562: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
1563: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
1564: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
1565: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
1566: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
1568: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1569: return(0);
1570: }