Actual source code: matnest.c
petsc-3.3-p7 2013-05-11
2: #include matnestimpl.h
4: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
5: static PetscErrorCode MatGetVecs_Nest(Mat A,Vec *right,Vec *left);
7: /* private functions */
10: static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
11: {
12: Mat_Nest *bA = (Mat_Nest*)A->data;
13: PetscInt i,j;
17: *m = *n = *M = *N = 0;
18: for (i=0; i<bA->nr; i++) { /* rows */
19: PetscInt sm,sM;
20: ISGetLocalSize(bA->isglobal.row[i],&sm);
21: ISGetSize(bA->isglobal.row[i],&sM);
22: *m += sm;
23: *M += sM;
24: }
25: for (j=0; j<bA->nc; j++) { /* cols */
26: PetscInt sn,sN;
27: ISGetLocalSize(bA->isglobal.col[j],&sn);
28: ISGetSize(bA->isglobal.col[j],&sN);
29: *n += sn;
30: *N += sN;
31: }
32: return(0);
33: }
35: /* operations */
38: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
39: {
40: Mat_Nest *bA = (Mat_Nest*)A->data;
41: Vec *bx = bA->right,*by = bA->left;
42: PetscInt i,j,nr = bA->nr,nc = bA->nc;
46: for (i=0; i<nr; i++) {VecGetSubVector(y,bA->isglobal.row[i],&by[i]);}
47: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
48: for (i=0; i<nr; i++) {
49: VecZeroEntries(by[i]);
50: for (j=0; j<nc; j++) {
51: if (!bA->m[i][j]) continue;
52: /* y[i] <- y[i] + A[i][j] * x[j] */
53: MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);
54: }
55: }
56: for (i=0; i<nr; i++) {VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);}
57: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
58: return(0);
59: }
63: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
64: {
65: Mat_Nest *bA = (Mat_Nest*)A->data;
66: Vec *bx = bA->right,*bz = bA->left;
67: PetscInt i,j,nr = bA->nr,nc = bA->nc;
71: for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
72: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
73: for (i=0; i<nr; i++) {
74: if (y != z) {
75: Vec by;
76: VecGetSubVector(y,bA->isglobal.row[i],&by);
77: VecCopy(by,bz[i]);
78: VecRestoreSubVector(y,bA->isglobal.row[i],&by);
79: }
80: for (j=0; j<nc; j++) {
81: if (!bA->m[i][j]) continue;
82: /* y[i] <- y[i] + A[i][j] * x[j] */
83: MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
84: }
85: }
86: for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
87: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
88: return(0);
89: }
93: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
94: {
95: Mat_Nest *bA = (Mat_Nest*)A->data;
96: Vec *bx = bA->left,*by = bA->right;
97: PetscInt i,j,nr = bA->nr,nc = bA->nc;
101: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
102: for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
103: for (j=0; j<nc; j++) {
104: VecZeroEntries(by[j]);
105: for (i=0; i<nr; i++) {
106: if (!bA->m[i][j]) continue;
107: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
108: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
109: }
110: }
111: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
112: for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
113: return(0);
114: }
118: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
119: {
120: Mat_Nest *bA = (Mat_Nest*)A->data;
121: Vec *bx = bA->left,*bz = bA->right;
122: PetscInt i,j,nr = bA->nr,nc = bA->nc;
126: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
127: for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
128: for (j=0; j<nc; j++) {
129: if (y != z) {
130: Vec by;
131: VecGetSubVector(y,bA->isglobal.col[j],&by);
132: VecCopy(by,bz[j]);
133: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
134: }
135: for (i=0; i<nr; i++) {
136: if (!bA->m[i][j]) continue;
137: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
138: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
139: }
140: }
141: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
142: for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
143: return(0);
144: }
148: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
149: {
151: IS *lst = *list;
152: PetscInt i;
155: if (!lst) return(0);
156: for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
157: PetscFree(lst);
158: *list = PETSC_NULL;
159: return(0);
160: }
164: static PetscErrorCode MatDestroy_Nest(Mat A)
165: {
166: Mat_Nest *vs = (Mat_Nest*)A->data;
167: PetscInt i,j;
171: /* release the matrices and the place holders */
172: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
173: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
174: MatNestDestroyISList(vs->nr,&vs->islocal.row);
175: MatNestDestroyISList(vs->nc,&vs->islocal.col);
177: PetscFree(vs->row_len);
178: PetscFree(vs->col_len);
180: PetscFree2(vs->left,vs->right);
182: /* release the matrices and the place holders */
183: if (vs->m) {
184: for (i=0; i<vs->nr; i++) {
185: for (j=0; j<vs->nc; j++) {
186: MatDestroy(&vs->m[i][j]);
187: }
188: PetscFree( vs->m[i] );
189: }
190: PetscFree(vs->m);
191: }
192: PetscFree(A->data);
194: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMat_C", "",0);
195: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestSetSubMat_C", "",0);
196: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMats_C", "",0);
197: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSize_C", "",0);
198: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetISs_C", "",0);
199: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetLocalISs_C", "",0);
200: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestSetVecType_C", "",0);
201: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestSetSubMats_C", "",0);
202: return(0);
203: }
207: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
208: {
209: Mat_Nest *vs = (Mat_Nest*)A->data;
210: PetscInt i,j;
214: for (i=0; i<vs->nr; i++) {
215: for (j=0; j<vs->nc; j++) {
216: if (vs->m[i][j]) {
217: MatAssemblyBegin(vs->m[i][j],type);
218: if (!vs->splitassembly) {
219: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
220: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
221: * already performing an assembly, but the result would by more complicated and appears to offer less
222: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
223: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
224: */
225: MatAssemblyEnd(vs->m[i][j],type);
226: }
227: }
228: }
229: }
230: return(0);
231: }
235: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
236: {
237: Mat_Nest *vs = (Mat_Nest*)A->data;
238: PetscInt i,j;
242: for (i=0; i<vs->nr; i++) {
243: for (j=0; j<vs->nc; j++) {
244: if (vs->m[i][j]) {
245: if (vs->splitassembly) {
246: MatAssemblyEnd(vs->m[i][j],type);
247: }
248: }
249: }
250: }
251: return(0);
252: }
256: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
257: {
259: Mat_Nest *vs = (Mat_Nest*)A->data;
260: PetscInt j;
261: Mat sub;
264: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)PETSC_NULL; /* Prefer to find on the diagonal */
265: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
266: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
267: *B = sub;
268: return(0);
269: }
273: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
274: {
276: Mat_Nest *vs = (Mat_Nest*)A->data;
277: PetscInt i;
278: Mat sub;
281: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)PETSC_NULL; /* Prefer to find on the diagonal */
282: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
283: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
284: *B = sub;
285: return(0);
286: }
290: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
291: {
293: PetscInt i;
294: PetscBool flg;
300: *found = -1;
301: for (i=0; i<n; i++) {
302: if (!list[i]) continue;
303: ISEqual(list[i],is,&flg);
304: if (flg) {
305: *found = i;
306: return(0);
307: }
308: }
309: SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_INCOMP,"Could not find index set");
310: return(0);
311: }
315: /* Get a block row as a new MatNest */
316: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
317: {
318: Mat_Nest *vs = (Mat_Nest*)A->data;
319: char keyname[256];
323: *B = PETSC_NULL;
324: PetscSNPrintf(keyname,sizeof keyname,"NestRow_%D",row);
325: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
326: if (*B) return(0);
328: MatCreateNest(((PetscObject)A)->comm,1,PETSC_NULL,vs->nc,vs->isglobal.col,vs->m[row],B);
329: (*B)->assembled = A->assembled;
330: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
331: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
332: return(0);
333: }
337: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
338: {
339: Mat_Nest *vs = (Mat_Nest*)A->data;
341: PetscInt row,col;
342: PetscBool same,isFullCol,isFullColGlobal;
345: /* Check if full column space. This is a hack */
346: isFullCol = PETSC_FALSE;
347: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
348: if (same) {
349: PetscInt n,first,step,i,an,am,afirst,astep;
350: ISStrideGetInfo(iscol,&first,&step);
351: ISGetLocalSize(iscol,&n);
352: isFullCol = PETSC_TRUE;
353: for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
354: ISStrideGetInfo(is->col[i],&afirst,&astep);
355: ISGetLocalSize(is->col[i],&am);
356: if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
357: an += am;
358: }
359: if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
360: }
361: MPI_Allreduce(&isFullCol,&isFullColGlobal,1,MPI_INT,MPI_LAND,((PetscObject)iscol)->comm);
363: if (isFullColGlobal) {
364: PetscInt row;
365: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
366: MatNestGetRow(A,row,B);
367: } else {
368: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
369: MatNestFindIS(A,vs->nc,is->col,iscol,&col);
370: *B = vs->m[row][col];
371: }
372: return(0);
373: }
377: static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
378: {
380: Mat_Nest *vs = (Mat_Nest*)A->data;
381: Mat sub;
384: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
385: switch (reuse) {
386: case MAT_INITIAL_MATRIX:
387: if (sub) { PetscObjectReference((PetscObject)sub); }
388: *B = sub;
389: break;
390: case MAT_REUSE_MATRIX:
391: if (sub != *B) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
392: break;
393: case MAT_IGNORE_MATRIX: /* Nothing to do */
394: break;
395: }
396: return(0);
397: }
401: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
402: {
404: Mat_Nest *vs = (Mat_Nest*)A->data;
405: Mat sub;
408: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
409: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
410: if (sub) {PetscObjectReference((PetscObject)sub);}
411: *B = sub;
412: return(0);
413: }
417: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
418: {
420: Mat_Nest *vs = (Mat_Nest*)A->data;
421: Mat sub;
424: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
425: if (*B != sub) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
426: if (sub) {
427: if (((PetscObject)sub)->refct <= 1) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
428: MatDestroy(B);
429: }
430: return(0);
431: }
435: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
436: {
437: Mat_Nest *bA = (Mat_Nest*)A->data;
438: PetscInt i;
442: for (i=0; i<bA->nr; i++) {
443: Vec bv;
444: VecGetSubVector(v,bA->isglobal.row[i],&bv);
445: if (bA->m[i][i]) {
446: MatGetDiagonal(bA->m[i][i],bv);
447: } else {
448: VecSet(bv,1.0);
449: }
450: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
451: }
452: return(0);
453: }
457: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
458: {
459: Mat_Nest *bA = (Mat_Nest*)A->data;
460: Vec bl,*br;
461: PetscInt i,j;
465: PetscMalloc(bA->nc*sizeof(Vec),&br);
466: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
467: for (i=0; i<bA->nr; i++) {
468: VecGetSubVector(l,bA->isglobal.row[i],&bl);
469: for (j=0; j<bA->nc; j++) {
470: if (bA->m[i][j]) {
471: MatDiagonalScale(bA->m[i][j],bl,br[j]);
472: }
473: }
474: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
475: }
476: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
477: PetscFree(br);
478: return(0);
479: }
483: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
484: {
485: Mat_Nest *bA = (Mat_Nest*)A->data;
486: PetscInt i,j;
490: for (i=0; i<bA->nr; i++) {
491: for (j=0; j<bA->nc; j++) {
492: if (bA->m[i][j]) {
493: MatScale(bA->m[i][j],a);
494: }
495: }
496: }
497: return(0);
498: }
502: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
503: {
504: Mat_Nest *bA = (Mat_Nest*)A->data;
505: PetscInt i;
509: for (i=0; i<bA->nr; i++) {
510: if (!bA->m[i][i]) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_SUP,"No support for shifting an empty diagonal block, insert a matrix in block (%D,%D)",i,i);
511: MatShift(bA->m[i][i],a);
512: }
513: return(0);
514: }
518: static PetscErrorCode MatGetVecs_Nest(Mat A,Vec *right,Vec *left)
519: {
520: Mat_Nest *bA = (Mat_Nest*)A->data;
521: Vec *L,*R;
522: MPI_Comm comm;
523: PetscInt i,j;
527: comm = ((PetscObject)A)->comm;
528: if (right) {
529: /* allocate R */
530: PetscMalloc( sizeof(Vec) * bA->nc, &R );
531: /* Create the right vectors */
532: for (j=0; j<bA->nc; j++) {
533: for (i=0; i<bA->nr; i++) {
534: if (bA->m[i][j]) {
535: MatGetVecs(bA->m[i][j],&R[j],PETSC_NULL);
536: break;
537: }
538: }
539: if (i==bA->nr) {
540: /* have an empty column */
541: SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
542: }
543: }
544: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
545: /* hand back control to the nest vector */
546: for (j=0; j<bA->nc; j++) {
547: VecDestroy(&R[j]);
548: }
549: PetscFree(R);
550: }
552: if (left) {
553: /* allocate L */
554: PetscMalloc( sizeof(Vec) * bA->nr, &L );
555: /* Create the left vectors */
556: for (i=0; i<bA->nr; i++) {
557: for (j=0; j<bA->nc; j++) {
558: if (bA->m[i][j]) {
559: MatGetVecs(bA->m[i][j],PETSC_NULL,&L[i]);
560: break;
561: }
562: }
563: if (j==bA->nc) {
564: /* have an empty row */
565: SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
566: }
567: }
569: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
570: for (i=0; i<bA->nr; i++) {
571: VecDestroy(&L[i]);
572: }
574: PetscFree(L);
575: }
576: return(0);
577: }
581: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
582: {
583: Mat_Nest *bA = (Mat_Nest*)A->data;
584: PetscBool isascii;
585: PetscInt i,j;
589: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
590: if (isascii) {
592: PetscViewerASCIIPrintf(viewer,"Matrix object: \n" );
593: PetscViewerASCIIPushTab( viewer ); /* push0 */
594: PetscViewerASCIIPrintf( viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc);
596: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n" );
597: for (i=0; i<bA->nr; i++) {
598: for (j=0; j<bA->nc; j++) {
599: const MatType type;
600: char name[256] = "",prefix[256] = "";
601: PetscInt NR,NC;
602: PetscBool isNest = PETSC_FALSE;
604: if (!bA->m[i][j]) {
605: PetscViewerASCIIPrintf( viewer, "(%D,%D) : PETSC_NULL \n",i,j);
606: continue;
607: }
608: MatGetSize(bA->m[i][j],&NR,&NC);
609: MatGetType( bA->m[i][j], &type );
610: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof name,"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
611: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof prefix,"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
612: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
614: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
616: if (isNest) {
617: PetscViewerASCIIPushTab(viewer); /* push1 */
618: MatView(bA->m[i][j],viewer);
619: PetscViewerASCIIPopTab(viewer); /* pop1 */
620: }
621: }
622: }
623: PetscViewerASCIIPopTab(viewer); /* pop0 */
624: }
625: return(0);
626: }
630: static PetscErrorCode MatZeroEntries_Nest(Mat A)
631: {
632: Mat_Nest *bA = (Mat_Nest*)A->data;
633: PetscInt i,j;
637: for (i=0; i<bA->nr; i++) {
638: for (j=0; j<bA->nc; j++) {
639: if (!bA->m[i][j]) continue;
640: MatZeroEntries(bA->m[i][j]);
641: }
642: }
643: return(0);
644: }
648: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
649: {
650: Mat_Nest *bA = (Mat_Nest*)A->data;
651: Mat *b;
652: PetscInt i,j,nr = bA->nr,nc = bA->nc;
656: PetscMalloc(nr*nc*sizeof(Mat),&b);
657: for (i=0; i<nr; i++) {
658: for (j=0; j<nc; j++) {
659: if (bA->m[i][j]) {
660: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
661: } else {
662: b[i*nc+j] = PETSC_NULL;
663: }
664: }
665: }
666: MatCreateNest(((PetscObject)A)->comm,nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
667: /* Give the new MatNest exclusive ownership */
668: for (i=0; i<nr*nc; i++) {
669: MatDestroy(&b[i]);
670: }
671: PetscFree(b);
673: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
674: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
675: return(0);
676: }
678: /* nest api */
679: EXTERN_C_BEGIN
682: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
683: {
684: Mat_Nest *bA = (Mat_Nest*)A->data;
686: if (idxm >= bA->nr) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
687: if (jdxm >= bA->nc) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
688: *mat = bA->m[idxm][jdxm];
689: return(0);
690: }
691: EXTERN_C_END
695: /*@C
696: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
698: Not collective
700: Input Parameters:
701: + A - nest matrix
702: . idxm - index of the matrix within the nest matrix
703: - jdxm - index of the matrix within the nest matrix
705: Output Parameter:
706: . sub - matrix at index idxm,jdxm within the nest matrix
708: Level: developer
710: .seealso: MatNestGetSize(), MatNestGetSubMats()
711: @*/
712: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
713: {
717: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
718: return(0);
719: }
721: EXTERN_C_BEGIN
724: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
725: {
726: Mat_Nest *bA = (Mat_Nest*)A->data;
727: PetscInt m,n,M,N,mi,ni,Mi,Ni;
731: if (idxm >= bA->nr) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
732: if (jdxm >= bA->nc) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
733: MatGetLocalSize(mat,&m,&n);
734: MatGetSize(mat,&M,&N);
735: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
736: ISGetSize(bA->isglobal.row[idxm],&Mi);
737: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
738: ISGetSize(bA->isglobal.col[jdxm],&Ni);
739: if (M != Mi || N != Ni) SETERRQ4(((PetscObject)mat)->comm,PETSC_ERR_ARG_INCOMP,"Submatrix dimension (%D,%D) incompatible with nest block (%D,%D)",M,N,Mi,Ni);
740: if (m != mi || n != ni) SETERRQ4(((PetscObject)mat)->comm,PETSC_ERR_ARG_INCOMP,"Submatrix local dimension (%D,%D) incompatible with nest block (%D,%D)",m,n,mi,ni);
741: PetscObjectReference((PetscObject)mat);
742: MatDestroy(&bA->m[idxm][jdxm]);
743: bA->m[idxm][jdxm] = mat;
744: return(0);
745: }
746: EXTERN_C_END
750: /*@C
751: MatNestSetSubMat - Set a single submatrix in the nest matrix.
753: Logically collective on the submatrix communicator
755: Input Parameters:
756: + A - nest matrix
757: . idxm - index of the matrix within the nest matrix
758: . jdxm - index of the matrix within the nest matrix
759: - sub - matrix at index idxm,jdxm within the nest matrix
761: Notes:
762: The new submatrix must have the same size and communicator as that block of the nest.
764: This increments the reference count of the submatrix.
766: Level: developer
768: .seealso: MatNestSetSubMats(), MatNestGetSubMat()
769: @*/
770: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
771: {
775: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
776: return(0);
777: }
779: EXTERN_C_BEGIN
782: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
783: {
784: Mat_Nest *bA = (Mat_Nest*)A->data;
786: if (M) { *M = bA->nr; }
787: if (N) { *N = bA->nc; }
788: if (mat) { *mat = bA->m; }
789: return(0);
790: }
791: EXTERN_C_END
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: }
825: EXTERN_C_BEGIN
828: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
829: {
830: Mat_Nest *bA = (Mat_Nest*)A->data;
833: if (M) { *M = bA->nr; }
834: if (N) { *N = bA->nc; }
835: return(0);
836: }
837: EXTERN_C_END
841: /*@C
842: MatNestGetSize - Returns the size of the nest matrix.
844: Not collective
846: Input Parameters:
847: . A - nest matrix
849: Output Parameter:
850: + M - number of rows in the nested mat
851: - N - number of cols in the nested mat
853: Notes:
855: Level: developer
857: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
858: @*/
859: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
860: {
864: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
865: return(0);
866: }
871: {
872: Mat_Nest *vs = (Mat_Nest*)A->data;
873: PetscInt i;
876: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
877: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
878: return(0);
879: }
883: /*@C
884: MatNestGetISs - Returns the index sets partitioning the row and column spaces
886: Not collective
888: Input Parameters:
889: . A - nest matrix
891: Output Parameter:
892: + rows - array of row index sets
893: - cols - array of column index sets
895: Level: advanced
897: Notes:
898: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
900: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
901: @*/
902: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
903: {
908: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
909: return(0);
910: }
915: {
916: Mat_Nest *vs = (Mat_Nest*)A->data;
917: PetscInt i;
920: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
921: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
922: return(0);
923: }
927: /*@C
928: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
930: Not collective
932: Input Parameters:
933: . A - nest matrix
935: Output Parameter:
936: + rows - array of row index sets (or PETSC_NULL to ignore)
937: - cols - array of column index sets (or PETSC_NULL to ignore)
939: Level: advanced
941: Notes:
942: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
944: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
945: @*/
946: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
947: {
952: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
953: return(0);
954: }
956: EXTERN_C_BEGIN
959: PetscErrorCode MatNestSetVecType_Nest(Mat A,const VecType vtype)
960: {
962: PetscBool flg;
965: PetscStrcmp(vtype,VECNEST,&flg);
966: /* In reality, this only distinguishes VECNEST and "other" */
967: if (flg) A->ops->getvecs = MatGetVecs_Nest;
968: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*))0;
969: return(0);
970: }
971: EXTERN_C_END
975: /*@C
976: MatNestSetVecType - Sets the type of Vec returned by MatGetVecs()
978: Not collective
980: Input Parameters:
981: + A - nest matrix
982: - vtype - type to use for creating vectors
984: Notes:
986: Level: developer
988: .seealso: MatGetVecs()
989: @*/
990: PetscErrorCode MatNestSetVecType(Mat A,const VecType vtype)
991: {
995: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,const VecType),(A,vtype));
996: return(0);
997: }
999: EXTERN_C_BEGIN
1002: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1003: {
1004: Mat_Nest *s = (Mat_Nest*)A->data;
1005: PetscInt i,j,m,n,M,N;
1009: s->nr = nr;
1010: s->nc = nc;
1012: /* Create space for submatrices */
1013: PetscMalloc(sizeof(Mat*)*nr,&s->m);
1014: for (i=0; i<nr; i++) {
1015: PetscMalloc(sizeof(Mat)*nc,&s->m[i]);
1016: }
1017: for (i=0; i<nr; i++) {
1018: for (j=0; j<nc; j++) {
1019: s->m[i][j] = a[i*nc+j];
1020: if (a[i*nc+j]) {
1021: PetscObjectReference((PetscObject)a[i*nc+j]);
1022: }
1023: }
1024: }
1026: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1028: PetscMalloc(sizeof(PetscInt)*nr,&s->row_len);
1029: PetscMalloc(sizeof(PetscInt)*nc,&s->col_len);
1030: for (i=0; i<nr; i++) s->row_len[i]=-1;
1031: for (j=0; j<nc; j++) s->col_len[j]=-1;
1033: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1035: PetscLayoutSetSize(A->rmap,M);
1036: PetscLayoutSetLocalSize(A->rmap,m);
1037: PetscLayoutSetSize(A->cmap,N);
1038: PetscLayoutSetLocalSize(A->cmap,n);
1040: PetscLayoutSetUp(A->rmap);
1041: PetscLayoutSetUp(A->cmap);
1043: PetscMalloc2(nr,Vec,&s->left,nc,Vec,&s->right);
1044: PetscMemzero(s->left,nr*sizeof(Vec));
1045: PetscMemzero(s->right,nc*sizeof(Vec));
1046: return(0);
1047: }
1048: EXTERN_C_END
1052: /*@
1053: MatNestSetSubMats - Sets the nested submatrices
1055: Collective on Mat
1057: Input Parameter:
1058: + N - nested matrix
1059: . nr - number of nested row blocks
1060: . is_row - index sets for each nested row block, or PETSC_NULL to make contiguous
1061: . nc - number of nested column blocks
1062: . is_col - index sets for each nested column block, or PETSC_NULL to make contiguous
1063: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using PETSC_NULL
1065: Level: advanced
1067: .seealso: MatCreateNest(), MATNEST
1068: @*/
1069: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1070: {
1072: PetscInt i;
1076: if (nr < 0) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1077: if (nr && is_row) {
1080: }
1081: if (nc < 0) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1082: if (nc && is_col) {
1085: }
1087: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1088: return(0);
1089: }
1093: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog,ISLocalToGlobalMapping *ltogb)
1094: {
1096: PetscBool flg;
1097: PetscInt i,j,m,mi,*ix;
1100: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1101: if (islocal[i]) {
1102: ISGetSize(islocal[i],&mi);
1103: flg = PETSC_TRUE; /* We found a non-trivial entry */
1104: } else {
1105: ISGetSize(isglobal[i],&mi);
1106: }
1107: m += mi;
1108: }
1109: if (flg) {
1110: PetscMalloc(m*sizeof(*ix),&ix);
1111: for (i=0,n=0; i<n; i++) {
1112: ISLocalToGlobalMapping smap = PETSC_NULL;
1113: VecScatter scat;
1114: IS isreq;
1115: Vec lvec,gvec;
1116: union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x;
1117: Mat sub;
1119: if (sizeof (*x) != sizeof(PetscScalar)) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No support when scalars smaller than integers");
1120: if (colflg) {
1121: MatNestFindNonzeroSubMatRow(A,i,&sub);
1122: } else {
1123: MatNestFindNonzeroSubMatCol(A,i,&sub);
1124: }
1125: if (sub) {MatGetLocalToGlobalMapping(sub,&smap,PETSC_NULL);}
1126: if (islocal[i]) {
1127: ISGetSize(islocal[i],&mi);
1128: } else {
1129: ISGetSize(isglobal[i],&mi);
1130: }
1131: for (j=0; j<mi; j++) ix[m+j] = j;
1132: if (smap) {ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);}
1133: /*
1134: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1135: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1136: The approach here is ugly because it uses VecScatter to move indices.
1137: */
1138: VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);
1139: VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);
1140: ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);
1141: VecScatterCreate(gvec,isreq,lvec,PETSC_NULL,&scat);
1142: VecGetArray(gvec,(PetscScalar**)&x);
1143: for (j=0; j<mi; j++) x[j].integer = ix[m+j];
1144: VecRestoreArray(gvec,(PetscScalar**)&x);
1145: VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1146: VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1147: VecGetArray(lvec,(PetscScalar**)&x);
1148: for (j=0; j<mi; j++) ix[m+j] = x[j].integer;
1149: VecRestoreArray(lvec,(PetscScalar**)&x);
1150: VecDestroy(&lvec);
1151: VecDestroy(&gvec);
1152: ISDestroy(&isreq);
1153: VecScatterDestroy(&scat);
1154: m += mi;
1155: }
1156: ISLocalToGlobalMappingCreate(((PetscObject)A)->comm,m,ix,PETSC_OWN_POINTER,ltog);
1157: *ltogb = PETSC_NULL;
1158: } else {
1159: *ltog = PETSC_NULL;
1160: *ltogb = PETSC_NULL;
1161: }
1162: return(0);
1163: }
1166: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1167: /*
1168: nprocessors = NP
1169: Nest x^T = ( (g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1) )
1170: proc 0: => (g_0,h_0,)
1171: proc 1: => (g_1,h_1,)
1172: ...
1173: proc nprocs-1: => (g_NP-1,h_NP-1,)
1175: proc 0: proc 1: proc nprocs-1:
1176: is[0] = ( 0,1,2,...,nlocal(g_0)-1 ) ( 0,1,...,nlocal(g_1)-1 ) ( 0,1,...,nlocal(g_NP-1) )
1178: proc 0:
1179: is[1] = ( nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1 )
1180: proc 1:
1181: is[1] = ( nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1 )
1183: proc NP-1:
1184: is[1] = ( nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1 )
1185: */
1188: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1189: {
1190: Mat_Nest *vs = (Mat_Nest*)A->data;
1191: PetscInt i,j,offset,n,nsum,bs;
1193: Mat sub;
1196: PetscMalloc(sizeof(IS)*nr,&vs->isglobal.row);
1197: PetscMalloc(sizeof(IS)*nc,&vs->isglobal.col);
1198: if (is_row) { /* valid IS is passed in */
1199: /* refs on is[] are incremeneted */
1200: for (i=0; i<vs->nr; i++) {
1201: PetscObjectReference((PetscObject)is_row[i]);
1202: vs->isglobal.row[i] = is_row[i];
1203: }
1204: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1205: nsum = 0;
1206: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1207: MatNestFindNonzeroSubMatRow(A,i,&sub);
1208: if (!sub) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1209: MatGetLocalSize(sub,&n,PETSC_NULL);
1210: if (n < 0) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1211: nsum += n;
1212: }
1213: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,((PetscObject)A)->comm);
1214: offset -= nsum;
1215: for (i=0; i<vs->nr; i++) {
1216: MatNestFindNonzeroSubMatRow(A,i,&sub);
1217: MatGetLocalSize(sub,&n,PETSC_NULL);
1218: MatGetBlockSize(sub,&bs);
1219: ISCreateStride(((PetscObject)sub)->comm,n,offset,1,&vs->isglobal.row[i]);
1220: ISSetBlockSize(vs->isglobal.row[i],bs);
1221: offset += n;
1222: }
1223: }
1225: if (is_col) { /* valid IS is passed in */
1226: /* refs on is[] are incremeneted */
1227: for (j=0; j<vs->nc; j++) {
1228: PetscObjectReference((PetscObject)is_col[j]);
1229: vs->isglobal.col[j] = is_col[j];
1230: }
1231: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1232: offset = A->cmap->rstart;
1233: nsum = 0;
1234: for (j=0; j<vs->nc; j++) {
1235: MatNestFindNonzeroSubMatCol(A,j,&sub);
1236: if (!sub) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1237: MatGetLocalSize(sub,PETSC_NULL,&n);
1238: if (n < 0) SETERRQ(((PetscObject)A)->comm,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,((PetscObject)A)->comm);
1242: offset -= nsum;
1243: for (j=0; j<vs->nc; j++) {
1244: MatNestFindNonzeroSubMatCol(A,j,&sub);
1245: MatGetLocalSize(sub,PETSC_NULL,&n);
1246: MatGetBlockSize(sub,&bs);
1247: ISCreateStride(((PetscObject)sub)->comm,n,offset,1,&vs->isglobal.col[j]);
1248: ISSetBlockSize(vs->isglobal.col[j],bs);
1249: offset += n;
1250: }
1251: }
1253: /* Set up the local ISs */
1254: PetscMalloc(vs->nr*sizeof(IS),&vs->islocal.row);
1255: PetscMalloc(vs->nc*sizeof(IS),&vs->islocal.col);
1256: for (i=0,offset=0; i<vs->nr; i++) {
1257: IS isloc;
1258: ISLocalToGlobalMapping rmap = PETSC_NULL;
1259: PetscInt nlocal,bs;
1260: MatNestFindNonzeroSubMatRow(A,i,&sub);
1261: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,PETSC_NULL);}
1262: if (rmap) {
1263: MatGetBlockSize(sub,&bs);
1264: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1265: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1266: ISSetBlockSize(isloc,bs);
1267: } else {
1268: nlocal = 0;
1269: isloc = PETSC_NULL;
1270: }
1271: vs->islocal.row[i] = isloc;
1272: offset += nlocal;
1273: }
1274: for (i=0,offset=0; i<vs->nc; i++) {
1275: IS isloc;
1276: ISLocalToGlobalMapping cmap = PETSC_NULL;
1277: PetscInt nlocal,bs;
1278: MatNestFindNonzeroSubMatCol(A,i,&sub);
1279: if (sub) {MatGetLocalToGlobalMapping(sub,PETSC_NULL,&cmap);}
1280: if (cmap) {
1281: MatGetBlockSize(sub,&bs);
1282: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1283: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1284: ISSetBlockSize(isloc,bs);
1285: } else {
1286: nlocal = 0;
1287: isloc = PETSC_NULL;
1288: }
1289: vs->islocal.col[i] = isloc;
1290: offset += nlocal;
1291: }
1293: /* Set up the aggregate ISLocalToGlobalMapping */
1294: {
1295: ISLocalToGlobalMapping rmap,rmapb,cmap,cmapb;
1296: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap,&rmapb);
1297: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap,&cmapb);
1298: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1299: if (rmapb && cmapb) {MatSetLocalToGlobalMappingBlock(A,rmapb,cmapb);}
1300: ISLocalToGlobalMappingDestroy(&rmap);
1301: ISLocalToGlobalMappingDestroy(&rmapb);
1302: ISLocalToGlobalMappingDestroy(&cmap);
1303: ISLocalToGlobalMappingDestroy(&cmapb);
1304: }
1306: #if defined(PETSC_USE_DEBUG)
1307: for (i=0; i<vs->nr; i++) {
1308: for (j=0; j<vs->nc; j++) {
1309: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1310: Mat B = vs->m[i][j];
1311: if (!B) continue;
1312: MatGetSize(B,&M,&N);
1313: MatGetLocalSize(B,&m,&n);
1314: ISGetSize(vs->isglobal.row[i],&Mi);
1315: ISGetSize(vs->isglobal.col[j],&Ni);
1316: ISGetLocalSize(vs->isglobal.row[i],&mi);
1317: ISGetLocalSize(vs->isglobal.col[j],&ni);
1318: if (M != Mi || N != Ni) SETERRQ6(((PetscObject)sub)->comm,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);
1319: if (m != mi || n != ni) SETERRQ6(((PetscObject)sub)->comm,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);
1320: }
1321: }
1322: #endif
1324: /* Set A->assembled if all non-null blocks are currently assembled */
1325: for (i=0; i<vs->nr; i++) {
1326: for (j=0; j<vs->nc; j++) {
1327: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1328: }
1329: }
1330: A->assembled = PETSC_TRUE;
1331: return(0);
1332: }
1336: /*@
1337: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1339: Collective on Mat
1341: Input Parameter:
1342: + comm - Communicator for the new Mat
1343: . nr - number of nested row blocks
1344: . is_row - index sets for each nested row block, or PETSC_NULL to make contiguous
1345: . nc - number of nested column blocks
1346: . is_col - index sets for each nested column block, or PETSC_NULL to make contiguous
1347: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using PETSC_NULL
1349: Output Parameter:
1350: . B - new matrix
1352: Level: advanced
1354: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1355: @*/
1356: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1357: {
1358: Mat A;
1362: *B = 0;
1363: MatCreate(comm,&A);
1364: MatSetType(A,MATNEST);
1365: MatSetUp(A);
1366: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1367: *B = A;
1368: return(0);
1369: }
1371: /*MC
1372: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
1374: Level: intermediate
1376: Notes:
1377: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
1378: It allows the use of symmetric and block formats for parts of multi-physics simulations.
1379: It is usually used with DMComposite and DMCreateMatrix()
1381: .seealso: MatCreate(), MatType, MatCreateNest()
1382: M*/
1383: EXTERN_C_BEGIN
1386: PetscErrorCode MatCreate_Nest(Mat A)
1387: {
1388: Mat_Nest *s;
1392: PetscNewLog(A,Mat_Nest,&s);
1393: A->data = (void*)s;
1395: s->nr = -1;
1396: s->nc = -1;
1397: s->m = PETSC_NULL;
1398: s->splitassembly = PETSC_FALSE;
1400: PetscMemzero(A->ops,sizeof(*A->ops));
1401: A->ops->mult = MatMult_Nest;
1402: A->ops->multadd = MatMultAdd_Nest;
1403: A->ops->multtranspose = MatMultTranspose_Nest;
1404: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
1405: A->ops->assemblybegin = MatAssemblyBegin_Nest;
1406: A->ops->assemblyend = MatAssemblyEnd_Nest;
1407: A->ops->zeroentries = MatZeroEntries_Nest;
1408: A->ops->duplicate = MatDuplicate_Nest;
1409: A->ops->getsubmatrix = MatGetSubMatrix_Nest;
1410: A->ops->destroy = MatDestroy_Nest;
1411: A->ops->view = MatView_Nest;
1412: A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1413: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
1414: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1415: A->ops->getdiagonal = MatGetDiagonal_Nest;
1416: A->ops->diagonalscale = MatDiagonalScale_Nest;
1417: A->ops->scale = MatScale_Nest;
1418: A->ops->shift = MatShift_Nest;
1420: A->spptr = 0;
1421: A->same_nonzero = PETSC_FALSE;
1422: A->assembled = PETSC_FALSE;
1424: /* expose Nest api's */
1425: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMat_C", "MatNestGetSubMat_Nest", MatNestGetSubMat_Nest);
1426: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestSetSubMat_C", "MatNestSetSubMat_Nest", MatNestSetSubMat_Nest);
1427: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMats_C", "MatNestGetSubMats_Nest", MatNestGetSubMats_Nest);
1428: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSize_C", "MatNestGetSize_Nest", MatNestGetSize_Nest);
1429: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetISs_C", "MatNestGetISs_Nest", MatNestGetISs_Nest);
1430: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetLocalISs_C", "MatNestGetLocalISs_Nest", MatNestGetLocalISs_Nest);
1431: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestSetVecType_C", "MatNestSetVecType_Nest", MatNestSetVecType_Nest);
1432: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestSetSubMats_C", "MatNestSetSubMats_Nest", MatNestSetSubMats_Nest);
1434: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1435: return(0);
1436: }
1437: EXTERN_C_END