Actual source code: matnest.c
petsc-3.12.5 2020-03-29
2: #include <../src/mat/impls/nest/matnestimpl.h>
3: #include <../src/mat/impls/aij/seq/aij.h>
4: #include <petscsf.h>
6: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
7: static PetscErrorCode MatCreateVecs_Nest(Mat,Vec*,Vec*);
8: static PetscErrorCode MatReset_Nest(Mat);
10: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*);
12: /* private functions */
13: static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
14: {
15: Mat_Nest *bA = (Mat_Nest*)A->data;
16: PetscInt i,j;
20: *m = *n = *M = *N = 0;
21: for (i=0; i<bA->nr; i++) { /* rows */
22: PetscInt sm,sM;
23: ISGetLocalSize(bA->isglobal.row[i],&sm);
24: ISGetSize(bA->isglobal.row[i],&sM);
25: *m += sm;
26: *M += sM;
27: }
28: for (j=0; j<bA->nc; j++) { /* cols */
29: PetscInt sn,sN;
30: ISGetLocalSize(bA->isglobal.col[j],&sn);
31: ISGetSize(bA->isglobal.col[j],&sN);
32: *n += sn;
33: *N += sN;
34: }
35: return(0);
36: }
38: /* 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: }
62: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
63: {
64: Mat_Nest *bA = (Mat_Nest*)A->data;
65: Vec *bx = bA->right,*bz = bA->left;
66: PetscInt i,j,nr = bA->nr,nc = bA->nc;
70: for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
71: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
72: for (i=0; i<nr; i++) {
73: if (y != z) {
74: Vec by;
75: VecGetSubVector(y,bA->isglobal.row[i],&by);
76: VecCopy(by,bz[i]);
77: VecRestoreSubVector(y,bA->isglobal.row[i],&by);
78: }
79: for (j=0; j<nc; j++) {
80: if (!bA->m[i][j]) continue;
81: /* y[i] <- y[i] + A[i][j] * x[j] */
82: MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
83: }
84: }
85: for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
86: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
87: return(0);
88: }
90: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
91: {
92: Mat_Nest *bA = (Mat_Nest*)A->data;
93: Vec *bx = bA->left,*by = bA->right;
94: PetscInt i,j,nr = bA->nr,nc = bA->nc;
98: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
99: for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
100: for (j=0; j<nc; j++) {
101: VecZeroEntries(by[j]);
102: for (i=0; i<nr; i++) {
103: if (!bA->m[i][j]) continue;
104: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
105: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
106: }
107: }
108: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
109: for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
110: return(0);
111: }
113: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
114: {
115: Mat_Nest *bA = (Mat_Nest*)A->data;
116: Vec *bx = bA->left,*bz = bA->right;
117: PetscInt i,j,nr = bA->nr,nc = bA->nc;
121: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
122: for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
123: for (j=0; j<nc; j++) {
124: if (y != z) {
125: Vec by;
126: VecGetSubVector(y,bA->isglobal.col[j],&by);
127: VecCopy(by,bz[j]);
128: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
129: }
130: for (i=0; i<nr; i++) {
131: if (!bA->m[i][j]) continue;
132: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
133: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
134: }
135: }
136: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
137: for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
138: return(0);
139: }
141: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
142: {
143: Mat_Nest *bA = (Mat_Nest*)A->data, *bC;
144: Mat C;
145: PetscInt i,j,nr = bA->nr,nc = bA->nc;
149: if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");
151: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
152: Mat *subs;
153: IS *is_row,*is_col;
155: PetscCalloc1(nr * nc,&subs);
156: PetscMalloc2(nr,&is_row,nc,&is_col);
157: MatNestGetISs(A,is_row,is_col);
158: if (reuse == MAT_INPLACE_MATRIX) {
159: for (i=0; i<nr; i++) {
160: for (j=0; j<nc; j++) {
161: subs[i + nr * j] = bA->m[i][j];
162: }
163: }
164: }
166: MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
167: PetscFree(subs);
168: PetscFree2(is_row,is_col);
169: } else {
170: C = *B;
171: }
173: bC = (Mat_Nest*)C->data;
174: for (i=0; i<nr; i++) {
175: for (j=0; j<nc; j++) {
176: if (bA->m[i][j]) {
177: MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
178: } else {
179: bC->m[j][i] = NULL;
180: }
181: }
182: }
184: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
185: *B = C;
186: } else {
187: MatHeaderMerge(A, &C);
188: }
189: return(0);
190: }
192: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
193: {
195: IS *lst = *list;
196: PetscInt i;
199: if (!lst) return(0);
200: for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
201: PetscFree(lst);
202: *list = NULL;
203: return(0);
204: }
206: static PetscErrorCode MatReset_Nest(Mat A)
207: {
208: Mat_Nest *vs = (Mat_Nest*)A->data;
209: PetscInt i,j;
213: /* release the matrices and the place holders */
214: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
215: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
216: MatNestDestroyISList(vs->nr,&vs->islocal.row);
217: MatNestDestroyISList(vs->nc,&vs->islocal.col);
219: PetscFree(vs->row_len);
220: PetscFree(vs->col_len);
221: PetscFree(vs->nnzstate);
223: PetscFree2(vs->left,vs->right);
225: /* release the matrices and the place holders */
226: if (vs->m) {
227: for (i=0; i<vs->nr; i++) {
228: for (j=0; j<vs->nc; j++) {
229: MatDestroy(&vs->m[i][j]);
230: }
231: PetscFree(vs->m[i]);
232: }
233: PetscFree(vs->m);
234: }
236: /* restore defaults */
237: vs->nr = 0;
238: vs->nc = 0;
239: vs->splitassembly = PETSC_FALSE;
240: return(0);
241: }
243: static PetscErrorCode MatDestroy_Nest(Mat A)
244: {
247: MatReset_Nest(A);
248: PetscFree(A->data);
250: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
251: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
252: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
253: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
254: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
255: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
256: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
257: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
258: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
259: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
260: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
261: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
262: return(0);
263: }
265: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
266: {
267: Mat_Nest *vs = (Mat_Nest*)A->data;
268: PetscInt i,j;
270: PetscBool nnzstate = PETSC_FALSE;
273: for (i=0; i<vs->nr; i++) {
274: for (j=0; j<vs->nc; j++) {
275: PetscObjectState subnnzstate = 0;
276: if (vs->m[i][j]) {
277: MatAssemblyBegin(vs->m[i][j],type);
278: if (!vs->splitassembly) {
279: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
280: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
281: * already performing an assembly, but the result would by more complicated and appears to offer less
282: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
283: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
284: */
285: MatAssemblyEnd(vs->m[i][j],type);
286: MatGetNonzeroState(vs->m[i][j],&subnnzstate);
287: }
288: }
289: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
290: vs->nnzstate[i*vs->nc+j] = subnnzstate;
291: }
292: }
293: if (nnzstate) A->nonzerostate++;
294: return(0);
295: }
297: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
298: {
299: Mat_Nest *vs = (Mat_Nest*)A->data;
300: PetscInt i,j;
304: for (i=0; i<vs->nr; i++) {
305: for (j=0; j<vs->nc; j++) {
306: if (vs->m[i][j]) {
307: if (vs->splitassembly) {
308: MatAssemblyEnd(vs->m[i][j],type);
309: }
310: }
311: }
312: }
313: return(0);
314: }
316: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
317: {
319: Mat_Nest *vs = (Mat_Nest*)A->data;
320: PetscInt j;
321: Mat sub;
324: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
325: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
326: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
327: *B = sub;
328: return(0);
329: }
331: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
332: {
334: Mat_Nest *vs = (Mat_Nest*)A->data;
335: PetscInt i;
336: Mat sub;
339: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
340: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
341: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
342: *B = sub;
343: return(0);
344: }
346: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
347: {
349: PetscInt i;
350: PetscBool flg;
356: *found = -1;
357: for (i=0; i<n; i++) {
358: if (!list[i]) continue;
359: ISEqualUnsorted(list[i],is,&flg);
360: if (flg) {
361: *found = i;
362: return(0);
363: }
364: }
365: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
366: return(0);
367: }
369: /* Get a block row as a new MatNest */
370: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
371: {
372: Mat_Nest *vs = (Mat_Nest*)A->data;
373: char keyname[256];
377: *B = NULL;
378: PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
379: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
380: if (*B) return(0);
382: MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);
384: (*B)->assembled = A->assembled;
386: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
387: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
388: return(0);
389: }
391: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
392: {
393: Mat_Nest *vs = (Mat_Nest*)A->data;
395: PetscInt row,col;
396: PetscBool same,isFullCol,isFullColGlobal;
399: /* Check if full column space. This is a hack */
400: isFullCol = PETSC_FALSE;
401: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
402: if (same) {
403: PetscInt n,first,step,i,an,am,afirst,astep;
404: ISStrideGetInfo(iscol,&first,&step);
405: ISGetLocalSize(iscol,&n);
406: isFullCol = PETSC_TRUE;
407: for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
408: PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
409: ISGetLocalSize(is->col[i],&am);
410: if (same) {
411: ISStrideGetInfo(is->col[i],&afirst,&astep);
412: if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
413: } else isFullCol = PETSC_FALSE;
414: an += am;
415: }
416: if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
417: }
418: MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));
420: if (isFullColGlobal && vs->nc > 1) {
421: PetscInt row;
422: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
423: MatNestGetRow(A,row,B);
424: } else {
425: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
426: MatNestFindIS(A,vs->nc,is->col,iscol,&col);
427: if (!vs->m[row][col]) {
428: PetscInt lr,lc;
430: MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
431: ISGetLocalSize(vs->isglobal.row[row],&lr);
432: ISGetLocalSize(vs->isglobal.col[col],&lc);
433: MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
434: MatSetType(vs->m[row][col],MATAIJ);
435: MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
436: MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
437: MatSetUp(vs->m[row][col]);
438: MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
439: MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
440: }
441: *B = vs->m[row][col];
442: }
443: return(0);
444: }
446: /*
447: TODO: This does not actually returns a submatrix we can modify
448: */
449: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
450: {
452: Mat_Nest *vs = (Mat_Nest*)A->data;
453: Mat sub;
456: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
457: switch (reuse) {
458: case MAT_INITIAL_MATRIX:
459: if (sub) { PetscObjectReference((PetscObject)sub); }
460: *B = sub;
461: break;
462: case MAT_REUSE_MATRIX:
463: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
464: break;
465: case MAT_IGNORE_MATRIX: /* Nothing to do */
466: break;
467: case MAT_INPLACE_MATRIX: /* Nothing to do */
468: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
469: break;
470: }
471: return(0);
472: }
474: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
475: {
477: Mat_Nest *vs = (Mat_Nest*)A->data;
478: Mat sub;
481: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
482: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
483: if (sub) {PetscObjectReference((PetscObject)sub);}
484: *B = sub;
485: return(0);
486: }
488: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
489: {
491: Mat_Nest *vs = (Mat_Nest*)A->data;
492: Mat sub;
495: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
496: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
497: if (sub) {
498: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
499: MatDestroy(B);
500: }
501: return(0);
502: }
504: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
505: {
506: Mat_Nest *bA = (Mat_Nest*)A->data;
507: PetscInt i;
511: for (i=0; i<bA->nr; i++) {
512: Vec bv;
513: VecGetSubVector(v,bA->isglobal.row[i],&bv);
514: if (bA->m[i][i]) {
515: MatGetDiagonal(bA->m[i][i],bv);
516: } else {
517: VecSet(bv,0.0);
518: }
519: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
520: }
521: return(0);
522: }
524: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
525: {
526: Mat_Nest *bA = (Mat_Nest*)A->data;
527: Vec bl,*br;
528: PetscInt i,j;
532: PetscCalloc1(bA->nc,&br);
533: if (r) {
534: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
535: }
536: bl = NULL;
537: for (i=0; i<bA->nr; i++) {
538: if (l) {
539: VecGetSubVector(l,bA->isglobal.row[i],&bl);
540: }
541: for (j=0; j<bA->nc; j++) {
542: if (bA->m[i][j]) {
543: MatDiagonalScale(bA->m[i][j],bl,br[j]);
544: }
545: }
546: if (l) {
547: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
548: }
549: }
550: if (r) {
551: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
552: }
553: PetscFree(br);
554: return(0);
555: }
557: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
558: {
559: Mat_Nest *bA = (Mat_Nest*)A->data;
560: PetscInt i,j;
564: for (i=0; i<bA->nr; i++) {
565: for (j=0; j<bA->nc; j++) {
566: if (bA->m[i][j]) {
567: MatScale(bA->m[i][j],a);
568: }
569: }
570: }
571: return(0);
572: }
574: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
575: {
576: Mat_Nest *bA = (Mat_Nest*)A->data;
577: PetscInt i;
579: PetscBool nnzstate = PETSC_FALSE;
582: for (i=0; i<bA->nr; i++) {
583: PetscObjectState subnnzstate = 0;
584: 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);
585: MatShift(bA->m[i][i],a);
586: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
587: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
588: bA->nnzstate[i*bA->nc+i] = subnnzstate;
589: }
590: if (nnzstate) A->nonzerostate++;
591: return(0);
592: }
594: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
595: {
596: Mat_Nest *bA = (Mat_Nest*)A->data;
597: PetscInt i;
599: PetscBool nnzstate = PETSC_FALSE;
602: for (i=0; i<bA->nr; i++) {
603: PetscObjectState subnnzstate = 0;
604: Vec bv;
605: VecGetSubVector(D,bA->isglobal.row[i],&bv);
606: if (bA->m[i][i]) {
607: MatDiagonalSet(bA->m[i][i],bv,is);
608: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
609: }
610: VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
611: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
612: bA->nnzstate[i*bA->nc+i] = subnnzstate;
613: }
614: if (nnzstate) A->nonzerostate++;
615: return(0);
616: }
618: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
619: {
620: Mat_Nest *bA = (Mat_Nest*)A->data;
621: PetscInt i,j;
625: for (i=0; i<bA->nr; i++) {
626: for (j=0; j<bA->nc; j++) {
627: if (bA->m[i][j]) {
628: MatSetRandom(bA->m[i][j],rctx);
629: }
630: }
631: }
632: return(0);
633: }
635: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
636: {
637: Mat_Nest *bA = (Mat_Nest*)A->data;
638: Vec *L,*R;
639: MPI_Comm comm;
640: PetscInt i,j;
644: PetscObjectGetComm((PetscObject)A,&comm);
645: if (right) {
646: /* allocate R */
647: PetscMalloc1(bA->nc, &R);
648: /* Create the right vectors */
649: for (j=0; j<bA->nc; j++) {
650: for (i=0; i<bA->nr; i++) {
651: if (bA->m[i][j]) {
652: MatCreateVecs(bA->m[i][j],&R[j],NULL);
653: break;
654: }
655: }
656: if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
657: }
658: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
659: /* hand back control to the nest vector */
660: for (j=0; j<bA->nc; j++) {
661: VecDestroy(&R[j]);
662: }
663: PetscFree(R);
664: }
666: if (left) {
667: /* allocate L */
668: PetscMalloc1(bA->nr, &L);
669: /* Create the left vectors */
670: for (i=0; i<bA->nr; i++) {
671: for (j=0; j<bA->nc; j++) {
672: if (bA->m[i][j]) {
673: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
674: break;
675: }
676: }
677: if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
678: }
680: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
681: for (i=0; i<bA->nr; i++) {
682: VecDestroy(&L[i]);
683: }
685: PetscFree(L);
686: }
687: return(0);
688: }
690: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
691: {
692: Mat_Nest *bA = (Mat_Nest*)A->data;
693: PetscBool isascii,viewSub = PETSC_FALSE;
694: PetscInt i,j;
698: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
699: if (isascii) {
701: PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
702: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
703: PetscViewerASCIIPushTab(viewer);
704: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);
706: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
707: for (i=0; i<bA->nr; i++) {
708: for (j=0; j<bA->nc; j++) {
709: MatType type;
710: char name[256] = "",prefix[256] = "";
711: PetscInt NR,NC;
712: PetscBool isNest = PETSC_FALSE;
714: if (!bA->m[i][j]) {
715: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
716: continue;
717: }
718: MatGetSize(bA->m[i][j],&NR,&NC);
719: MatGetType(bA->m[i][j], &type);
720: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
721: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
722: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
724: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
726: if (isNest || viewSub) {
727: PetscViewerASCIIPushTab(viewer); /* push1 */
728: MatView(bA->m[i][j],viewer);
729: PetscViewerASCIIPopTab(viewer); /* pop1 */
730: }
731: }
732: }
733: PetscViewerASCIIPopTab(viewer); /* pop0 */
734: }
735: return(0);
736: }
738: static PetscErrorCode MatZeroEntries_Nest(Mat A)
739: {
740: Mat_Nest *bA = (Mat_Nest*)A->data;
741: PetscInt i,j;
745: for (i=0; i<bA->nr; i++) {
746: for (j=0; j<bA->nc; j++) {
747: if (!bA->m[i][j]) continue;
748: MatZeroEntries(bA->m[i][j]);
749: }
750: }
751: return(0);
752: }
754: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
755: {
756: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
757: PetscInt i,j,nr = bA->nr,nc = bA->nc;
759: PetscBool nnzstate = PETSC_FALSE;
762: 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);
763: for (i=0; i<nr; i++) {
764: for (j=0; j<nc; j++) {
765: PetscObjectState subnnzstate = 0;
766: if (bA->m[i][j] && bB->m[i][j]) {
767: MatCopy(bA->m[i][j],bB->m[i][j],str);
768: } 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);
769: MatGetNonzeroState(bB->m[i][j],&subnnzstate);
770: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
771: bB->nnzstate[i*nc+j] = subnnzstate;
772: }
773: }
774: if (nnzstate) B->nonzerostate++;
775: return(0);
776: }
778: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
779: {
780: Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
781: PetscInt i,j,nr = bY->nr,nc = bY->nc;
783: PetscBool nnzstate = PETSC_FALSE;
786: if (nr != bX->nr || nc != bX->nc) SETERRQ4(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Cannot AXPY a MatNest of block size (%D,%D) with a MatNest of block size (%D,%D)",bX->nr,bX->nc,nr,nc);
787: for (i=0; i<nr; i++) {
788: for (j=0; j<nc; j++) {
789: PetscObjectState subnnzstate = 0;
790: if (bY->m[i][j] && bX->m[i][j]) {
791: MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
792: } else if (bX->m[i][j]) {
793: Mat M;
795: if (str != DIFFERENT_NONZERO_PATTERN) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D. Use DIFFERENT_NONZERO_PATTERN",i,j);
796: MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
797: MatNestSetSubMat(Y,i,j,M);
798: MatDestroy(&M);
799: }
800: if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
801: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
802: bY->nnzstate[i*nc+j] = subnnzstate;
803: }
804: }
805: if (nnzstate) Y->nonzerostate++;
806: return(0);
807: }
809: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
810: {
811: Mat_Nest *bA = (Mat_Nest*)A->data;
812: Mat *b;
813: PetscInt i,j,nr = bA->nr,nc = bA->nc;
817: PetscMalloc1(nr*nc,&b);
818: for (i=0; i<nr; i++) {
819: for (j=0; j<nc; j++) {
820: if (bA->m[i][j]) {
821: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
822: } else {
823: b[i*nc+j] = NULL;
824: }
825: }
826: }
827: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
828: /* Give the new MatNest exclusive ownership */
829: for (i=0; i<nr*nc; i++) {
830: MatDestroy(&b[i]);
831: }
832: PetscFree(b);
834: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
835: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
836: return(0);
837: }
839: /* nest api */
840: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
841: {
842: Mat_Nest *bA = (Mat_Nest*)A->data;
845: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
846: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
847: *mat = bA->m[idxm][jdxm];
848: return(0);
849: }
851: /*@
852: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
854: Not collective
856: Input Parameters:
857: + A - nest matrix
858: . idxm - index of the matrix within the nest matrix
859: - jdxm - index of the matrix within the nest matrix
861: Output Parameter:
862: . sub - matrix at index idxm,jdxm within the nest matrix
864: Level: developer
866: .seealso: MatNestGetSize(), MatNestGetSubMats()
867: @*/
868: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
869: {
873: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
874: return(0);
875: }
877: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
878: {
879: Mat_Nest *bA = (Mat_Nest*)A->data;
880: PetscInt m,n,M,N,mi,ni,Mi,Ni;
884: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
885: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
886: MatGetLocalSize(mat,&m,&n);
887: MatGetSize(mat,&M,&N);
888: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
889: ISGetSize(bA->isglobal.row[idxm],&Mi);
890: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
891: ISGetSize(bA->isglobal.col[jdxm],&Ni);
892: 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);
893: 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);
895: /* do not increase object state */
896: if (mat == bA->m[idxm][jdxm]) return(0);
898: PetscObjectReference((PetscObject)mat);
899: MatDestroy(&bA->m[idxm][jdxm]);
900: bA->m[idxm][jdxm] = mat;
901: PetscObjectStateIncrease((PetscObject)A);
902: MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
903: A->nonzerostate++;
904: return(0);
905: }
907: /*@
908: MatNestSetSubMat - Set a single submatrix in the nest matrix.
910: Logically collective on the submatrix communicator
912: Input Parameters:
913: + A - nest matrix
914: . idxm - index of the matrix within the nest matrix
915: . jdxm - index of the matrix within the nest matrix
916: - sub - matrix at index idxm,jdxm within the nest matrix
918: Notes:
919: The new submatrix must have the same size and communicator as that block of the nest.
921: This increments the reference count of the submatrix.
923: Level: developer
925: .seealso: MatNestSetSubMats(), MatNestGetSubMats()
926: @*/
927: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
928: {
932: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
933: return(0);
934: }
936: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
937: {
938: Mat_Nest *bA = (Mat_Nest*)A->data;
941: if (M) *M = bA->nr;
942: if (N) *N = bA->nc;
943: if (mat) *mat = bA->m;
944: return(0);
945: }
947: /*@C
948: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
950: Not collective
952: Input Parameters:
953: . A - nest matrix
955: Output Parameter:
956: + M - number of rows in the nest matrix
957: . N - number of cols in the nest matrix
958: - mat - 2d array of matrices
960: Notes:
962: The user should not free the array mat.
964: In Fortran, this routine has a calling sequence
965: $ call MatNestGetSubMats(A, M, N, mat, ierr)
966: where the space allocated for the optional argument mat is assumed large enough (if provided).
968: Level: developer
970: .seealso: MatNestGetSize(), MatNestGetSubMat()
971: @*/
972: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
973: {
977: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
978: return(0);
979: }
981: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
982: {
983: Mat_Nest *bA = (Mat_Nest*)A->data;
986: if (M) *M = bA->nr;
987: if (N) *N = bA->nc;
988: return(0);
989: }
991: /*@
992: MatNestGetSize - Returns the size of the nest matrix.
994: Not collective
996: Input Parameters:
997: . A - nest matrix
999: Output Parameter:
1000: + M - number of rows in the nested mat
1001: - N - number of cols in the nested mat
1003: Notes:
1005: Level: developer
1007: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
1008: @*/
1009: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1010: {
1014: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1015: return(0);
1016: }
1018: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1019: {
1020: Mat_Nest *vs = (Mat_Nest*)A->data;
1021: PetscInt i;
1024: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1025: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1026: return(0);
1027: }
1029: /*@C
1030: MatNestGetISs - Returns the index sets partitioning the row and column spaces
1032: Not collective
1034: Input Parameters:
1035: . A - nest matrix
1037: Output Parameter:
1038: + rows - array of row index sets
1039: - cols - array of column index sets
1041: Level: advanced
1043: Notes:
1044: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1046: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
1047: @*/
1048: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1049: {
1054: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1055: return(0);
1056: }
1058: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1059: {
1060: Mat_Nest *vs = (Mat_Nest*)A->data;
1061: PetscInt i;
1064: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1065: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1066: return(0);
1067: }
1069: /*@C
1070: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1072: Not collective
1074: Input Parameters:
1075: . A - nest matrix
1077: Output Parameter:
1078: + rows - array of row index sets (or NULL to ignore)
1079: - cols - array of column index sets (or NULL to ignore)
1081: Level: advanced
1083: Notes:
1084: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1086: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
1087: @*/
1088: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1089: {
1094: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1095: return(0);
1096: }
1098: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1099: {
1101: PetscBool flg;
1104: PetscStrcmp(vtype,VECNEST,&flg);
1105: /* In reality, this only distinguishes VECNEST and "other" */
1106: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1107: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1108: return(0);
1109: }
1111: /*@C
1112: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1114: Not collective
1116: Input Parameters:
1117: + A - nest matrix
1118: - vtype - type to use for creating vectors
1120: Notes:
1122: Level: developer
1124: .seealso: MatCreateVecs()
1125: @*/
1126: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1127: {
1131: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1132: return(0);
1133: }
1135: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1136: {
1137: Mat_Nest *s = (Mat_Nest*)A->data;
1138: PetscInt i,j,m,n,M,N;
1140: PetscBool cong;
1143: MatReset_Nest(A);
1145: s->nr = nr;
1146: s->nc = nc;
1148: /* Create space for submatrices */
1149: PetscMalloc1(nr,&s->m);
1150: for (i=0; i<nr; i++) {
1151: PetscMalloc1(nc,&s->m[i]);
1152: }
1153: for (i=0; i<nr; i++) {
1154: for (j=0; j<nc; j++) {
1155: s->m[i][j] = a[i*nc+j];
1156: if (a[i*nc+j]) {
1157: PetscObjectReference((PetscObject)a[i*nc+j]);
1158: }
1159: }
1160: }
1162: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1164: PetscMalloc1(nr,&s->row_len);
1165: PetscMalloc1(nc,&s->col_len);
1166: for (i=0; i<nr; i++) s->row_len[i]=-1;
1167: for (j=0; j<nc; j++) s->col_len[j]=-1;
1169: PetscCalloc1(nr*nc,&s->nnzstate);
1170: for (i=0; i<nr; i++) {
1171: for (j=0; j<nc; j++) {
1172: if (s->m[i][j]) {
1173: MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1174: }
1175: }
1176: }
1178: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1180: PetscLayoutSetSize(A->rmap,M);
1181: PetscLayoutSetLocalSize(A->rmap,m);
1182: PetscLayoutSetSize(A->cmap,N);
1183: PetscLayoutSetLocalSize(A->cmap,n);
1185: PetscLayoutSetUp(A->rmap);
1186: PetscLayoutSetUp(A->cmap);
1188: /* disable operations that are not supported for non-square matrices,
1189: or matrices for which is_row != is_col */
1190: MatHasCongruentLayouts(A,&cong);
1191: if (cong && nr != nc) cong = PETSC_FALSE;
1192: if (cong) {
1193: for (i = 0; cong && i < nr; i++) {
1194: ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1195: }
1196: }
1197: if (!cong) {
1198: A->ops->getdiagonal = NULL;
1199: A->ops->shift = NULL;
1200: A->ops->diagonalset = NULL;
1201: }
1203: PetscCalloc2(nr,&s->left,nc,&s->right);
1204: PetscObjectStateIncrease((PetscObject)A);
1205: A->nonzerostate++;
1206: return(0);
1207: }
1209: /*@
1210: MatNestSetSubMats - Sets the nested submatrices
1212: Collective on Mat
1214: Input Parameter:
1215: + A - nested matrix
1216: . nr - number of nested row blocks
1217: . is_row - index sets for each nested row block, or NULL to make contiguous
1218: . nc - number of nested column blocks
1219: . is_col - index sets for each nested column block, or NULL to make contiguous
1220: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1222: Notes: this always resets any submatrix information previously set
1224: Level: advanced
1226: .seealso: MatCreateNest(), MATNEST
1227: @*/
1228: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1229: {
1231: PetscInt i;
1235: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1236: if (nr && is_row) {
1239: }
1240: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1241: if (nc && is_col) {
1244: }
1246: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1247: return(0);
1248: }
1250: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1251: {
1253: PetscBool flg;
1254: PetscInt i,j,m,mi,*ix;
1257: *ltog = NULL;
1258: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1259: if (islocal[i]) {
1260: ISGetLocalSize(islocal[i],&mi);
1261: flg = PETSC_TRUE; /* We found a non-trivial entry */
1262: } else {
1263: ISGetLocalSize(isglobal[i],&mi);
1264: }
1265: m += mi;
1266: }
1267: if (!flg) return(0);
1269: PetscMalloc1(m,&ix);
1270: for (i=0,m=0; i<n; i++) {
1271: ISLocalToGlobalMapping smap = NULL;
1272: Mat sub = NULL;
1273: PetscSF sf;
1274: PetscLayout map;
1275: const PetscInt *ix2;
1277: if (!colflg) {
1278: MatNestFindNonzeroSubMatRow(A,i,&sub);
1279: } else {
1280: MatNestFindNonzeroSubMatCol(A,i,&sub);
1281: }
1282: if (sub) {
1283: if (!colflg) {
1284: MatGetLocalToGlobalMapping(sub,&smap,NULL);
1285: } else {
1286: MatGetLocalToGlobalMapping(sub,NULL,&smap);
1287: }
1288: }
1289: /*
1290: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1291: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1292: */
1293: ISGetIndices(isglobal[i],&ix2);
1294: if (islocal[i]) {
1295: PetscInt *ilocal,*iremote;
1296: PetscInt mil,nleaves;
1298: ISGetLocalSize(islocal[i],&mi);
1299: if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1300: for (j=0; j<mi; j++) ix[m+j] = j;
1301: ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);
1303: /* PetscSFSetGraphLayout does not like negative indices */
1304: PetscMalloc2(mi,&ilocal,mi,&iremote);
1305: for (j=0, nleaves = 0; j<mi; j++) {
1306: if (ix[m+j] < 0) continue;
1307: ilocal[nleaves] = j;
1308: iremote[nleaves] = ix[m+j];
1309: nleaves++;
1310: }
1311: ISGetLocalSize(isglobal[i],&mil);
1312: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1313: PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1314: PetscLayoutSetLocalSize(map,mil);
1315: PetscLayoutSetUp(map);
1316: PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1317: PetscLayoutDestroy(&map);
1318: PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1319: PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1320: PetscSFDestroy(&sf);
1321: PetscFree2(ilocal,iremote);
1322: } else {
1323: ISGetLocalSize(isglobal[i],&mi);
1324: for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1325: }
1326: ISRestoreIndices(isglobal[i],&ix2);
1327: m += mi;
1328: }
1329: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1330: return(0);
1331: }
1334: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1335: /*
1336: nprocessors = NP
1337: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1338: proc 0: => (g_0,h_0,)
1339: proc 1: => (g_1,h_1,)
1340: ...
1341: proc nprocs-1: => (g_NP-1,h_NP-1,)
1343: proc 0: proc 1: proc nprocs-1:
1344: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1346: proc 0:
1347: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1348: proc 1:
1349: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1351: proc NP-1:
1352: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1353: */
1354: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1355: {
1356: Mat_Nest *vs = (Mat_Nest*)A->data;
1357: PetscInt i,j,offset,n,nsum,bs;
1359: Mat sub = NULL;
1362: PetscMalloc1(nr,&vs->isglobal.row);
1363: PetscMalloc1(nc,&vs->isglobal.col);
1364: if (is_row) { /* valid IS is passed in */
1365: /* refs on is[] are incremeneted */
1366: for (i=0; i<vs->nr; i++) {
1367: PetscObjectReference((PetscObject)is_row[i]);
1369: vs->isglobal.row[i] = is_row[i];
1370: }
1371: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1372: nsum = 0;
1373: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1374: MatNestFindNonzeroSubMatRow(A,i,&sub);
1375: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1376: MatGetLocalSize(sub,&n,NULL);
1377: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1378: nsum += n;
1379: }
1380: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1381: offset -= nsum;
1382: for (i=0; i<vs->nr; i++) {
1383: MatNestFindNonzeroSubMatRow(A,i,&sub);
1384: MatGetLocalSize(sub,&n,NULL);
1385: MatGetBlockSize(sub,&bs);
1386: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1387: ISSetBlockSize(vs->isglobal.row[i],bs);
1388: offset += n;
1389: }
1390: }
1392: if (is_col) { /* valid IS is passed in */
1393: /* refs on is[] are incremeneted */
1394: for (j=0; j<vs->nc; j++) {
1395: PetscObjectReference((PetscObject)is_col[j]);
1397: vs->isglobal.col[j] = is_col[j];
1398: }
1399: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1400: offset = A->cmap->rstart;
1401: nsum = 0;
1402: for (j=0; j<vs->nc; j++) {
1403: MatNestFindNonzeroSubMatCol(A,j,&sub);
1404: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1405: MatGetLocalSize(sub,NULL,&n);
1406: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1407: nsum += n;
1408: }
1409: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1410: offset -= nsum;
1411: for (j=0; j<vs->nc; j++) {
1412: MatNestFindNonzeroSubMatCol(A,j,&sub);
1413: MatGetLocalSize(sub,NULL,&n);
1414: MatGetBlockSize(sub,&bs);
1415: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1416: ISSetBlockSize(vs->isglobal.col[j],bs);
1417: offset += n;
1418: }
1419: }
1421: /* Set up the local ISs */
1422: PetscMalloc1(vs->nr,&vs->islocal.row);
1423: PetscMalloc1(vs->nc,&vs->islocal.col);
1424: for (i=0,offset=0; i<vs->nr; i++) {
1425: IS isloc;
1426: ISLocalToGlobalMapping rmap = NULL;
1427: PetscInt nlocal,bs;
1428: MatNestFindNonzeroSubMatRow(A,i,&sub);
1429: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1430: if (rmap) {
1431: MatGetBlockSize(sub,&bs);
1432: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1433: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1434: ISSetBlockSize(isloc,bs);
1435: } else {
1436: nlocal = 0;
1437: isloc = NULL;
1438: }
1439: vs->islocal.row[i] = isloc;
1440: offset += nlocal;
1441: }
1442: for (i=0,offset=0; i<vs->nc; i++) {
1443: IS isloc;
1444: ISLocalToGlobalMapping cmap = NULL;
1445: PetscInt nlocal,bs;
1446: MatNestFindNonzeroSubMatCol(A,i,&sub);
1447: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1448: if (cmap) {
1449: MatGetBlockSize(sub,&bs);
1450: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1451: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1452: ISSetBlockSize(isloc,bs);
1453: } else {
1454: nlocal = 0;
1455: isloc = NULL;
1456: }
1457: vs->islocal.col[i] = isloc;
1458: offset += nlocal;
1459: }
1461: /* Set up the aggregate ISLocalToGlobalMapping */
1462: {
1463: ISLocalToGlobalMapping rmap,cmap;
1464: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1465: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1466: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1467: ISLocalToGlobalMappingDestroy(&rmap);
1468: ISLocalToGlobalMappingDestroy(&cmap);
1469: }
1471: #if defined(PETSC_USE_DEBUG)
1472: for (i=0; i<vs->nr; i++) {
1473: for (j=0; j<vs->nc; j++) {
1474: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1475: Mat B = vs->m[i][j];
1476: if (!B) continue;
1477: MatGetSize(B,&M,&N);
1478: MatGetLocalSize(B,&m,&n);
1479: ISGetSize(vs->isglobal.row[i],&Mi);
1480: ISGetSize(vs->isglobal.col[j],&Ni);
1481: ISGetLocalSize(vs->isglobal.row[i],&mi);
1482: ISGetLocalSize(vs->isglobal.col[j],&ni);
1483: 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);
1484: 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);
1485: }
1486: }
1487: #endif
1489: /* Set A->assembled if all non-null blocks are currently assembled */
1490: for (i=0; i<vs->nr; i++) {
1491: for (j=0; j<vs->nc; j++) {
1492: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1493: }
1494: }
1495: A->assembled = PETSC_TRUE;
1496: return(0);
1497: }
1499: /*@C
1500: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1502: Collective on Mat
1504: Input Parameter:
1505: + comm - Communicator for the new Mat
1506: . nr - number of nested row blocks
1507: . is_row - index sets for each nested row block, or NULL to make contiguous
1508: . nc - number of nested column blocks
1509: . is_col - index sets for each nested column block, or NULL to make contiguous
1510: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1512: Output Parameter:
1513: . B - new matrix
1515: Level: advanced
1517: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1518: @*/
1519: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1520: {
1521: Mat A;
1525: *B = 0;
1526: MatCreate(comm,&A);
1527: MatSetType(A,MATNEST);
1528: A->preallocated = PETSC_TRUE;
1529: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1530: *B = A;
1531: return(0);
1532: }
1534: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1535: {
1536: Mat_Nest *nest = (Mat_Nest*)A->data;
1537: Mat *trans;
1538: PetscScalar **avv;
1539: PetscScalar *vv;
1540: PetscInt **aii,**ajj;
1541: PetscInt *ii,*jj,*ci;
1542: PetscInt nr,nc,nnz,i,j;
1543: PetscBool done;
1547: MatGetSize(A,&nr,&nc);
1548: if (reuse == MAT_REUSE_MATRIX) {
1549: PetscInt rnr;
1551: MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1552: if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1553: if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1554: MatSeqAIJGetArray(*newmat,&vv);
1555: }
1556: /* extract CSR for nested SeqAIJ matrices */
1557: nnz = 0;
1558: PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1559: for (i=0; i<nest->nr; ++i) {
1560: for (j=0; j<nest->nc; ++j) {
1561: Mat B = nest->m[i][j];
1562: if (B) {
1563: PetscScalar *naa;
1564: PetscInt *nii,*njj,nnr;
1565: PetscBool istrans;
1567: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1568: if (istrans) {
1569: Mat Bt;
1571: MatTransposeGetMat(B,&Bt);
1572: MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1573: B = trans[i*nest->nc+j];
1574: }
1575: MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1576: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1577: MatSeqAIJGetArray(B,&naa);
1578: nnz += nii[nnr];
1580: aii[i*nest->nc+j] = nii;
1581: ajj[i*nest->nc+j] = njj;
1582: avv[i*nest->nc+j] = naa;
1583: }
1584: }
1585: }
1586: if (reuse != MAT_REUSE_MATRIX) {
1587: PetscMalloc1(nr+1,&ii);
1588: PetscMalloc1(nnz,&jj);
1589: PetscMalloc1(nnz,&vv);
1590: } else {
1591: if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1592: }
1594: /* new row pointer */
1595: PetscArrayzero(ii,nr+1);
1596: for (i=0; i<nest->nr; ++i) {
1597: PetscInt ncr,rst;
1599: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1600: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1601: for (j=0; j<nest->nc; ++j) {
1602: if (aii[i*nest->nc+j]) {
1603: PetscInt *nii = aii[i*nest->nc+j];
1604: PetscInt ir;
1606: for (ir=rst; ir<ncr+rst; ++ir) {
1607: ii[ir+1] += nii[1]-nii[0];
1608: nii++;
1609: }
1610: }
1611: }
1612: }
1613: for (i=0; i<nr; i++) ii[i+1] += ii[i];
1615: /* construct CSR for the new matrix */
1616: PetscCalloc1(nr,&ci);
1617: for (i=0; i<nest->nr; ++i) {
1618: PetscInt ncr,rst;
1620: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1621: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1622: for (j=0; j<nest->nc; ++j) {
1623: if (aii[i*nest->nc+j]) {
1624: PetscScalar *nvv = avv[i*nest->nc+j];
1625: PetscInt *nii = aii[i*nest->nc+j];
1626: PetscInt *njj = ajj[i*nest->nc+j];
1627: PetscInt ir,cst;
1629: ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1630: for (ir=rst; ir<ncr+rst; ++ir) {
1631: PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1633: for (ij=0;ij<rsize;ij++) {
1634: jj[ist+ij] = *njj+cst;
1635: vv[ist+ij] = *nvv;
1636: njj++;
1637: nvv++;
1638: }
1639: ci[ir] += rsize;
1640: nii++;
1641: }
1642: }
1643: }
1644: }
1645: PetscFree(ci);
1647: /* restore info */
1648: for (i=0; i<nest->nr; ++i) {
1649: for (j=0; j<nest->nc; ++j) {
1650: Mat B = nest->m[i][j];
1651: if (B) {
1652: PetscInt nnr = 0, k = i*nest->nc+j;
1654: B = (trans[k] ? trans[k] : B);
1655: MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1656: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1657: MatSeqAIJRestoreArray(B,&avv[k]);
1658: MatDestroy(&trans[k]);
1659: }
1660: }
1661: }
1662: PetscFree4(aii,ajj,avv,trans);
1664: /* finalize newmat */
1665: if (reuse == MAT_INITIAL_MATRIX) {
1666: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1667: } else if (reuse == MAT_INPLACE_MATRIX) {
1668: Mat B;
1670: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1671: MatHeaderReplace(A,&B);
1672: }
1673: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1674: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1675: {
1676: Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1677: a->free_a = PETSC_TRUE;
1678: a->free_ij = PETSC_TRUE;
1679: }
1680: return(0);
1681: }
1683: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1684: {
1686: Mat_Nest *nest = (Mat_Nest*)A->data;
1687: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1688: PetscInt cstart,cend;
1689: PetscMPIInt size;
1690: Mat C;
1693: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1694: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1695: PetscInt nf;
1696: PetscBool fast;
1698: PetscStrcmp(newtype,MATAIJ,&fast);
1699: if (!fast) {
1700: PetscStrcmp(newtype,MATSEQAIJ,&fast);
1701: }
1702: for (i=0; i<nest->nr && fast; ++i) {
1703: for (j=0; j<nest->nc && fast; ++j) {
1704: Mat B = nest->m[i][j];
1705: if (B) {
1706: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1707: if (!fast) {
1708: PetscBool istrans;
1710: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1711: if (istrans) {
1712: Mat Bt;
1714: MatTransposeGetMat(B,&Bt);
1715: PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1716: }
1717: }
1718: }
1719: }
1720: }
1721: for (i=0, nf=0; i<nest->nr && fast; ++i) {
1722: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1723: if (fast) {
1724: PetscInt f,s;
1726: ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1727: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1728: else {
1729: ISGetSize(nest->isglobal.row[i],&f);
1730: nf += f;
1731: }
1732: }
1733: }
1734: for (i=0, nf=0; i<nest->nc && fast; ++i) {
1735: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1736: if (fast) {
1737: PetscInt f,s;
1739: ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1740: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1741: else {
1742: ISGetSize(nest->isglobal.col[i],&f);
1743: nf += f;
1744: }
1745: }
1746: }
1747: if (fast) {
1748: MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1749: return(0);
1750: }
1751: }
1752: MatGetSize(A,&M,&N);
1753: MatGetLocalSize(A,&m,&n);
1754: MatGetOwnershipRangeColumn(A,&cstart,&cend);
1755: switch (reuse) {
1756: case MAT_INITIAL_MATRIX:
1757: MatCreate(PetscObjectComm((PetscObject)A),&C);
1758: MatSetType(C,newtype);
1759: MatSetSizes(C,m,n,M,N);
1760: *newmat = C;
1761: break;
1762: case MAT_REUSE_MATRIX:
1763: C = *newmat;
1764: break;
1765: default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1766: }
1767: PetscMalloc1(2*m,&dnnz);
1768: onnz = dnnz + m;
1769: for (k=0; k<m; k++) {
1770: dnnz[k] = 0;
1771: onnz[k] = 0;
1772: }
1773: for (j=0; j<nest->nc; ++j) {
1774: IS bNis;
1775: PetscInt bN;
1776: const PetscInt *bNindices;
1777: /* Using global column indices and ISAllGather() is not scalable. */
1778: ISAllGather(nest->isglobal.col[j], &bNis);
1779: ISGetSize(bNis, &bN);
1780: ISGetIndices(bNis,&bNindices);
1781: for (i=0; i<nest->nr; ++i) {
1782: PetscSF bmsf;
1783: PetscSFNode *iremote;
1784: Mat B;
1785: PetscInt bm, *sub_dnnz,*sub_onnz, br;
1786: const PetscInt *bmindices;
1787: B = nest->m[i][j];
1788: if (!B) continue;
1789: ISGetLocalSize(nest->isglobal.row[i],&bm);
1790: ISGetIndices(nest->isglobal.row[i],&bmindices);
1791: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1792: PetscMalloc1(bm,&iremote);
1793: PetscMalloc1(bm,&sub_dnnz);
1794: PetscMalloc1(bm,&sub_onnz);
1795: for (k = 0; k < bm; ++k){
1796: sub_dnnz[k] = 0;
1797: sub_onnz[k] = 0;
1798: }
1799: /*
1800: Locate the owners for all of the locally-owned global row indices for this row block.
1801: These determine the roots of PetscSF used to communicate preallocation data to row owners.
1802: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1803: */
1804: MatGetOwnershipRange(B,&rstart,NULL);
1805: for (br = 0; br < bm; ++br) {
1806: PetscInt row = bmindices[br], rowowner = 0, brncols, col;
1807: const PetscInt *brcols;
1808: PetscInt rowrel = 0; /* row's relative index on its owner rank */
1809: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1810: /* how many roots */
1811: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
1812: /* get nonzero pattern */
1813: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
1814: for (k=0; k<brncols; k++) {
1815: col = bNindices[brcols[k]];
1816: if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
1817: sub_dnnz[br]++;
1818: } else {
1819: sub_onnz[br]++;
1820: }
1821: }
1822: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
1823: }
1824: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1825: /* bsf will have to take care of disposing of bedges. */
1826: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1827: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1828: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1829: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1830: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1831: PetscFree(sub_dnnz);
1832: PetscFree(sub_onnz);
1833: PetscSFDestroy(&bmsf);
1834: }
1835: ISRestoreIndices(bNis,&bNindices);
1836: ISDestroy(&bNis);
1837: }
1838: /* Resize preallocation if overestimated */
1839: for (i=0;i<m;i++) {
1840: dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
1841: onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
1842: }
1843: MatSeqAIJSetPreallocation(C,0,dnnz);
1844: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1845: PetscFree(dnnz);
1847: /* Fill by row */
1848: for (j=0; j<nest->nc; ++j) {
1849: /* Using global column indices and ISAllGather() is not scalable. */
1850: IS bNis;
1851: PetscInt bN;
1852: const PetscInt *bNindices;
1853: ISAllGather(nest->isglobal.col[j], &bNis);
1854: ISGetSize(bNis,&bN);
1855: ISGetIndices(bNis,&bNindices);
1856: for (i=0; i<nest->nr; ++i) {
1857: Mat B;
1858: PetscInt bm, br;
1859: const PetscInt *bmindices;
1860: B = nest->m[i][j];
1861: if (!B) continue;
1862: ISGetLocalSize(nest->isglobal.row[i],&bm);
1863: ISGetIndices(nest->isglobal.row[i],&bmindices);
1864: MatGetOwnershipRange(B,&rstart,NULL);
1865: for (br = 0; br < bm; ++br) {
1866: PetscInt row = bmindices[br], brncols, *cols;
1867: const PetscInt *brcols;
1868: const PetscScalar *brcoldata;
1869: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1870: PetscMalloc1(brncols,&cols);
1871: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1872: /*
1873: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1874: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1875: */
1876: MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1877: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1878: PetscFree(cols);
1879: }
1880: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1881: }
1882: ISRestoreIndices(bNis,&bNindices);
1883: ISDestroy(&bNis);
1884: }
1885: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1886: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1887: return(0);
1888: }
1890: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
1891: {
1893: *has = PETSC_FALSE;
1894: if (op == MATOP_MULT_TRANSPOSE) {
1895: Mat_Nest *bA = (Mat_Nest*)mat->data;
1896: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1898: PetscBool flg;
1900: for (j=0; j<nc; j++) {
1901: for (i=0; i<nr; i++) {
1902: if (!bA->m[i][j]) continue;
1903: MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);
1904: if (!flg) return(0);
1905: }
1906: }
1907: }
1908: if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
1909: return(0);
1910: }
1912: /*MC
1913: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
1915: Level: intermediate
1917: Notes:
1918: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
1919: It allows the use of symmetric and block formats for parts of multi-physics simulations.
1920: It is usually used with DMComposite and DMCreateMatrix()
1922: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
1923: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
1924: than the nest matrix.
1926: .seealso: MatCreate(), MatType, MatCreateNest()
1927: M*/
1928: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1929: {
1930: Mat_Nest *s;
1934: PetscNewLog(A,&s);
1935: A->data = (void*)s;
1937: s->nr = -1;
1938: s->nc = -1;
1939: s->m = NULL;
1940: s->splitassembly = PETSC_FALSE;
1942: PetscMemzero(A->ops,sizeof(*A->ops));
1944: A->ops->mult = MatMult_Nest;
1945: A->ops->multadd = MatMultAdd_Nest;
1946: A->ops->multtranspose = MatMultTranspose_Nest;
1947: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
1948: A->ops->transpose = MatTranspose_Nest;
1949: A->ops->assemblybegin = MatAssemblyBegin_Nest;
1950: A->ops->assemblyend = MatAssemblyEnd_Nest;
1951: A->ops->zeroentries = MatZeroEntries_Nest;
1952: A->ops->copy = MatCopy_Nest;
1953: A->ops->axpy = MatAXPY_Nest;
1954: A->ops->duplicate = MatDuplicate_Nest;
1955: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
1956: A->ops->destroy = MatDestroy_Nest;
1957: A->ops->view = MatView_Nest;
1958: A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1959: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
1960: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1961: A->ops->getdiagonal = MatGetDiagonal_Nest;
1962: A->ops->diagonalscale = MatDiagonalScale_Nest;
1963: A->ops->scale = MatScale_Nest;
1964: A->ops->shift = MatShift_Nest;
1965: A->ops->diagonalset = MatDiagonalSet_Nest;
1966: A->ops->setrandom = MatSetRandom_Nest;
1967: A->ops->hasoperation = MatHasOperation_Nest;
1969: A->spptr = 0;
1970: A->assembled = PETSC_FALSE;
1972: /* expose Nest api's */
1973: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
1974: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
1975: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
1976: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
1977: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
1978: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
1979: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
1980: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
1981: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",MatConvert_Nest_AIJ);
1982: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",MatConvert_Nest_AIJ);
1983: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
1984: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);
1986: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1987: return(0);
1988: }