Actual source code: matnest.c
petsc-3.14.6 2021-03-30
1: #include <../src/mat/impls/nest/matnestimpl.h>
2: #include <../src/mat/impls/aij/seq/aij.h>
3: #include <petscsf.h>
5: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
6: static PetscErrorCode MatCreateVecs_Nest(Mat,Vec*,Vec*);
7: static PetscErrorCode MatReset_Nest(Mat);
9: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*);
11: /* private functions */
12: static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
13: {
14: Mat_Nest *bA = (Mat_Nest*)A->data;
15: PetscInt i,j;
19: *m = *n = *M = *N = 0;
20: for (i=0; i<bA->nr; i++) { /* rows */
21: PetscInt sm,sM;
22: ISGetLocalSize(bA->isglobal.row[i],&sm);
23: ISGetSize(bA->isglobal.row[i],&sM);
24: *m += sm;
25: *M += sM;
26: }
27: for (j=0; j<bA->nc; j++) { /* cols */
28: PetscInt sn,sN;
29: ISGetLocalSize(bA->isglobal.col[j],&sn);
30: ISGetSize(bA->isglobal.col[j],&sN);
31: *n += sn;
32: *N += sN;
33: }
34: return(0);
35: }
37: /* 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: }
61: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
62: {
63: Mat_Nest *bA = (Mat_Nest*)A->data;
64: Vec *bx = bA->right,*bz = bA->left;
65: PetscInt i,j,nr = bA->nr,nc = bA->nc;
69: for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
70: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
71: for (i=0; i<nr; i++) {
72: if (y != z) {
73: Vec by;
74: VecGetSubVector(y,bA->isglobal.row[i],&by);
75: VecCopy(by,bz[i]);
76: VecRestoreSubVector(y,bA->isglobal.row[i],&by);
77: }
78: for (j=0; j<nc; j++) {
79: if (!bA->m[i][j]) continue;
80: /* y[i] <- y[i] + A[i][j] * x[j] */
81: MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
82: }
83: }
84: for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
85: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
86: return(0);
87: }
89: typedef struct {
90: Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
91: PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */
92: PetscInt *dm,*dn,k; /* displacements and number of submatrices */
93: } Nest_Dense;
95: PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
96: {
97: Mat_Nest *bA;
98: Nest_Dense *contents;
99: Mat viewB,viewC,productB,workC;
100: const PetscScalar *barray;
101: PetscScalar *carray;
102: PetscInt i,j,M,N,nr,nc,ldb,ldc;
103: PetscErrorCode ierr;
104: Mat A,B;
107: MatCheckProduct(C,3);
108: A = C->product->A;
109: B = C->product->B;
110: MatGetSize(B,NULL,&N);
111: if (!N) {
112: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
114: return(0);
115: }
116: contents = (Nest_Dense*)C->product->data;
117: if (!contents) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
118: bA = (Mat_Nest*)A->data;
119: nr = bA->nr;
120: nc = bA->nc;
121: MatDenseGetLDA(B,&ldb);
122: MatDenseGetLDA(C,&ldc);
123: MatZeroEntries(C);
124: MatDenseGetArrayRead(B,&barray);
125: MatDenseGetArrayWrite(C,&carray);
126: for (i=0; i<nr; i++) {
127: ISGetSize(bA->isglobal.row[i],&M);
128: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
129: MatDenseSetLDA(viewC,ldc);
130: for (j=0; j<nc; j++) {
131: if (!bA->m[i][j]) continue;
132: ISGetSize(bA->isglobal.col[j],&M);
133: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
134: MatDenseSetLDA(viewB,ldb);
136: /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
137: workC = contents->workC[i*nc + j];
138: productB = workC->product->B;
139: workC->product->B = viewB; /* use newly created dense matrix viewB */
140: MatProductNumeric(workC);
141: MatDestroy(&viewB);
142: workC->product->B = productB; /* resume original B */
144: /* C[i] <- workC + C[i] */
145: MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
146: }
147: MatDestroy(&viewC);
148: }
149: MatDenseRestoreArrayWrite(C,&carray);
150: MatDenseRestoreArrayRead(B,&barray);
152: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
154: return(0);
155: }
157: PetscErrorCode MatNest_DenseDestroy(void *ctx)
158: {
159: Nest_Dense *contents = (Nest_Dense*)ctx;
160: PetscInt i;
164: PetscFree(contents->tarray);
165: for (i=0; i<contents->k; i++) {
166: MatDestroy(contents->workC + i);
167: }
168: PetscFree3(contents->dm,contents->dn,contents->workC);
169: PetscFree(contents);
170: return(0);
171: }
173: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
174: {
175: Mat_Nest *bA;
176: Mat viewB,workC;
177: const PetscScalar *barray;
178: PetscInt i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
179: Nest_Dense *contents=NULL;
180: PetscBool cisdense;
181: PetscErrorCode ierr;
182: Mat A,B;
183: PetscReal fill;
186: MatCheckProduct(C,4);
187: if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
188: A = C->product->A;
189: B = C->product->B;
190: fill = C->product->fill;
191: bA = (Mat_Nest*)A->data;
192: nr = bA->nr;
193: nc = bA->nc;
194: MatGetLocalSize(B,NULL,&n);
195: MatGetSize(B,NULL,&N);
196: MatGetLocalSize(A,&m,NULL);
197: MatGetSize(A,&M,NULL);
198: MatSetSizes(C,m,n,M,N);
199: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"");
200: if (!cisdense) {
201: MatSetType(C,((PetscObject)B)->type_name);
202: }
203: MatSetUp(C);
204: if (!N) {
205: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
206: return(0);
207: }
209: PetscNew(&contents);
210: C->product->data = contents;
211: C->product->destroy = MatNest_DenseDestroy;
212: PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
213: contents->k = nr*nc;
214: for (i=0; i<nr; i++) {
215: ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
216: maxm = PetscMax(maxm,contents->dm[i+1]);
217: contents->dm[i+1] += contents->dm[i];
218: }
219: for (i=0; i<nc; i++) {
220: ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
221: contents->dn[i+1] += contents->dn[i];
222: }
223: PetscMalloc1(maxm*N,&contents->tarray);
224: MatDenseGetLDA(B,&ldb);
225: MatGetSize(B,NULL,&N);
226: MatDenseGetArrayRead(B,&barray);
227: /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
228: for (j=0; j<nc; j++) {
229: ISGetSize(bA->isglobal.col[j],&M);
230: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
231: MatDenseSetLDA(viewB,ldb);
232: for (i=0; i<nr; i++) {
233: if (!bA->m[i][j]) continue;
234: /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
236: MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
237: workC = contents->workC[i*nc + j];
238: MatProductSetType(workC,MATPRODUCT_AB);
239: MatProductSetAlgorithm(workC,"default");
240: MatProductSetFill(workC,fill);
241: MatProductSetFromOptions(workC);
242: MatProductSymbolic(workC);
244: /* since tarray will be shared by all Mat */
245: MatSeqDenseSetPreallocation(workC,contents->tarray);
246: MatMPIDenseSetPreallocation(workC,contents->tarray);
247: }
248: MatDestroy(&viewB);
249: }
250: MatDenseRestoreArrayRead(B,&barray);
252: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
253: return(0);
254: }
256: /* --------------------------------------------------------- */
257: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
258: {
260: C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
261: return(0);
262: }
264: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
265: {
267: Mat_Product *product = C->product;
270: if (product->type == MATPRODUCT_AB) {
271: MatProductSetFromOptions_Nest_Dense_AB(C);
272: }
273: return(0);
274: }
275: /* --------------------------------------------------------- */
277: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
278: {
279: Mat_Nest *bA = (Mat_Nest*)A->data;
280: Vec *bx = bA->left,*by = bA->right;
281: PetscInt i,j,nr = bA->nr,nc = bA->nc;
285: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
286: for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
287: for (j=0; j<nc; j++) {
288: VecZeroEntries(by[j]);
289: for (i=0; i<nr; i++) {
290: if (!bA->m[i][j]) continue;
291: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
292: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
293: }
294: }
295: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
296: for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
297: return(0);
298: }
300: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
301: {
302: Mat_Nest *bA = (Mat_Nest*)A->data;
303: Vec *bx = bA->left,*bz = bA->right;
304: PetscInt i,j,nr = bA->nr,nc = bA->nc;
308: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
309: for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
310: for (j=0; j<nc; j++) {
311: if (y != z) {
312: Vec by;
313: VecGetSubVector(y,bA->isglobal.col[j],&by);
314: VecCopy(by,bz[j]);
315: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
316: }
317: for (i=0; i<nr; i++) {
318: if (!bA->m[i][j]) continue;
319: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
320: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
321: }
322: }
323: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
324: for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
325: return(0);
326: }
328: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
329: {
330: Mat_Nest *bA = (Mat_Nest*)A->data, *bC;
331: Mat C;
332: PetscInt i,j,nr = bA->nr,nc = bA->nc;
336: if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");
338: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
339: Mat *subs;
340: IS *is_row,*is_col;
342: PetscCalloc1(nr * nc,&subs);
343: PetscMalloc2(nr,&is_row,nc,&is_col);
344: MatNestGetISs(A,is_row,is_col);
345: if (reuse == MAT_INPLACE_MATRIX) {
346: for (i=0; i<nr; i++) {
347: for (j=0; j<nc; j++) {
348: subs[i + nr * j] = bA->m[i][j];
349: }
350: }
351: }
353: MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
354: PetscFree(subs);
355: PetscFree2(is_row,is_col);
356: } else {
357: C = *B;
358: }
360: bC = (Mat_Nest*)C->data;
361: for (i=0; i<nr; i++) {
362: for (j=0; j<nc; j++) {
363: if (bA->m[i][j]) {
364: MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
365: } else {
366: bC->m[j][i] = NULL;
367: }
368: }
369: }
371: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
372: *B = C;
373: } else {
374: MatHeaderMerge(A, &C);
375: }
376: return(0);
377: }
379: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
380: {
382: IS *lst = *list;
383: PetscInt i;
386: if (!lst) return(0);
387: for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
388: PetscFree(lst);
389: *list = NULL;
390: return(0);
391: }
393: static PetscErrorCode MatReset_Nest(Mat A)
394: {
395: Mat_Nest *vs = (Mat_Nest*)A->data;
396: PetscInt i,j;
400: /* release the matrices and the place holders */
401: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
402: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
403: MatNestDestroyISList(vs->nr,&vs->islocal.row);
404: MatNestDestroyISList(vs->nc,&vs->islocal.col);
406: PetscFree(vs->row_len);
407: PetscFree(vs->col_len);
408: PetscFree(vs->nnzstate);
410: PetscFree2(vs->left,vs->right);
412: /* release the matrices and the place holders */
413: if (vs->m) {
414: for (i=0; i<vs->nr; i++) {
415: for (j=0; j<vs->nc; j++) {
416: MatDestroy(&vs->m[i][j]);
417: }
418: PetscFree(vs->m[i]);
419: }
420: PetscFree(vs->m);
421: }
423: /* restore defaults */
424: vs->nr = 0;
425: vs->nc = 0;
426: vs->splitassembly = PETSC_FALSE;
427: return(0);
428: }
430: static PetscErrorCode MatDestroy_Nest(Mat A)
431: {
434: MatReset_Nest(A);
435: PetscFree(A->data);
437: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
438: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
439: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
440: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
441: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
442: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
443: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
444: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
445: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
446: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
447: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
448: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
449: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
450: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
451: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
452: return(0);
453: }
455: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
456: {
457: Mat_Nest *vs = (Mat_Nest*)mat->data;
458: PetscInt i;
462: if (dd) *dd = 0;
463: if (!vs->nr) {
464: *missing = PETSC_TRUE;
465: return(0);
466: }
467: *missing = PETSC_FALSE;
468: for (i = 0; i < vs->nr && !(*missing); i++) {
469: *missing = PETSC_TRUE;
470: if (vs->m[i][i]) {
471: MatMissingDiagonal(vs->m[i][i],missing,NULL);
472: if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
473: }
474: }
475: return(0);
476: }
478: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
479: {
480: Mat_Nest *vs = (Mat_Nest*)A->data;
481: PetscInt i,j;
483: PetscBool nnzstate = PETSC_FALSE;
486: for (i=0; i<vs->nr; i++) {
487: for (j=0; j<vs->nc; j++) {
488: PetscObjectState subnnzstate = 0;
489: if (vs->m[i][j]) {
490: MatAssemblyBegin(vs->m[i][j],type);
491: if (!vs->splitassembly) {
492: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
493: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
494: * already performing an assembly, but the result would by more complicated and appears to offer less
495: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
496: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
497: */
498: MatAssemblyEnd(vs->m[i][j],type);
499: MatGetNonzeroState(vs->m[i][j],&subnnzstate);
500: }
501: }
502: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
503: vs->nnzstate[i*vs->nc+j] = subnnzstate;
504: }
505: }
506: if (nnzstate) A->nonzerostate++;
507: return(0);
508: }
510: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
511: {
512: Mat_Nest *vs = (Mat_Nest*)A->data;
513: PetscInt i,j;
517: for (i=0; i<vs->nr; i++) {
518: for (j=0; j<vs->nc; j++) {
519: if (vs->m[i][j]) {
520: if (vs->splitassembly) {
521: MatAssemblyEnd(vs->m[i][j],type);
522: }
523: }
524: }
525: }
526: return(0);
527: }
529: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
530: {
532: Mat_Nest *vs = (Mat_Nest*)A->data;
533: PetscInt j;
534: Mat sub;
537: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
538: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
539: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
540: *B = sub;
541: return(0);
542: }
544: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
545: {
547: Mat_Nest *vs = (Mat_Nest*)A->data;
548: PetscInt i;
549: Mat sub;
552: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
553: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
554: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
555: *B = sub;
556: return(0);
557: }
559: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
560: {
562: PetscInt i;
563: PetscBool flg;
569: *found = -1;
570: for (i=0; i<n; i++) {
571: if (!list[i]) continue;
572: ISEqualUnsorted(list[i],is,&flg);
573: if (flg) {
574: *found = i;
575: return(0);
576: }
577: }
578: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
579: }
581: /* Get a block row as a new MatNest */
582: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
583: {
584: Mat_Nest *vs = (Mat_Nest*)A->data;
585: char keyname[256];
589: *B = NULL;
590: PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
591: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
592: if (*B) return(0);
594: MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);
596: (*B)->assembled = A->assembled;
598: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
599: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
600: return(0);
601: }
603: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
604: {
605: Mat_Nest *vs = (Mat_Nest*)A->data;
607: PetscInt row,col;
608: PetscBool same,isFullCol,isFullColGlobal;
611: /* Check if full column space. This is a hack */
612: isFullCol = PETSC_FALSE;
613: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
614: if (same) {
615: PetscInt n,first,step,i,an,am,afirst,astep;
616: ISStrideGetInfo(iscol,&first,&step);
617: ISGetLocalSize(iscol,&n);
618: isFullCol = PETSC_TRUE;
619: for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
620: PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
621: ISGetLocalSize(is->col[i],&am);
622: if (same) {
623: ISStrideGetInfo(is->col[i],&afirst,&astep);
624: if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
625: } else isFullCol = PETSC_FALSE;
626: an += am;
627: }
628: if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
629: }
630: MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));
632: if (isFullColGlobal && vs->nc > 1) {
633: PetscInt row;
634: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
635: MatNestGetRow(A,row,B);
636: } else {
637: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
638: MatNestFindIS(A,vs->nc,is->col,iscol,&col);
639: if (!vs->m[row][col]) {
640: PetscInt lr,lc;
642: MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
643: ISGetLocalSize(vs->isglobal.row[row],&lr);
644: ISGetLocalSize(vs->isglobal.col[col],&lc);
645: MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
646: MatSetType(vs->m[row][col],MATAIJ);
647: MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
648: MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
649: MatSetUp(vs->m[row][col]);
650: MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
651: MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
652: }
653: *B = vs->m[row][col];
654: }
655: return(0);
656: }
658: /*
659: TODO: This does not actually returns a submatrix we can modify
660: */
661: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
662: {
664: Mat_Nest *vs = (Mat_Nest*)A->data;
665: Mat sub;
668: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
669: switch (reuse) {
670: case MAT_INITIAL_MATRIX:
671: if (sub) { PetscObjectReference((PetscObject)sub); }
672: *B = sub;
673: break;
674: case MAT_REUSE_MATRIX:
675: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
676: break;
677: case MAT_IGNORE_MATRIX: /* Nothing to do */
678: break;
679: case MAT_INPLACE_MATRIX: /* Nothing to do */
680: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
681: }
682: return(0);
683: }
685: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
686: {
688: Mat_Nest *vs = (Mat_Nest*)A->data;
689: Mat sub;
692: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
693: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
694: if (sub) {PetscObjectReference((PetscObject)sub);}
695: *B = sub;
696: return(0);
697: }
699: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
700: {
702: Mat_Nest *vs = (Mat_Nest*)A->data;
703: Mat sub;
706: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
707: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
708: if (sub) {
709: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
710: MatDestroy(B);
711: }
712: return(0);
713: }
715: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
716: {
717: Mat_Nest *bA = (Mat_Nest*)A->data;
718: PetscInt i;
722: for (i=0; i<bA->nr; i++) {
723: Vec bv;
724: VecGetSubVector(v,bA->isglobal.row[i],&bv);
725: if (bA->m[i][i]) {
726: MatGetDiagonal(bA->m[i][i],bv);
727: } else {
728: VecSet(bv,0.0);
729: }
730: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
731: }
732: return(0);
733: }
735: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
736: {
737: Mat_Nest *bA = (Mat_Nest*)A->data;
738: Vec bl,*br;
739: PetscInt i,j;
743: PetscCalloc1(bA->nc,&br);
744: if (r) {
745: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
746: }
747: bl = NULL;
748: for (i=0; i<bA->nr; i++) {
749: if (l) {
750: VecGetSubVector(l,bA->isglobal.row[i],&bl);
751: }
752: for (j=0; j<bA->nc; j++) {
753: if (bA->m[i][j]) {
754: MatDiagonalScale(bA->m[i][j],bl,br[j]);
755: }
756: }
757: if (l) {
758: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
759: }
760: }
761: if (r) {
762: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
763: }
764: PetscFree(br);
765: return(0);
766: }
768: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
769: {
770: Mat_Nest *bA = (Mat_Nest*)A->data;
771: PetscInt i,j;
775: for (i=0; i<bA->nr; i++) {
776: for (j=0; j<bA->nc; j++) {
777: if (bA->m[i][j]) {
778: MatScale(bA->m[i][j],a);
779: }
780: }
781: }
782: return(0);
783: }
785: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
786: {
787: Mat_Nest *bA = (Mat_Nest*)A->data;
788: PetscInt i;
790: PetscBool nnzstate = PETSC_FALSE;
793: for (i=0; i<bA->nr; i++) {
794: PetscObjectState subnnzstate = 0;
795: 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);
796: MatShift(bA->m[i][i],a);
797: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
798: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
799: bA->nnzstate[i*bA->nc+i] = subnnzstate;
800: }
801: if (nnzstate) A->nonzerostate++;
802: return(0);
803: }
805: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
806: {
807: Mat_Nest *bA = (Mat_Nest*)A->data;
808: PetscInt i;
810: PetscBool nnzstate = PETSC_FALSE;
813: for (i=0; i<bA->nr; i++) {
814: PetscObjectState subnnzstate = 0;
815: Vec bv;
816: VecGetSubVector(D,bA->isglobal.row[i],&bv);
817: if (bA->m[i][i]) {
818: MatDiagonalSet(bA->m[i][i],bv,is);
819: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
820: }
821: VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
822: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
823: bA->nnzstate[i*bA->nc+i] = subnnzstate;
824: }
825: if (nnzstate) A->nonzerostate++;
826: return(0);
827: }
829: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
830: {
831: Mat_Nest *bA = (Mat_Nest*)A->data;
832: PetscInt i,j;
836: for (i=0; i<bA->nr; i++) {
837: for (j=0; j<bA->nc; j++) {
838: if (bA->m[i][j]) {
839: MatSetRandom(bA->m[i][j],rctx);
840: }
841: }
842: }
843: return(0);
844: }
846: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
847: {
848: Mat_Nest *bA = (Mat_Nest*)A->data;
849: Vec *L,*R;
850: MPI_Comm comm;
851: PetscInt i,j;
855: PetscObjectGetComm((PetscObject)A,&comm);
856: if (right) {
857: /* allocate R */
858: PetscMalloc1(bA->nc, &R);
859: /* Create the right vectors */
860: for (j=0; j<bA->nc; j++) {
861: for (i=0; i<bA->nr; i++) {
862: if (bA->m[i][j]) {
863: MatCreateVecs(bA->m[i][j],&R[j],NULL);
864: break;
865: }
866: }
867: if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
868: }
869: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
870: /* hand back control to the nest vector */
871: for (j=0; j<bA->nc; j++) {
872: VecDestroy(&R[j]);
873: }
874: PetscFree(R);
875: }
877: if (left) {
878: /* allocate L */
879: PetscMalloc1(bA->nr, &L);
880: /* Create the left vectors */
881: for (i=0; i<bA->nr; i++) {
882: for (j=0; j<bA->nc; j++) {
883: if (bA->m[i][j]) {
884: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
885: break;
886: }
887: }
888: if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
889: }
891: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
892: for (i=0; i<bA->nr; i++) {
893: VecDestroy(&L[i]);
894: }
896: PetscFree(L);
897: }
898: return(0);
899: }
901: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
902: {
903: Mat_Nest *bA = (Mat_Nest*)A->data;
904: PetscBool isascii,viewSub = PETSC_FALSE;
905: PetscInt i,j;
909: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
910: if (isascii) {
912: PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
913: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
914: PetscViewerASCIIPushTab(viewer);
915: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);
917: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
918: for (i=0; i<bA->nr; i++) {
919: for (j=0; j<bA->nc; j++) {
920: MatType type;
921: char name[256] = "",prefix[256] = "";
922: PetscInt NR,NC;
923: PetscBool isNest = PETSC_FALSE;
925: if (!bA->m[i][j]) {
926: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
927: continue;
928: }
929: MatGetSize(bA->m[i][j],&NR,&NC);
930: MatGetType(bA->m[i][j], &type);
931: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
932: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
933: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
935: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
937: if (isNest || viewSub) {
938: PetscViewerASCIIPushTab(viewer); /* push1 */
939: MatView(bA->m[i][j],viewer);
940: PetscViewerASCIIPopTab(viewer); /* pop1 */
941: }
942: }
943: }
944: PetscViewerASCIIPopTab(viewer); /* pop0 */
945: }
946: return(0);
947: }
949: static PetscErrorCode MatZeroEntries_Nest(Mat A)
950: {
951: Mat_Nest *bA = (Mat_Nest*)A->data;
952: PetscInt i,j;
956: for (i=0; i<bA->nr; i++) {
957: for (j=0; j<bA->nc; j++) {
958: if (!bA->m[i][j]) continue;
959: MatZeroEntries(bA->m[i][j]);
960: }
961: }
962: return(0);
963: }
965: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
966: {
967: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
968: PetscInt i,j,nr = bA->nr,nc = bA->nc;
970: PetscBool nnzstate = PETSC_FALSE;
973: 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);
974: for (i=0; i<nr; i++) {
975: for (j=0; j<nc; j++) {
976: PetscObjectState subnnzstate = 0;
977: if (bA->m[i][j] && bB->m[i][j]) {
978: MatCopy(bA->m[i][j],bB->m[i][j],str);
979: } 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);
980: MatGetNonzeroState(bB->m[i][j],&subnnzstate);
981: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
982: bB->nnzstate[i*nc+j] = subnnzstate;
983: }
984: }
985: if (nnzstate) B->nonzerostate++;
986: return(0);
987: }
989: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
990: {
991: Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
992: PetscInt i,j,nr = bY->nr,nc = bY->nc;
994: PetscBool nnzstate = PETSC_FALSE;
997: 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);
998: for (i=0; i<nr; i++) {
999: for (j=0; j<nc; j++) {
1000: PetscObjectState subnnzstate = 0;
1001: if (bY->m[i][j] && bX->m[i][j]) {
1002: MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
1003: } else if (bX->m[i][j]) {
1004: Mat M;
1006: 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);
1007: MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
1008: MatNestSetSubMat(Y,i,j,M);
1009: MatDestroy(&M);
1010: }
1011: if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
1012: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1013: bY->nnzstate[i*nc+j] = subnnzstate;
1014: }
1015: }
1016: if (nnzstate) Y->nonzerostate++;
1017: return(0);
1018: }
1020: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1021: {
1022: Mat_Nest *bA = (Mat_Nest*)A->data;
1023: Mat *b;
1024: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1028: PetscMalloc1(nr*nc,&b);
1029: for (i=0; i<nr; i++) {
1030: for (j=0; j<nc; j++) {
1031: if (bA->m[i][j]) {
1032: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1033: } else {
1034: b[i*nc+j] = NULL;
1035: }
1036: }
1037: }
1038: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1039: /* Give the new MatNest exclusive ownership */
1040: for (i=0; i<nr*nc; i++) {
1041: MatDestroy(&b[i]);
1042: }
1043: PetscFree(b);
1045: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1046: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1047: return(0);
1048: }
1050: /* nest api */
1051: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1052: {
1053: Mat_Nest *bA = (Mat_Nest*)A->data;
1056: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1057: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1058: *mat = bA->m[idxm][jdxm];
1059: return(0);
1060: }
1062: /*@
1063: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
1065: Not collective
1067: Input Parameters:
1068: + A - nest matrix
1069: . idxm - index of the matrix within the nest matrix
1070: - jdxm - index of the matrix within the nest matrix
1072: Output Parameter:
1073: . sub - matrix at index idxm,jdxm within the nest matrix
1075: Level: developer
1077: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1078: MatNestGetLocalISs(), MatNestGetISs()
1079: @*/
1080: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1081: {
1085: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1086: return(0);
1087: }
1089: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1090: {
1091: Mat_Nest *bA = (Mat_Nest*)A->data;
1092: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1096: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1097: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1098: MatGetLocalSize(mat,&m,&n);
1099: MatGetSize(mat,&M,&N);
1100: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1101: ISGetSize(bA->isglobal.row[idxm],&Mi);
1102: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1103: ISGetSize(bA->isglobal.col[jdxm],&Ni);
1104: 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);
1105: 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);
1107: /* do not increase object state */
1108: if (mat == bA->m[idxm][jdxm]) return(0);
1110: PetscObjectReference((PetscObject)mat);
1111: MatDestroy(&bA->m[idxm][jdxm]);
1112: bA->m[idxm][jdxm] = mat;
1113: PetscObjectStateIncrease((PetscObject)A);
1114: MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1115: A->nonzerostate++;
1116: return(0);
1117: }
1119: /*@
1120: MatNestSetSubMat - Set a single submatrix in the nest matrix.
1122: Logically collective on the submatrix communicator
1124: Input Parameters:
1125: + A - nest matrix
1126: . idxm - index of the matrix within the nest matrix
1127: . jdxm - index of the matrix within the nest matrix
1128: - sub - matrix at index idxm,jdxm within the nest matrix
1130: Notes:
1131: The new submatrix must have the same size and communicator as that block of the nest.
1133: This increments the reference count of the submatrix.
1135: Level: developer
1137: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1138: MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1139: @*/
1140: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1141: {
1145: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1146: return(0);
1147: }
1149: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1150: {
1151: Mat_Nest *bA = (Mat_Nest*)A->data;
1154: if (M) *M = bA->nr;
1155: if (N) *N = bA->nc;
1156: if (mat) *mat = bA->m;
1157: return(0);
1158: }
1160: /*@C
1161: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
1163: Not collective
1165: Input Parameters:
1166: . A - nest matrix
1168: Output Parameter:
1169: + M - number of rows in the nest matrix
1170: . N - number of cols in the nest matrix
1171: - mat - 2d array of matrices
1173: Notes:
1175: The user should not free the array mat.
1177: In Fortran, this routine has a calling sequence
1178: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1179: where the space allocated for the optional argument mat is assumed large enough (if provided).
1181: Level: developer
1183: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1184: MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1185: @*/
1186: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1187: {
1191: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1192: return(0);
1193: }
1195: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1196: {
1197: Mat_Nest *bA = (Mat_Nest*)A->data;
1200: if (M) *M = bA->nr;
1201: if (N) *N = bA->nc;
1202: return(0);
1203: }
1205: /*@
1206: MatNestGetSize - Returns the size of the nest matrix.
1208: Not collective
1210: Input Parameters:
1211: . A - nest matrix
1213: Output Parameter:
1214: + M - number of rows in the nested mat
1215: - N - number of cols in the nested mat
1217: Notes:
1219: Level: developer
1221: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1222: MatNestGetISs()
1223: @*/
1224: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1225: {
1229: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1230: return(0);
1231: }
1233: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1234: {
1235: Mat_Nest *vs = (Mat_Nest*)A->data;
1236: PetscInt i;
1239: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1240: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1241: return(0);
1242: }
1244: /*@C
1245: MatNestGetISs - Returns the index sets partitioning the row and column spaces
1247: Not collective
1249: Input Parameters:
1250: . A - nest matrix
1252: Output Parameter:
1253: + rows - array of row index sets
1254: - cols - array of column index sets
1256: Level: advanced
1258: Notes:
1259: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1261: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1262: MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1263: @*/
1264: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1265: {
1270: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1271: return(0);
1272: }
1274: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1275: {
1276: Mat_Nest *vs = (Mat_Nest*)A->data;
1277: PetscInt i;
1280: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1281: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1282: return(0);
1283: }
1285: /*@C
1286: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1288: Not collective
1290: Input Parameters:
1291: . A - nest matrix
1293: Output Parameter:
1294: + rows - array of row index sets (or NULL to ignore)
1295: - cols - array of column index sets (or NULL to ignore)
1297: Level: advanced
1299: Notes:
1300: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1302: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1303: MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1304: @*/
1305: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1306: {
1311: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1312: return(0);
1313: }
1315: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1316: {
1318: PetscBool flg;
1321: PetscStrcmp(vtype,VECNEST,&flg);
1322: /* In reality, this only distinguishes VECNEST and "other" */
1323: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1324: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1325: return(0);
1326: }
1328: /*@C
1329: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1331: Not collective
1333: Input Parameters:
1334: + A - nest matrix
1335: - vtype - type to use for creating vectors
1337: Notes:
1339: Level: developer
1341: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1342: @*/
1343: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1344: {
1348: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1349: return(0);
1350: }
1352: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1353: {
1354: Mat_Nest *s = (Mat_Nest*)A->data;
1355: PetscInt i,j,m,n,M,N;
1357: PetscBool cong;
1360: MatReset_Nest(A);
1362: s->nr = nr;
1363: s->nc = nc;
1365: /* Create space for submatrices */
1366: PetscMalloc1(nr,&s->m);
1367: for (i=0; i<nr; i++) {
1368: PetscMalloc1(nc,&s->m[i]);
1369: }
1370: for (i=0; i<nr; i++) {
1371: for (j=0; j<nc; j++) {
1372: s->m[i][j] = a[i*nc+j];
1373: if (a[i*nc+j]) {
1374: PetscObjectReference((PetscObject)a[i*nc+j]);
1375: }
1376: }
1377: }
1379: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1381: PetscMalloc1(nr,&s->row_len);
1382: PetscMalloc1(nc,&s->col_len);
1383: for (i=0; i<nr; i++) s->row_len[i]=-1;
1384: for (j=0; j<nc; j++) s->col_len[j]=-1;
1386: PetscCalloc1(nr*nc,&s->nnzstate);
1387: for (i=0; i<nr; i++) {
1388: for (j=0; j<nc; j++) {
1389: if (s->m[i][j]) {
1390: MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1391: }
1392: }
1393: }
1395: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1397: PetscLayoutSetSize(A->rmap,M);
1398: PetscLayoutSetLocalSize(A->rmap,m);
1399: PetscLayoutSetSize(A->cmap,N);
1400: PetscLayoutSetLocalSize(A->cmap,n);
1402: PetscLayoutSetUp(A->rmap);
1403: PetscLayoutSetUp(A->cmap);
1405: /* disable operations that are not supported for non-square matrices,
1406: or matrices for which is_row != is_col */
1407: MatHasCongruentLayouts(A,&cong);
1408: if (cong && nr != nc) cong = PETSC_FALSE;
1409: if (cong) {
1410: for (i = 0; cong && i < nr; i++) {
1411: ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1412: }
1413: }
1414: if (!cong) {
1415: A->ops->missingdiagonal = NULL;
1416: A->ops->getdiagonal = NULL;
1417: A->ops->shift = NULL;
1418: A->ops->diagonalset = NULL;
1419: }
1421: PetscCalloc2(nr,&s->left,nc,&s->right);
1422: PetscObjectStateIncrease((PetscObject)A);
1423: A->nonzerostate++;
1424: return(0);
1425: }
1427: /*@
1428: MatNestSetSubMats - Sets the nested submatrices
1430: Collective on Mat
1432: Input Parameter:
1433: + A - nested matrix
1434: . nr - number of nested row blocks
1435: . is_row - index sets for each nested row block, or NULL to make contiguous
1436: . nc - number of nested column blocks
1437: . is_col - index sets for each nested column block, or NULL to make contiguous
1438: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1440: Notes: this always resets any submatrix information previously set
1442: Level: advanced
1444: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1445: @*/
1446: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1447: {
1449: PetscInt i;
1453: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1454: if (nr && is_row) {
1457: }
1458: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1459: if (nc && is_col) {
1462: }
1464: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1465: return(0);
1466: }
1468: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1469: {
1471: PetscBool flg;
1472: PetscInt i,j,m,mi,*ix;
1475: *ltog = NULL;
1476: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1477: if (islocal[i]) {
1478: ISGetLocalSize(islocal[i],&mi);
1479: flg = PETSC_TRUE; /* We found a non-trivial entry */
1480: } else {
1481: ISGetLocalSize(isglobal[i],&mi);
1482: }
1483: m += mi;
1484: }
1485: if (!flg) return(0);
1487: PetscMalloc1(m,&ix);
1488: for (i=0,m=0; i<n; i++) {
1489: ISLocalToGlobalMapping smap = NULL;
1490: Mat sub = NULL;
1491: PetscSF sf;
1492: PetscLayout map;
1493: const PetscInt *ix2;
1495: if (!colflg) {
1496: MatNestFindNonzeroSubMatRow(A,i,&sub);
1497: } else {
1498: MatNestFindNonzeroSubMatCol(A,i,&sub);
1499: }
1500: if (sub) {
1501: if (!colflg) {
1502: MatGetLocalToGlobalMapping(sub,&smap,NULL);
1503: } else {
1504: MatGetLocalToGlobalMapping(sub,NULL,&smap);
1505: }
1506: }
1507: /*
1508: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1509: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1510: */
1511: ISGetIndices(isglobal[i],&ix2);
1512: if (islocal[i]) {
1513: PetscInt *ilocal,*iremote;
1514: PetscInt mil,nleaves;
1516: ISGetLocalSize(islocal[i],&mi);
1517: if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1518: for (j=0; j<mi; j++) ix[m+j] = j;
1519: ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);
1521: /* PetscSFSetGraphLayout does not like negative indices */
1522: PetscMalloc2(mi,&ilocal,mi,&iremote);
1523: for (j=0, nleaves = 0; j<mi; j++) {
1524: if (ix[m+j] < 0) continue;
1525: ilocal[nleaves] = j;
1526: iremote[nleaves] = ix[m+j];
1527: nleaves++;
1528: }
1529: ISGetLocalSize(isglobal[i],&mil);
1530: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1531: PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1532: PetscLayoutSetLocalSize(map,mil);
1533: PetscLayoutSetUp(map);
1534: PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1535: PetscLayoutDestroy(&map);
1536: PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1537: PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1538: PetscSFDestroy(&sf);
1539: PetscFree2(ilocal,iremote);
1540: } else {
1541: ISGetLocalSize(isglobal[i],&mi);
1542: for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1543: }
1544: ISRestoreIndices(isglobal[i],&ix2);
1545: m += mi;
1546: }
1547: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1548: return(0);
1549: }
1552: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1553: /*
1554: nprocessors = NP
1555: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1556: proc 0: => (g_0,h_0,)
1557: proc 1: => (g_1,h_1,)
1558: ...
1559: proc nprocs-1: => (g_NP-1,h_NP-1,)
1561: proc 0: proc 1: proc nprocs-1:
1562: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1564: proc 0:
1565: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1566: proc 1:
1567: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1569: proc NP-1:
1570: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1571: */
1572: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1573: {
1574: Mat_Nest *vs = (Mat_Nest*)A->data;
1575: PetscInt i,j,offset,n,nsum,bs;
1577: Mat sub = NULL;
1580: PetscMalloc1(nr,&vs->isglobal.row);
1581: PetscMalloc1(nc,&vs->isglobal.col);
1582: if (is_row) { /* valid IS is passed in */
1583: /* refs on is[] are incremeneted */
1584: for (i=0; i<vs->nr; i++) {
1585: PetscObjectReference((PetscObject)is_row[i]);
1587: vs->isglobal.row[i] = is_row[i];
1588: }
1589: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1590: nsum = 0;
1591: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1592: MatNestFindNonzeroSubMatRow(A,i,&sub);
1593: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1594: MatGetLocalSize(sub,&n,NULL);
1595: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1596: nsum += n;
1597: }
1598: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1599: offset -= nsum;
1600: for (i=0; i<vs->nr; i++) {
1601: MatNestFindNonzeroSubMatRow(A,i,&sub);
1602: MatGetLocalSize(sub,&n,NULL);
1603: MatGetBlockSizes(sub,&bs,NULL);
1604: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1605: ISSetBlockSize(vs->isglobal.row[i],bs);
1606: offset += n;
1607: }
1608: }
1610: if (is_col) { /* valid IS is passed in */
1611: /* refs on is[] are incremeneted */
1612: for (j=0; j<vs->nc; j++) {
1613: PetscObjectReference((PetscObject)is_col[j]);
1615: vs->isglobal.col[j] = is_col[j];
1616: }
1617: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1618: offset = A->cmap->rstart;
1619: nsum = 0;
1620: for (j=0; j<vs->nc; j++) {
1621: MatNestFindNonzeroSubMatCol(A,j,&sub);
1622: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1623: MatGetLocalSize(sub,NULL,&n);
1624: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1625: nsum += n;
1626: }
1627: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1628: offset -= nsum;
1629: for (j=0; j<vs->nc; j++) {
1630: MatNestFindNonzeroSubMatCol(A,j,&sub);
1631: MatGetLocalSize(sub,NULL,&n);
1632: MatGetBlockSizes(sub,NULL,&bs);
1633: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1634: ISSetBlockSize(vs->isglobal.col[j],bs);
1635: offset += n;
1636: }
1637: }
1639: /* Set up the local ISs */
1640: PetscMalloc1(vs->nr,&vs->islocal.row);
1641: PetscMalloc1(vs->nc,&vs->islocal.col);
1642: for (i=0,offset=0; i<vs->nr; i++) {
1643: IS isloc;
1644: ISLocalToGlobalMapping rmap = NULL;
1645: PetscInt nlocal,bs;
1646: MatNestFindNonzeroSubMatRow(A,i,&sub);
1647: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1648: if (rmap) {
1649: MatGetBlockSizes(sub,&bs,NULL);
1650: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1651: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1652: ISSetBlockSize(isloc,bs);
1653: } else {
1654: nlocal = 0;
1655: isloc = NULL;
1656: }
1657: vs->islocal.row[i] = isloc;
1658: offset += nlocal;
1659: }
1660: for (i=0,offset=0; i<vs->nc; i++) {
1661: IS isloc;
1662: ISLocalToGlobalMapping cmap = NULL;
1663: PetscInt nlocal,bs;
1664: MatNestFindNonzeroSubMatCol(A,i,&sub);
1665: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1666: if (cmap) {
1667: MatGetBlockSizes(sub,NULL,&bs);
1668: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1669: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1670: ISSetBlockSize(isloc,bs);
1671: } else {
1672: nlocal = 0;
1673: isloc = NULL;
1674: }
1675: vs->islocal.col[i] = isloc;
1676: offset += nlocal;
1677: }
1679: /* Set up the aggregate ISLocalToGlobalMapping */
1680: {
1681: ISLocalToGlobalMapping rmap,cmap;
1682: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1683: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1684: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1685: ISLocalToGlobalMappingDestroy(&rmap);
1686: ISLocalToGlobalMappingDestroy(&cmap);
1687: }
1689: if (PetscDefined(USE_DEBUG)) {
1690: for (i=0; i<vs->nr; i++) {
1691: for (j=0; j<vs->nc; j++) {
1692: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1693: Mat B = vs->m[i][j];
1694: if (!B) continue;
1695: MatGetSize(B,&M,&N);
1696: MatGetLocalSize(B,&m,&n);
1697: ISGetSize(vs->isglobal.row[i],&Mi);
1698: ISGetSize(vs->isglobal.col[j],&Ni);
1699: ISGetLocalSize(vs->isglobal.row[i],&mi);
1700: ISGetLocalSize(vs->isglobal.col[j],&ni);
1701: 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);
1702: 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);
1703: }
1704: }
1705: }
1707: /* Set A->assembled if all non-null blocks are currently assembled */
1708: for (i=0; i<vs->nr; i++) {
1709: for (j=0; j<vs->nc; j++) {
1710: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1711: }
1712: }
1713: A->assembled = PETSC_TRUE;
1714: return(0);
1715: }
1717: /*@C
1718: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1720: Collective on Mat
1722: Input Parameter:
1723: + comm - Communicator for the new Mat
1724: . nr - number of nested row blocks
1725: . is_row - index sets for each nested row block, or NULL to make contiguous
1726: . nc - number of nested column blocks
1727: . is_col - index sets for each nested column block, or NULL to make contiguous
1728: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1730: Output Parameter:
1731: . B - new matrix
1733: Level: advanced
1735: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1736: MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1737: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1738: @*/
1739: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1740: {
1741: Mat A;
1745: *B = NULL;
1746: MatCreate(comm,&A);
1747: MatSetType(A,MATNEST);
1748: A->preallocated = PETSC_TRUE;
1749: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1750: *B = A;
1751: return(0);
1752: }
1754: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1755: {
1756: Mat_Nest *nest = (Mat_Nest*)A->data;
1757: Mat *trans;
1758: PetscScalar **avv;
1759: PetscScalar *vv;
1760: PetscInt **aii,**ajj;
1761: PetscInt *ii,*jj,*ci;
1762: PetscInt nr,nc,nnz,i,j;
1763: PetscBool done;
1767: MatGetSize(A,&nr,&nc);
1768: if (reuse == MAT_REUSE_MATRIX) {
1769: PetscInt rnr;
1771: MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1772: if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1773: if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1774: MatSeqAIJGetArray(*newmat,&vv);
1775: }
1776: /* extract CSR for nested SeqAIJ matrices */
1777: nnz = 0;
1778: PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1779: for (i=0; i<nest->nr; ++i) {
1780: for (j=0; j<nest->nc; ++j) {
1781: Mat B = nest->m[i][j];
1782: if (B) {
1783: PetscScalar *naa;
1784: PetscInt *nii,*njj,nnr;
1785: PetscBool istrans;
1787: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1788: if (istrans) {
1789: Mat Bt;
1791: MatTransposeGetMat(B,&Bt);
1792: MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1793: B = trans[i*nest->nc+j];
1794: }
1795: MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1796: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1797: MatSeqAIJGetArray(B,&naa);
1798: nnz += nii[nnr];
1800: aii[i*nest->nc+j] = nii;
1801: ajj[i*nest->nc+j] = njj;
1802: avv[i*nest->nc+j] = naa;
1803: }
1804: }
1805: }
1806: if (reuse != MAT_REUSE_MATRIX) {
1807: PetscMalloc1(nr+1,&ii);
1808: PetscMalloc1(nnz,&jj);
1809: PetscMalloc1(nnz,&vv);
1810: } else {
1811: if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1812: }
1814: /* new row pointer */
1815: PetscArrayzero(ii,nr+1);
1816: for (i=0; i<nest->nr; ++i) {
1817: PetscInt ncr,rst;
1819: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1820: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1821: for (j=0; j<nest->nc; ++j) {
1822: if (aii[i*nest->nc+j]) {
1823: PetscInt *nii = aii[i*nest->nc+j];
1824: PetscInt ir;
1826: for (ir=rst; ir<ncr+rst; ++ir) {
1827: ii[ir+1] += nii[1]-nii[0];
1828: nii++;
1829: }
1830: }
1831: }
1832: }
1833: for (i=0; i<nr; i++) ii[i+1] += ii[i];
1835: /* construct CSR for the new matrix */
1836: PetscCalloc1(nr,&ci);
1837: for (i=0; i<nest->nr; ++i) {
1838: PetscInt ncr,rst;
1840: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1841: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1842: for (j=0; j<nest->nc; ++j) {
1843: if (aii[i*nest->nc+j]) {
1844: PetscScalar *nvv = avv[i*nest->nc+j];
1845: PetscInt *nii = aii[i*nest->nc+j];
1846: PetscInt *njj = ajj[i*nest->nc+j];
1847: PetscInt ir,cst;
1849: ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1850: for (ir=rst; ir<ncr+rst; ++ir) {
1851: PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1853: for (ij=0;ij<rsize;ij++) {
1854: jj[ist+ij] = *njj+cst;
1855: vv[ist+ij] = *nvv;
1856: njj++;
1857: nvv++;
1858: }
1859: ci[ir] += rsize;
1860: nii++;
1861: }
1862: }
1863: }
1864: }
1865: PetscFree(ci);
1867: /* restore info */
1868: for (i=0; i<nest->nr; ++i) {
1869: for (j=0; j<nest->nc; ++j) {
1870: Mat B = nest->m[i][j];
1871: if (B) {
1872: PetscInt nnr = 0, k = i*nest->nc+j;
1874: B = (trans[k] ? trans[k] : B);
1875: MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1876: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1877: MatSeqAIJRestoreArray(B,&avv[k]);
1878: MatDestroy(&trans[k]);
1879: }
1880: }
1881: }
1882: PetscFree4(aii,ajj,avv,trans);
1884: /* finalize newmat */
1885: if (reuse == MAT_INITIAL_MATRIX) {
1886: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1887: } else if (reuse == MAT_INPLACE_MATRIX) {
1888: Mat B;
1890: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1891: MatHeaderReplace(A,&B);
1892: }
1893: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1894: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1895: {
1896: Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1897: a->free_a = PETSC_TRUE;
1898: a->free_ij = PETSC_TRUE;
1899: }
1900: return(0);
1901: }
1903: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1904: {
1906: Mat_Nest *nest = (Mat_Nest*)A->data;
1907: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1908: PetscInt cstart,cend;
1909: PetscMPIInt size;
1910: Mat C;
1913: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1914: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1915: PetscInt nf;
1916: PetscBool fast;
1918: PetscStrcmp(newtype,MATAIJ,&fast);
1919: if (!fast) {
1920: PetscStrcmp(newtype,MATSEQAIJ,&fast);
1921: }
1922: for (i=0; i<nest->nr && fast; ++i) {
1923: for (j=0; j<nest->nc && fast; ++j) {
1924: Mat B = nest->m[i][j];
1925: if (B) {
1926: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1927: if (!fast) {
1928: PetscBool istrans;
1930: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1931: if (istrans) {
1932: Mat Bt;
1934: MatTransposeGetMat(B,&Bt);
1935: PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1936: }
1937: }
1938: }
1939: }
1940: }
1941: for (i=0, nf=0; i<nest->nr && fast; ++i) {
1942: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1943: if (fast) {
1944: PetscInt f,s;
1946: ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1947: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1948: else {
1949: ISGetSize(nest->isglobal.row[i],&f);
1950: nf += f;
1951: }
1952: }
1953: }
1954: for (i=0, nf=0; i<nest->nc && fast; ++i) {
1955: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1956: if (fast) {
1957: PetscInt f,s;
1959: ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1960: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1961: else {
1962: ISGetSize(nest->isglobal.col[i],&f);
1963: nf += f;
1964: }
1965: }
1966: }
1967: if (fast) {
1968: MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1969: return(0);
1970: }
1971: }
1972: MatGetSize(A,&M,&N);
1973: MatGetLocalSize(A,&m,&n);
1974: MatGetOwnershipRangeColumn(A,&cstart,&cend);
1975: switch (reuse) {
1976: case MAT_INITIAL_MATRIX:
1977: MatCreate(PetscObjectComm((PetscObject)A),&C);
1978: MatSetType(C,newtype);
1979: MatSetSizes(C,m,n,M,N);
1980: *newmat = C;
1981: break;
1982: case MAT_REUSE_MATRIX:
1983: C = *newmat;
1984: break;
1985: default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1986: }
1987: PetscMalloc1(2*m,&dnnz);
1988: onnz = dnnz + m;
1989: for (k=0; k<m; k++) {
1990: dnnz[k] = 0;
1991: onnz[k] = 0;
1992: }
1993: for (j=0; j<nest->nc; ++j) {
1994: IS bNis;
1995: PetscInt bN;
1996: const PetscInt *bNindices;
1997: /* Using global column indices and ISAllGather() is not scalable. */
1998: ISAllGather(nest->isglobal.col[j], &bNis);
1999: ISGetSize(bNis, &bN);
2000: ISGetIndices(bNis,&bNindices);
2001: for (i=0; i<nest->nr; ++i) {
2002: PetscSF bmsf;
2003: PetscSFNode *iremote;
2004: Mat B;
2005: PetscInt bm, *sub_dnnz,*sub_onnz, br;
2006: const PetscInt *bmindices;
2007: B = nest->m[i][j];
2008: if (!B) continue;
2009: ISGetLocalSize(nest->isglobal.row[i],&bm);
2010: ISGetIndices(nest->isglobal.row[i],&bmindices);
2011: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2012: PetscMalloc1(bm,&iremote);
2013: PetscMalloc1(bm,&sub_dnnz);
2014: PetscMalloc1(bm,&sub_onnz);
2015: for (k = 0; k < bm; ++k){
2016: sub_dnnz[k] = 0;
2017: sub_onnz[k] = 0;
2018: }
2019: /*
2020: Locate the owners for all of the locally-owned global row indices for this row block.
2021: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2022: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2023: */
2024: MatGetOwnershipRange(B,&rstart,NULL);
2025: for (br = 0; br < bm; ++br) {
2026: PetscInt row = bmindices[br], brncols, col;
2027: const PetscInt *brcols;
2028: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2029: PetscMPIInt rowowner = 0;
2030: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2031: /* how many roots */
2032: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2033: /* get nonzero pattern */
2034: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2035: for (k=0; k<brncols; k++) {
2036: col = bNindices[brcols[k]];
2037: if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2038: sub_dnnz[br]++;
2039: } else {
2040: sub_onnz[br]++;
2041: }
2042: }
2043: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2044: }
2045: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2046: /* bsf will have to take care of disposing of bedges. */
2047: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2048: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2049: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2050: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2051: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2052: PetscFree(sub_dnnz);
2053: PetscFree(sub_onnz);
2054: PetscSFDestroy(&bmsf);
2055: }
2056: ISRestoreIndices(bNis,&bNindices);
2057: ISDestroy(&bNis);
2058: }
2059: /* Resize preallocation if overestimated */
2060: for (i=0;i<m;i++) {
2061: dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2062: onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2063: }
2064: MatSeqAIJSetPreallocation(C,0,dnnz);
2065: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2066: PetscFree(dnnz);
2068: /* Fill by row */
2069: for (j=0; j<nest->nc; ++j) {
2070: /* Using global column indices and ISAllGather() is not scalable. */
2071: IS bNis;
2072: PetscInt bN;
2073: const PetscInt *bNindices;
2074: ISAllGather(nest->isglobal.col[j], &bNis);
2075: ISGetSize(bNis,&bN);
2076: ISGetIndices(bNis,&bNindices);
2077: for (i=0; i<nest->nr; ++i) {
2078: Mat B;
2079: PetscInt bm, br;
2080: const PetscInt *bmindices;
2081: B = nest->m[i][j];
2082: if (!B) continue;
2083: ISGetLocalSize(nest->isglobal.row[i],&bm);
2084: ISGetIndices(nest->isglobal.row[i],&bmindices);
2085: MatGetOwnershipRange(B,&rstart,NULL);
2086: for (br = 0; br < bm; ++br) {
2087: PetscInt row = bmindices[br], brncols, *cols;
2088: const PetscInt *brcols;
2089: const PetscScalar *brcoldata;
2090: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2091: PetscMalloc1(brncols,&cols);
2092: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2093: /*
2094: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2095: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2096: */
2097: MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2098: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2099: PetscFree(cols);
2100: }
2101: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2102: }
2103: ISRestoreIndices(bNis,&bNindices);
2104: ISDestroy(&bNis);
2105: }
2106: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2107: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2108: return(0);
2109: }
2111: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2112: {
2113: Mat_Nest *bA = (Mat_Nest*)mat->data;
2114: MatOperation opAdd;
2115: PetscInt i,j,nr = bA->nr,nc = bA->nc;
2116: PetscBool flg;
2120: *has = PETSC_FALSE;
2121: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2122: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2123: for (j=0; j<nc; j++) {
2124: for (i=0; i<nr; i++) {
2125: if (!bA->m[i][j]) continue;
2126: MatHasOperation(bA->m[i][j],opAdd,&flg);
2127: if (!flg) return(0);
2128: }
2129: }
2130: }
2131: if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2132: return(0);
2133: }
2135: /*MC
2136: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2138: Level: intermediate
2140: Notes:
2141: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2142: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2143: It is usually used with DMComposite and DMCreateMatrix()
2145: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2146: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2147: than the nest matrix.
2149: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2150: VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2151: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2152: M*/
2153: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2154: {
2155: Mat_Nest *s;
2159: PetscNewLog(A,&s);
2160: A->data = (void*)s;
2162: s->nr = -1;
2163: s->nc = -1;
2164: s->m = NULL;
2165: s->splitassembly = PETSC_FALSE;
2167: PetscMemzero(A->ops,sizeof(*A->ops));
2169: A->ops->mult = MatMult_Nest;
2170: A->ops->multadd = MatMultAdd_Nest;
2171: A->ops->multtranspose = MatMultTranspose_Nest;
2172: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2173: A->ops->transpose = MatTranspose_Nest;
2174: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2175: A->ops->assemblyend = MatAssemblyEnd_Nest;
2176: A->ops->zeroentries = MatZeroEntries_Nest;
2177: A->ops->copy = MatCopy_Nest;
2178: A->ops->axpy = MatAXPY_Nest;
2179: A->ops->duplicate = MatDuplicate_Nest;
2180: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2181: A->ops->destroy = MatDestroy_Nest;
2182: A->ops->view = MatView_Nest;
2183: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2184: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2185: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2186: A->ops->getdiagonal = MatGetDiagonal_Nest;
2187: A->ops->diagonalscale = MatDiagonalScale_Nest;
2188: A->ops->scale = MatScale_Nest;
2189: A->ops->shift = MatShift_Nest;
2190: A->ops->diagonalset = MatDiagonalSet_Nest;
2191: A->ops->setrandom = MatSetRandom_Nest;
2192: A->ops->hasoperation = MatHasOperation_Nest;
2193: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2195: A->spptr = NULL;
2196: A->assembled = PETSC_FALSE;
2198: /* expose Nest api's */
2199: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2200: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2201: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2202: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
2203: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
2204: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2205: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
2206: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2207: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2208: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2209: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2210: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);
2211: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2212: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2213: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);
2215: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2216: return(0);
2217: }