Actual source code: matnest.c
petsc-3.13.6 2020-09-29
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 MatMatMultNumeric_Nest_Dense(Mat A,Mat B,Mat C)
96: {
97: Mat_Nest *bA = (Mat_Nest*)A->data;
98: PetscContainer container;
99: Nest_Dense *contents;
100: Mat viewB,viewC,seq,productB,workC;
101: const PetscScalar *barray;
102: PetscScalar *carray;
103: PetscInt i,j,M,N,nr = bA->nr,nc = bA->nc,ldb,ldc;
104: PetscErrorCode ierr;
107: PetscObjectQuery((PetscObject)C,"workC",(PetscObject*)&container);
108: if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exist");
109: PetscContainerGetPointer(container,(void**)&contents);
110: MatDenseGetLDA(B,&ldb);
111: MatDenseGetLDA(C,&ldc);
112: MatGetSize(B,NULL,&N);
113: MatZeroEntries(C);
114: MatDenseGetArrayRead(B,&barray);
115: MatDenseGetArray(C,&carray);
116: for (i=0; i<nr; i++) {
117: ISGetSize(bA->isglobal.row[i],&M);
118: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
119: MatDenseGetLocalMatrix(viewC,&seq);
120: MatSeqDenseSetLDA(seq,ldc);
121: for (j=0; j<nc; j++) {
122: if (!bA->m[i][j]) continue;
123: ISGetSize(bA->isglobal.col[j],&M);
124: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
125: MatDenseGetLocalMatrix(viewB,&seq);
126: MatSeqDenseSetLDA(seq,ldb);
128: /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
129: workC = contents->workC[i*nc + j];
130: productB = workC->product->B;
131: workC->product->B = viewB; /* use newly created dense matrix viewB */
132: (*workC->ops->productnumeric)(workC);
133: MatDestroy(&viewB);
134: workC->product->B = productB; /* resume original B */
136: /* C[i] <- workC + C[i] */
137: MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
138: }
139: MatDestroy(&viewC);
140: }
141: MatDenseRestoreArray(C,&carray);
142: MatDenseRestoreArrayRead(B,&barray);
144: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
145: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
146: return(0);
147: }
149: PetscErrorCode MatNest_DenseDestroy(void *ctx)
150: {
151: Nest_Dense *contents = (Nest_Dense*)ctx;
152: PetscInt i;
156: PetscFree(contents->tarray);
157: for (i=0; i<contents->k; i++) {
158: MatDestroy(contents->workC + i);
159: }
160: PetscFree3(contents->dm,contents->dn,contents->workC);
161: PetscFree(contents);
162: return(0);
163: }
165: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_Nest_Dense(Mat A,Mat B,PetscReal fill,Mat C)
166: {
167: Mat_Nest *bA = (Mat_Nest*)A->data;
168: Mat viewB,viewSeq,workC;
169: const PetscScalar *barray;
170: PetscInt i,j,M,N,m,nr = bA->nr,nc = bA->nc,maxm = 0,ldb;
171: PetscContainer container;
172: Nest_Dense *contents=NULL;
173: PetscErrorCode ierr;
176: MatGetSize(B,NULL,&N);
177: if (!C->assembled) {
178: MatGetLocalSize(A,&m,NULL);
179: MatGetSize(A,&M,NULL);
181: MatSetSizes(C,m,PETSC_DECIDE,M,N);
182: MatSetType(C,MATDENSE);
183: MatSetUp(C);
184: }
186: PetscNew(&contents);
187: PetscContainerCreate(PetscObjectComm((PetscObject)A),&container);
188: PetscContainerSetPointer(container,contents);
189: PetscContainerSetUserDestroy(container,MatNest_DenseDestroy);
190: PetscObjectCompose((PetscObject)C,"workC",(PetscObject)container);
191: PetscContainerDestroy(&container);
192: PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
193: contents->k = nr*nc;
194: for (i=0; i<nr; i++) {
195: ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
196: maxm = PetscMax(maxm,contents->dm[i+1]);
197: contents->dm[i+1] += contents->dm[i];
198: }
199: for (i=0; i<nc; i++) {
200: ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
201: contents->dn[i+1] += contents->dn[i];
202: }
203: PetscMalloc1(maxm*N,&contents->tarray);
204: MatDenseGetLDA(B,&ldb);
205: MatGetSize(B,NULL,&N);
206: MatDenseGetArrayRead(B,&barray);
207: /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
208: for (j=0; j<nc; j++) {
209: ISGetSize(bA->isglobal.col[j],&M);
210: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
211: MatDenseGetLocalMatrix(viewB,&viewSeq);
212: MatSeqDenseSetLDA(viewSeq,ldb);
213: for (i=0; i<nr; i++) {
214: if (!bA->m[i][j]) continue;
215: /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
217: MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
218: workC = contents->workC[i*nc + j];
219: MatProductSetType(workC,MATPRODUCT_AB);
220: MatProductSetAlgorithm(workC,"default");
221: MatProductSetFill(workC,fill);
222: MatProductSetFromOptions(workC);
223: MatProductSymbolic(workC);
225: MatDenseGetLocalMatrix(workC,&viewSeq);
226: /* free the memory allocated in MatMatMultSymbolic, since tarray will be shared by all Mat */
227: MatSeqDenseSetPreallocation(viewSeq,contents->tarray);
228: }
229: MatDestroy(&viewB);
230: }
231: MatDenseRestoreArrayRead(B,&barray);
233: C->ops->matmultnumeric = MatMatMultNumeric_Nest_Dense;
234: return(0);
235: }
237: /* --------------------------------------------------------- */
238: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
239: {
241: C->ops->matmultsymbolic = MatMatMultSymbolic_Nest_Dense;
242: C->ops->productsymbolic = MatProductSymbolic_AB;
243: return(0);
244: }
246: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
247: {
249: Mat_Product *product = C->product;
252: if (product->type == MATPRODUCT_AB) {
253: MatProductSetFromOptions_Nest_Dense_AB(C);
254: } else SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for Nest and Dense matrices",MatProductTypes[product->type]);
255: return(0);
256: }
257: /* --------------------------------------------------------- */
259: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
260: {
261: Mat_Nest *bA = (Mat_Nest*)A->data;
262: Vec *bx = bA->left,*by = bA->right;
263: PetscInt i,j,nr = bA->nr,nc = bA->nc;
267: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
268: for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
269: for (j=0; j<nc; j++) {
270: VecZeroEntries(by[j]);
271: for (i=0; i<nr; i++) {
272: if (!bA->m[i][j]) continue;
273: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
274: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
275: }
276: }
277: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
278: for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
279: return(0);
280: }
282: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
283: {
284: Mat_Nest *bA = (Mat_Nest*)A->data;
285: Vec *bx = bA->left,*bz = bA->right;
286: PetscInt i,j,nr = bA->nr,nc = bA->nc;
290: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
291: for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
292: for (j=0; j<nc; j++) {
293: if (y != z) {
294: Vec by;
295: VecGetSubVector(y,bA->isglobal.col[j],&by);
296: VecCopy(by,bz[j]);
297: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
298: }
299: for (i=0; i<nr; i++) {
300: if (!bA->m[i][j]) continue;
301: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
302: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
303: }
304: }
305: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
306: for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
307: return(0);
308: }
310: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
311: {
312: Mat_Nest *bA = (Mat_Nest*)A->data, *bC;
313: Mat C;
314: PetscInt i,j,nr = bA->nr,nc = bA->nc;
318: if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");
320: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
321: Mat *subs;
322: IS *is_row,*is_col;
324: PetscCalloc1(nr * nc,&subs);
325: PetscMalloc2(nr,&is_row,nc,&is_col);
326: MatNestGetISs(A,is_row,is_col);
327: if (reuse == MAT_INPLACE_MATRIX) {
328: for (i=0; i<nr; i++) {
329: for (j=0; j<nc; j++) {
330: subs[i + nr * j] = bA->m[i][j];
331: }
332: }
333: }
335: MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
336: PetscFree(subs);
337: PetscFree2(is_row,is_col);
338: } else {
339: C = *B;
340: }
342: bC = (Mat_Nest*)C->data;
343: for (i=0; i<nr; i++) {
344: for (j=0; j<nc; j++) {
345: if (bA->m[i][j]) {
346: MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
347: } else {
348: bC->m[j][i] = NULL;
349: }
350: }
351: }
353: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
354: *B = C;
355: } else {
356: MatHeaderMerge(A, &C);
357: }
358: return(0);
359: }
361: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
362: {
364: IS *lst = *list;
365: PetscInt i;
368: if (!lst) return(0);
369: for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
370: PetscFree(lst);
371: *list = NULL;
372: return(0);
373: }
375: static PetscErrorCode MatReset_Nest(Mat A)
376: {
377: Mat_Nest *vs = (Mat_Nest*)A->data;
378: PetscInt i,j;
382: /* release the matrices and the place holders */
383: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
384: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
385: MatNestDestroyISList(vs->nr,&vs->islocal.row);
386: MatNestDestroyISList(vs->nc,&vs->islocal.col);
388: PetscFree(vs->row_len);
389: PetscFree(vs->col_len);
390: PetscFree(vs->nnzstate);
392: PetscFree2(vs->left,vs->right);
394: /* release the matrices and the place holders */
395: if (vs->m) {
396: for (i=0; i<vs->nr; i++) {
397: for (j=0; j<vs->nc; j++) {
398: MatDestroy(&vs->m[i][j]);
399: }
400: PetscFree(vs->m[i]);
401: }
402: PetscFree(vs->m);
403: }
405: /* restore defaults */
406: vs->nr = 0;
407: vs->nc = 0;
408: vs->splitassembly = PETSC_FALSE;
409: return(0);
410: }
412: static PetscErrorCode MatDestroy_Nest(Mat A)
413: {
416: MatReset_Nest(A);
417: PetscFree(A->data);
419: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
420: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
421: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
422: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
423: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
424: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
425: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
426: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
427: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
428: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
429: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
430: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
431: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
432: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
433: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
434: return(0);
435: }
437: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
438: {
439: Mat_Nest *vs = (Mat_Nest*)mat->data;
440: PetscInt i;
444: if (dd) *dd = 0;
445: if (!vs->nr) {
446: *missing = PETSC_TRUE;
447: return(0);
448: }
449: *missing = PETSC_FALSE;
450: for (i = 0; i < vs->nr && !(*missing); i++) {
451: *missing = PETSC_TRUE;
452: if (vs->m[i][i]) {
453: MatMissingDiagonal(vs->m[i][i],missing,NULL);
454: if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
455: }
456: }
457: return(0);
458: }
460: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
461: {
462: Mat_Nest *vs = (Mat_Nest*)A->data;
463: PetscInt i,j;
465: PetscBool nnzstate = PETSC_FALSE;
468: for (i=0; i<vs->nr; i++) {
469: for (j=0; j<vs->nc; j++) {
470: PetscObjectState subnnzstate = 0;
471: if (vs->m[i][j]) {
472: MatAssemblyBegin(vs->m[i][j],type);
473: if (!vs->splitassembly) {
474: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
475: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
476: * already performing an assembly, but the result would by more complicated and appears to offer less
477: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
478: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
479: */
480: MatAssemblyEnd(vs->m[i][j],type);
481: MatGetNonzeroState(vs->m[i][j],&subnnzstate);
482: }
483: }
484: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
485: vs->nnzstate[i*vs->nc+j] = subnnzstate;
486: }
487: }
488: if (nnzstate) A->nonzerostate++;
489: return(0);
490: }
492: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
493: {
494: Mat_Nest *vs = (Mat_Nest*)A->data;
495: PetscInt i,j;
499: for (i=0; i<vs->nr; i++) {
500: for (j=0; j<vs->nc; j++) {
501: if (vs->m[i][j]) {
502: if (vs->splitassembly) {
503: MatAssemblyEnd(vs->m[i][j],type);
504: }
505: }
506: }
507: }
508: return(0);
509: }
511: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
512: {
514: Mat_Nest *vs = (Mat_Nest*)A->data;
515: PetscInt j;
516: Mat sub;
519: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
520: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
521: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
522: *B = sub;
523: return(0);
524: }
526: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
527: {
529: Mat_Nest *vs = (Mat_Nest*)A->data;
530: PetscInt i;
531: Mat sub;
534: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
535: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
536: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
537: *B = sub;
538: return(0);
539: }
541: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
542: {
544: PetscInt i;
545: PetscBool flg;
551: *found = -1;
552: for (i=0; i<n; i++) {
553: if (!list[i]) continue;
554: ISEqualUnsorted(list[i],is,&flg);
555: if (flg) {
556: *found = i;
557: return(0);
558: }
559: }
560: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
561: return(0);
562: }
564: /* Get a block row as a new MatNest */
565: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
566: {
567: Mat_Nest *vs = (Mat_Nest*)A->data;
568: char keyname[256];
572: *B = NULL;
573: PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
574: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
575: if (*B) return(0);
577: MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);
579: (*B)->assembled = A->assembled;
581: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
582: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
583: return(0);
584: }
586: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
587: {
588: Mat_Nest *vs = (Mat_Nest*)A->data;
590: PetscInt row,col;
591: PetscBool same,isFullCol,isFullColGlobal;
594: /* Check if full column space. This is a hack */
595: isFullCol = PETSC_FALSE;
596: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
597: if (same) {
598: PetscInt n,first,step,i,an,am,afirst,astep;
599: ISStrideGetInfo(iscol,&first,&step);
600: ISGetLocalSize(iscol,&n);
601: isFullCol = PETSC_TRUE;
602: for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
603: PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
604: ISGetLocalSize(is->col[i],&am);
605: if (same) {
606: ISStrideGetInfo(is->col[i],&afirst,&astep);
607: if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
608: } else isFullCol = PETSC_FALSE;
609: an += am;
610: }
611: if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
612: }
613: MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));
615: if (isFullColGlobal && vs->nc > 1) {
616: PetscInt row;
617: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
618: MatNestGetRow(A,row,B);
619: } else {
620: MatNestFindIS(A,vs->nr,is->row,isrow,&row);
621: MatNestFindIS(A,vs->nc,is->col,iscol,&col);
622: if (!vs->m[row][col]) {
623: PetscInt lr,lc;
625: MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
626: ISGetLocalSize(vs->isglobal.row[row],&lr);
627: ISGetLocalSize(vs->isglobal.col[col],&lc);
628: MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
629: MatSetType(vs->m[row][col],MATAIJ);
630: MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
631: MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
632: MatSetUp(vs->m[row][col]);
633: MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
634: MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
635: }
636: *B = vs->m[row][col];
637: }
638: return(0);
639: }
641: /*
642: TODO: This does not actually returns a submatrix we can modify
643: */
644: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
645: {
647: Mat_Nest *vs = (Mat_Nest*)A->data;
648: Mat sub;
651: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
652: switch (reuse) {
653: case MAT_INITIAL_MATRIX:
654: if (sub) { PetscObjectReference((PetscObject)sub); }
655: *B = sub;
656: break;
657: case MAT_REUSE_MATRIX:
658: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
659: break;
660: case MAT_IGNORE_MATRIX: /* Nothing to do */
661: break;
662: case MAT_INPLACE_MATRIX: /* Nothing to do */
663: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
664: break;
665: }
666: return(0);
667: }
669: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
670: {
672: Mat_Nest *vs = (Mat_Nest*)A->data;
673: Mat sub;
676: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
677: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
678: if (sub) {PetscObjectReference((PetscObject)sub);}
679: *B = sub;
680: return(0);
681: }
683: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
684: {
686: Mat_Nest *vs = (Mat_Nest*)A->data;
687: Mat sub;
690: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
691: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
692: if (sub) {
693: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
694: MatDestroy(B);
695: }
696: return(0);
697: }
699: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
700: {
701: Mat_Nest *bA = (Mat_Nest*)A->data;
702: PetscInt i;
706: for (i=0; i<bA->nr; i++) {
707: Vec bv;
708: VecGetSubVector(v,bA->isglobal.row[i],&bv);
709: if (bA->m[i][i]) {
710: MatGetDiagonal(bA->m[i][i],bv);
711: } else {
712: VecSet(bv,0.0);
713: }
714: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
715: }
716: return(0);
717: }
719: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
720: {
721: Mat_Nest *bA = (Mat_Nest*)A->data;
722: Vec bl,*br;
723: PetscInt i,j;
727: PetscCalloc1(bA->nc,&br);
728: if (r) {
729: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
730: }
731: bl = NULL;
732: for (i=0; i<bA->nr; i++) {
733: if (l) {
734: VecGetSubVector(l,bA->isglobal.row[i],&bl);
735: }
736: for (j=0; j<bA->nc; j++) {
737: if (bA->m[i][j]) {
738: MatDiagonalScale(bA->m[i][j],bl,br[j]);
739: }
740: }
741: if (l) {
742: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
743: }
744: }
745: if (r) {
746: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
747: }
748: PetscFree(br);
749: return(0);
750: }
752: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
753: {
754: Mat_Nest *bA = (Mat_Nest*)A->data;
755: PetscInt i,j;
759: for (i=0; i<bA->nr; i++) {
760: for (j=0; j<bA->nc; j++) {
761: if (bA->m[i][j]) {
762: MatScale(bA->m[i][j],a);
763: }
764: }
765: }
766: return(0);
767: }
769: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
770: {
771: Mat_Nest *bA = (Mat_Nest*)A->data;
772: PetscInt i;
774: PetscBool nnzstate = PETSC_FALSE;
777: for (i=0; i<bA->nr; i++) {
778: PetscObjectState subnnzstate = 0;
779: 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);
780: MatShift(bA->m[i][i],a);
781: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
782: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
783: bA->nnzstate[i*bA->nc+i] = subnnzstate;
784: }
785: if (nnzstate) A->nonzerostate++;
786: return(0);
787: }
789: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
790: {
791: Mat_Nest *bA = (Mat_Nest*)A->data;
792: PetscInt i;
794: PetscBool nnzstate = PETSC_FALSE;
797: for (i=0; i<bA->nr; i++) {
798: PetscObjectState subnnzstate = 0;
799: Vec bv;
800: VecGetSubVector(D,bA->isglobal.row[i],&bv);
801: if (bA->m[i][i]) {
802: MatDiagonalSet(bA->m[i][i],bv,is);
803: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
804: }
805: VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
806: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
807: bA->nnzstate[i*bA->nc+i] = subnnzstate;
808: }
809: if (nnzstate) A->nonzerostate++;
810: return(0);
811: }
813: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
814: {
815: Mat_Nest *bA = (Mat_Nest*)A->data;
816: PetscInt i,j;
820: for (i=0; i<bA->nr; i++) {
821: for (j=0; j<bA->nc; j++) {
822: if (bA->m[i][j]) {
823: MatSetRandom(bA->m[i][j],rctx);
824: }
825: }
826: }
827: return(0);
828: }
830: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
831: {
832: Mat_Nest *bA = (Mat_Nest*)A->data;
833: Vec *L,*R;
834: MPI_Comm comm;
835: PetscInt i,j;
839: PetscObjectGetComm((PetscObject)A,&comm);
840: if (right) {
841: /* allocate R */
842: PetscMalloc1(bA->nc, &R);
843: /* Create the right vectors */
844: for (j=0; j<bA->nc; j++) {
845: for (i=0; i<bA->nr; i++) {
846: if (bA->m[i][j]) {
847: MatCreateVecs(bA->m[i][j],&R[j],NULL);
848: break;
849: }
850: }
851: if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
852: }
853: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
854: /* hand back control to the nest vector */
855: for (j=0; j<bA->nc; j++) {
856: VecDestroy(&R[j]);
857: }
858: PetscFree(R);
859: }
861: if (left) {
862: /* allocate L */
863: PetscMalloc1(bA->nr, &L);
864: /* Create the left vectors */
865: for (i=0; i<bA->nr; i++) {
866: for (j=0; j<bA->nc; j++) {
867: if (bA->m[i][j]) {
868: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
869: break;
870: }
871: }
872: if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
873: }
875: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
876: for (i=0; i<bA->nr; i++) {
877: VecDestroy(&L[i]);
878: }
880: PetscFree(L);
881: }
882: return(0);
883: }
885: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
886: {
887: Mat_Nest *bA = (Mat_Nest*)A->data;
888: PetscBool isascii,viewSub = PETSC_FALSE;
889: PetscInt i,j;
893: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
894: if (isascii) {
896: PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
897: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
898: PetscViewerASCIIPushTab(viewer);
899: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);
901: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
902: for (i=0; i<bA->nr; i++) {
903: for (j=0; j<bA->nc; j++) {
904: MatType type;
905: char name[256] = "",prefix[256] = "";
906: PetscInt NR,NC;
907: PetscBool isNest = PETSC_FALSE;
909: if (!bA->m[i][j]) {
910: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
911: continue;
912: }
913: MatGetSize(bA->m[i][j],&NR,&NC);
914: MatGetType(bA->m[i][j], &type);
915: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
916: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
917: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
919: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
921: if (isNest || viewSub) {
922: PetscViewerASCIIPushTab(viewer); /* push1 */
923: MatView(bA->m[i][j],viewer);
924: PetscViewerASCIIPopTab(viewer); /* pop1 */
925: }
926: }
927: }
928: PetscViewerASCIIPopTab(viewer); /* pop0 */
929: }
930: return(0);
931: }
933: static PetscErrorCode MatZeroEntries_Nest(Mat A)
934: {
935: Mat_Nest *bA = (Mat_Nest*)A->data;
936: PetscInt i,j;
940: for (i=0; i<bA->nr; i++) {
941: for (j=0; j<bA->nc; j++) {
942: if (!bA->m[i][j]) continue;
943: MatZeroEntries(bA->m[i][j]);
944: }
945: }
946: return(0);
947: }
949: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
950: {
951: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
952: PetscInt i,j,nr = bA->nr,nc = bA->nc;
954: PetscBool nnzstate = PETSC_FALSE;
957: 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);
958: for (i=0; i<nr; i++) {
959: for (j=0; j<nc; j++) {
960: PetscObjectState subnnzstate = 0;
961: if (bA->m[i][j] && bB->m[i][j]) {
962: MatCopy(bA->m[i][j],bB->m[i][j],str);
963: } 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);
964: MatGetNonzeroState(bB->m[i][j],&subnnzstate);
965: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
966: bB->nnzstate[i*nc+j] = subnnzstate;
967: }
968: }
969: if (nnzstate) B->nonzerostate++;
970: return(0);
971: }
973: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
974: {
975: Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
976: PetscInt i,j,nr = bY->nr,nc = bY->nc;
978: PetscBool nnzstate = PETSC_FALSE;
981: 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);
982: for (i=0; i<nr; i++) {
983: for (j=0; j<nc; j++) {
984: PetscObjectState subnnzstate = 0;
985: if (bY->m[i][j] && bX->m[i][j]) {
986: MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
987: } else if (bX->m[i][j]) {
988: Mat M;
990: 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);
991: MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
992: MatNestSetSubMat(Y,i,j,M);
993: MatDestroy(&M);
994: }
995: if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
996: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
997: bY->nnzstate[i*nc+j] = subnnzstate;
998: }
999: }
1000: if (nnzstate) Y->nonzerostate++;
1001: return(0);
1002: }
1004: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1005: {
1006: Mat_Nest *bA = (Mat_Nest*)A->data;
1007: Mat *b;
1008: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1012: PetscMalloc1(nr*nc,&b);
1013: for (i=0; i<nr; i++) {
1014: for (j=0; j<nc; j++) {
1015: if (bA->m[i][j]) {
1016: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1017: } else {
1018: b[i*nc+j] = NULL;
1019: }
1020: }
1021: }
1022: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1023: /* Give the new MatNest exclusive ownership */
1024: for (i=0; i<nr*nc; i++) {
1025: MatDestroy(&b[i]);
1026: }
1027: PetscFree(b);
1029: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1030: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1031: return(0);
1032: }
1034: /* nest api */
1035: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1036: {
1037: Mat_Nest *bA = (Mat_Nest*)A->data;
1040: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1041: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1042: *mat = bA->m[idxm][jdxm];
1043: return(0);
1044: }
1046: /*@
1047: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
1049: Not collective
1051: Input Parameters:
1052: + A - nest matrix
1053: . idxm - index of the matrix within the nest matrix
1054: - jdxm - index of the matrix within the nest matrix
1056: Output Parameter:
1057: . sub - matrix at index idxm,jdxm within the nest matrix
1059: Level: developer
1061: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1062: MatNestGetLocalISs(), MatNestGetISs()
1063: @*/
1064: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1065: {
1069: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1070: return(0);
1071: }
1073: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1074: {
1075: Mat_Nest *bA = (Mat_Nest*)A->data;
1076: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1080: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1081: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1082: MatGetLocalSize(mat,&m,&n);
1083: MatGetSize(mat,&M,&N);
1084: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1085: ISGetSize(bA->isglobal.row[idxm],&Mi);
1086: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1087: ISGetSize(bA->isglobal.col[jdxm],&Ni);
1088: 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);
1089: 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);
1091: /* do not increase object state */
1092: if (mat == bA->m[idxm][jdxm]) return(0);
1094: PetscObjectReference((PetscObject)mat);
1095: MatDestroy(&bA->m[idxm][jdxm]);
1096: bA->m[idxm][jdxm] = mat;
1097: PetscObjectStateIncrease((PetscObject)A);
1098: MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1099: A->nonzerostate++;
1100: return(0);
1101: }
1103: /*@
1104: MatNestSetSubMat - Set a single submatrix in the nest matrix.
1106: Logically collective on the submatrix communicator
1108: Input Parameters:
1109: + A - nest matrix
1110: . idxm - index of the matrix within the nest matrix
1111: . jdxm - index of the matrix within the nest matrix
1112: - sub - matrix at index idxm,jdxm within the nest matrix
1114: Notes:
1115: The new submatrix must have the same size and communicator as that block of the nest.
1117: This increments the reference count of the submatrix.
1119: Level: developer
1121: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1122: MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1123: @*/
1124: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1125: {
1129: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1130: return(0);
1131: }
1133: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1134: {
1135: Mat_Nest *bA = (Mat_Nest*)A->data;
1138: if (M) *M = bA->nr;
1139: if (N) *N = bA->nc;
1140: if (mat) *mat = bA->m;
1141: return(0);
1142: }
1144: /*@C
1145: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
1147: Not collective
1149: Input Parameters:
1150: . A - nest matrix
1152: Output Parameter:
1153: + M - number of rows in the nest matrix
1154: . N - number of cols in the nest matrix
1155: - mat - 2d array of matrices
1157: Notes:
1159: The user should not free the array mat.
1161: In Fortran, this routine has a calling sequence
1162: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1163: where the space allocated for the optional argument mat is assumed large enough (if provided).
1165: Level: developer
1167: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1168: MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1169: @*/
1170: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1171: {
1175: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1176: return(0);
1177: }
1179: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1180: {
1181: Mat_Nest *bA = (Mat_Nest*)A->data;
1184: if (M) *M = bA->nr;
1185: if (N) *N = bA->nc;
1186: return(0);
1187: }
1189: /*@
1190: MatNestGetSize - Returns the size of the nest matrix.
1192: Not collective
1194: Input Parameters:
1195: . A - nest matrix
1197: Output Parameter:
1198: + M - number of rows in the nested mat
1199: - N - number of cols in the nested mat
1201: Notes:
1203: Level: developer
1205: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1206: MatNestGetISs()
1207: @*/
1208: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1209: {
1213: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1214: return(0);
1215: }
1217: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1218: {
1219: Mat_Nest *vs = (Mat_Nest*)A->data;
1220: PetscInt i;
1223: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1224: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1225: return(0);
1226: }
1228: /*@C
1229: MatNestGetISs - Returns the index sets partitioning the row and column spaces
1231: Not collective
1233: Input Parameters:
1234: . A - nest matrix
1236: Output Parameter:
1237: + rows - array of row index sets
1238: - cols - array of column index sets
1240: Level: advanced
1242: Notes:
1243: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1245: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1246: MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1247: @*/
1248: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1249: {
1254: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1255: return(0);
1256: }
1258: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1259: {
1260: Mat_Nest *vs = (Mat_Nest*)A->data;
1261: PetscInt i;
1264: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1265: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1266: return(0);
1267: }
1269: /*@C
1270: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1272: Not collective
1274: Input Parameters:
1275: . A - nest matrix
1277: Output Parameter:
1278: + rows - array of row index sets (or NULL to ignore)
1279: - cols - array of column index sets (or NULL to ignore)
1281: Level: advanced
1283: Notes:
1284: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1286: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1287: MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1288: @*/
1289: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1290: {
1295: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1296: return(0);
1297: }
1299: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1300: {
1302: PetscBool flg;
1305: PetscStrcmp(vtype,VECNEST,&flg);
1306: /* In reality, this only distinguishes VECNEST and "other" */
1307: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1308: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1309: return(0);
1310: }
1312: /*@C
1313: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1315: Not collective
1317: Input Parameters:
1318: + A - nest matrix
1319: - vtype - type to use for creating vectors
1321: Notes:
1323: Level: developer
1325: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1326: @*/
1327: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1328: {
1332: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1333: return(0);
1334: }
1336: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1337: {
1338: Mat_Nest *s = (Mat_Nest*)A->data;
1339: PetscInt i,j,m,n,M,N;
1341: PetscBool cong;
1344: MatReset_Nest(A);
1346: s->nr = nr;
1347: s->nc = nc;
1349: /* Create space for submatrices */
1350: PetscMalloc1(nr,&s->m);
1351: for (i=0; i<nr; i++) {
1352: PetscMalloc1(nc,&s->m[i]);
1353: }
1354: for (i=0; i<nr; i++) {
1355: for (j=0; j<nc; j++) {
1356: s->m[i][j] = a[i*nc+j];
1357: if (a[i*nc+j]) {
1358: PetscObjectReference((PetscObject)a[i*nc+j]);
1359: }
1360: }
1361: }
1363: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1365: PetscMalloc1(nr,&s->row_len);
1366: PetscMalloc1(nc,&s->col_len);
1367: for (i=0; i<nr; i++) s->row_len[i]=-1;
1368: for (j=0; j<nc; j++) s->col_len[j]=-1;
1370: PetscCalloc1(nr*nc,&s->nnzstate);
1371: for (i=0; i<nr; i++) {
1372: for (j=0; j<nc; j++) {
1373: if (s->m[i][j]) {
1374: MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1375: }
1376: }
1377: }
1379: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1381: PetscLayoutSetSize(A->rmap,M);
1382: PetscLayoutSetLocalSize(A->rmap,m);
1383: PetscLayoutSetSize(A->cmap,N);
1384: PetscLayoutSetLocalSize(A->cmap,n);
1386: PetscLayoutSetUp(A->rmap);
1387: PetscLayoutSetUp(A->cmap);
1389: /* disable operations that are not supported for non-square matrices,
1390: or matrices for which is_row != is_col */
1391: MatHasCongruentLayouts(A,&cong);
1392: if (cong && nr != nc) cong = PETSC_FALSE;
1393: if (cong) {
1394: for (i = 0; cong && i < nr; i++) {
1395: ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1396: }
1397: }
1398: if (!cong) {
1399: A->ops->missingdiagonal = NULL;
1400: A->ops->getdiagonal = NULL;
1401: A->ops->shift = NULL;
1402: A->ops->diagonalset = NULL;
1403: }
1405: PetscCalloc2(nr,&s->left,nc,&s->right);
1406: PetscObjectStateIncrease((PetscObject)A);
1407: A->nonzerostate++;
1408: return(0);
1409: }
1411: /*@
1412: MatNestSetSubMats - Sets the nested submatrices
1414: Collective on Mat
1416: Input Parameter:
1417: + A - nested matrix
1418: . nr - number of nested row blocks
1419: . is_row - index sets for each nested row block, or NULL to make contiguous
1420: . nc - number of nested column blocks
1421: . is_col - index sets for each nested column block, or NULL to make contiguous
1422: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1424: Notes: this always resets any submatrix information previously set
1426: Level: advanced
1428: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1429: @*/
1430: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1431: {
1433: PetscInt i;
1437: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1438: if (nr && is_row) {
1441: }
1442: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1443: if (nc && is_col) {
1446: }
1448: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1449: return(0);
1450: }
1452: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1453: {
1455: PetscBool flg;
1456: PetscInt i,j,m,mi,*ix;
1459: *ltog = NULL;
1460: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1461: if (islocal[i]) {
1462: ISGetLocalSize(islocal[i],&mi);
1463: flg = PETSC_TRUE; /* We found a non-trivial entry */
1464: } else {
1465: ISGetLocalSize(isglobal[i],&mi);
1466: }
1467: m += mi;
1468: }
1469: if (!flg) return(0);
1471: PetscMalloc1(m,&ix);
1472: for (i=0,m=0; i<n; i++) {
1473: ISLocalToGlobalMapping smap = NULL;
1474: Mat sub = NULL;
1475: PetscSF sf;
1476: PetscLayout map;
1477: const PetscInt *ix2;
1479: if (!colflg) {
1480: MatNestFindNonzeroSubMatRow(A,i,&sub);
1481: } else {
1482: MatNestFindNonzeroSubMatCol(A,i,&sub);
1483: }
1484: if (sub) {
1485: if (!colflg) {
1486: MatGetLocalToGlobalMapping(sub,&smap,NULL);
1487: } else {
1488: MatGetLocalToGlobalMapping(sub,NULL,&smap);
1489: }
1490: }
1491: /*
1492: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1493: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1494: */
1495: ISGetIndices(isglobal[i],&ix2);
1496: if (islocal[i]) {
1497: PetscInt *ilocal,*iremote;
1498: PetscInt mil,nleaves;
1500: ISGetLocalSize(islocal[i],&mi);
1501: if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1502: for (j=0; j<mi; j++) ix[m+j] = j;
1503: ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);
1505: /* PetscSFSetGraphLayout does not like negative indices */
1506: PetscMalloc2(mi,&ilocal,mi,&iremote);
1507: for (j=0, nleaves = 0; j<mi; j++) {
1508: if (ix[m+j] < 0) continue;
1509: ilocal[nleaves] = j;
1510: iremote[nleaves] = ix[m+j];
1511: nleaves++;
1512: }
1513: ISGetLocalSize(isglobal[i],&mil);
1514: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1515: PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1516: PetscLayoutSetLocalSize(map,mil);
1517: PetscLayoutSetUp(map);
1518: PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1519: PetscLayoutDestroy(&map);
1520: PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1521: PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1522: PetscSFDestroy(&sf);
1523: PetscFree2(ilocal,iremote);
1524: } else {
1525: ISGetLocalSize(isglobal[i],&mi);
1526: for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1527: }
1528: ISRestoreIndices(isglobal[i],&ix2);
1529: m += mi;
1530: }
1531: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1532: return(0);
1533: }
1536: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1537: /*
1538: nprocessors = NP
1539: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1540: proc 0: => (g_0,h_0,)
1541: proc 1: => (g_1,h_1,)
1542: ...
1543: proc nprocs-1: => (g_NP-1,h_NP-1,)
1545: proc 0: proc 1: proc nprocs-1:
1546: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1548: proc 0:
1549: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1550: proc 1:
1551: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1553: proc NP-1:
1554: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1555: */
1556: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1557: {
1558: Mat_Nest *vs = (Mat_Nest*)A->data;
1559: PetscInt i,j,offset,n,nsum,bs;
1561: Mat sub = NULL;
1564: PetscMalloc1(nr,&vs->isglobal.row);
1565: PetscMalloc1(nc,&vs->isglobal.col);
1566: if (is_row) { /* valid IS is passed in */
1567: /* refs on is[] are incremeneted */
1568: for (i=0; i<vs->nr; i++) {
1569: PetscObjectReference((PetscObject)is_row[i]);
1571: vs->isglobal.row[i] = is_row[i];
1572: }
1573: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1574: nsum = 0;
1575: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1576: MatNestFindNonzeroSubMatRow(A,i,&sub);
1577: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1578: MatGetLocalSize(sub,&n,NULL);
1579: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1580: nsum += n;
1581: }
1582: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1583: offset -= nsum;
1584: for (i=0; i<vs->nr; i++) {
1585: MatNestFindNonzeroSubMatRow(A,i,&sub);
1586: MatGetLocalSize(sub,&n,NULL);
1587: MatGetBlockSizes(sub,&bs,NULL);
1588: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1589: ISSetBlockSize(vs->isglobal.row[i],bs);
1590: offset += n;
1591: }
1592: }
1594: if (is_col) { /* valid IS is passed in */
1595: /* refs on is[] are incremeneted */
1596: for (j=0; j<vs->nc; j++) {
1597: PetscObjectReference((PetscObject)is_col[j]);
1599: vs->isglobal.col[j] = is_col[j];
1600: }
1601: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1602: offset = A->cmap->rstart;
1603: nsum = 0;
1604: for (j=0; j<vs->nc; j++) {
1605: MatNestFindNonzeroSubMatCol(A,j,&sub);
1606: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1607: MatGetLocalSize(sub,NULL,&n);
1608: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1609: nsum += n;
1610: }
1611: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1612: offset -= nsum;
1613: for (j=0; j<vs->nc; j++) {
1614: MatNestFindNonzeroSubMatCol(A,j,&sub);
1615: MatGetLocalSize(sub,NULL,&n);
1616: MatGetBlockSizes(sub,NULL,&bs);
1617: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1618: ISSetBlockSize(vs->isglobal.col[j],bs);
1619: offset += n;
1620: }
1621: }
1623: /* Set up the local ISs */
1624: PetscMalloc1(vs->nr,&vs->islocal.row);
1625: PetscMalloc1(vs->nc,&vs->islocal.col);
1626: for (i=0,offset=0; i<vs->nr; i++) {
1627: IS isloc;
1628: ISLocalToGlobalMapping rmap = NULL;
1629: PetscInt nlocal,bs;
1630: MatNestFindNonzeroSubMatRow(A,i,&sub);
1631: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1632: if (rmap) {
1633: MatGetBlockSizes(sub,&bs,NULL);
1634: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1635: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1636: ISSetBlockSize(isloc,bs);
1637: } else {
1638: nlocal = 0;
1639: isloc = NULL;
1640: }
1641: vs->islocal.row[i] = isloc;
1642: offset += nlocal;
1643: }
1644: for (i=0,offset=0; i<vs->nc; i++) {
1645: IS isloc;
1646: ISLocalToGlobalMapping cmap = NULL;
1647: PetscInt nlocal,bs;
1648: MatNestFindNonzeroSubMatCol(A,i,&sub);
1649: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1650: if (cmap) {
1651: MatGetBlockSizes(sub,NULL,&bs);
1652: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1653: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1654: ISSetBlockSize(isloc,bs);
1655: } else {
1656: nlocal = 0;
1657: isloc = NULL;
1658: }
1659: vs->islocal.col[i] = isloc;
1660: offset += nlocal;
1661: }
1663: /* Set up the aggregate ISLocalToGlobalMapping */
1664: {
1665: ISLocalToGlobalMapping rmap,cmap;
1666: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1667: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1668: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1669: ISLocalToGlobalMappingDestroy(&rmap);
1670: ISLocalToGlobalMappingDestroy(&cmap);
1671: }
1673: #if defined(PETSC_USE_DEBUG)
1674: for (i=0; i<vs->nr; i++) {
1675: for (j=0; j<vs->nc; j++) {
1676: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1677: Mat B = vs->m[i][j];
1678: if (!B) continue;
1679: MatGetSize(B,&M,&N);
1680: MatGetLocalSize(B,&m,&n);
1681: ISGetSize(vs->isglobal.row[i],&Mi);
1682: ISGetSize(vs->isglobal.col[j],&Ni);
1683: ISGetLocalSize(vs->isglobal.row[i],&mi);
1684: ISGetLocalSize(vs->isglobal.col[j],&ni);
1685: 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);
1686: 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);
1687: }
1688: }
1689: #endif
1691: /* Set A->assembled if all non-null blocks are currently assembled */
1692: for (i=0; i<vs->nr; i++) {
1693: for (j=0; j<vs->nc; j++) {
1694: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1695: }
1696: }
1697: A->assembled = PETSC_TRUE;
1698: return(0);
1699: }
1701: /*@C
1702: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1704: Collective on Mat
1706: Input Parameter:
1707: + comm - Communicator for the new Mat
1708: . nr - number of nested row blocks
1709: . is_row - index sets for each nested row block, or NULL to make contiguous
1710: . nc - number of nested column blocks
1711: . is_col - index sets for each nested column block, or NULL to make contiguous
1712: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1714: Output Parameter:
1715: . B - new matrix
1717: Level: advanced
1719: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1720: MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1721: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1722: @*/
1723: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1724: {
1725: Mat A;
1729: *B = 0;
1730: MatCreate(comm,&A);
1731: MatSetType(A,MATNEST);
1732: A->preallocated = PETSC_TRUE;
1733: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1734: *B = A;
1735: return(0);
1736: }
1738: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1739: {
1740: Mat_Nest *nest = (Mat_Nest*)A->data;
1741: Mat *trans;
1742: PetscScalar **avv;
1743: PetscScalar *vv;
1744: PetscInt **aii,**ajj;
1745: PetscInt *ii,*jj,*ci;
1746: PetscInt nr,nc,nnz,i,j;
1747: PetscBool done;
1751: MatGetSize(A,&nr,&nc);
1752: if (reuse == MAT_REUSE_MATRIX) {
1753: PetscInt rnr;
1755: MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1756: if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1757: if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1758: MatSeqAIJGetArray(*newmat,&vv);
1759: }
1760: /* extract CSR for nested SeqAIJ matrices */
1761: nnz = 0;
1762: PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1763: for (i=0; i<nest->nr; ++i) {
1764: for (j=0; j<nest->nc; ++j) {
1765: Mat B = nest->m[i][j];
1766: if (B) {
1767: PetscScalar *naa;
1768: PetscInt *nii,*njj,nnr;
1769: PetscBool istrans;
1771: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1772: if (istrans) {
1773: Mat Bt;
1775: MatTransposeGetMat(B,&Bt);
1776: MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1777: B = trans[i*nest->nc+j];
1778: }
1779: MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1780: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1781: MatSeqAIJGetArray(B,&naa);
1782: nnz += nii[nnr];
1784: aii[i*nest->nc+j] = nii;
1785: ajj[i*nest->nc+j] = njj;
1786: avv[i*nest->nc+j] = naa;
1787: }
1788: }
1789: }
1790: if (reuse != MAT_REUSE_MATRIX) {
1791: PetscMalloc1(nr+1,&ii);
1792: PetscMalloc1(nnz,&jj);
1793: PetscMalloc1(nnz,&vv);
1794: } else {
1795: if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1796: }
1798: /* new row pointer */
1799: PetscArrayzero(ii,nr+1);
1800: for (i=0; i<nest->nr; ++i) {
1801: PetscInt ncr,rst;
1803: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1804: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1805: for (j=0; j<nest->nc; ++j) {
1806: if (aii[i*nest->nc+j]) {
1807: PetscInt *nii = aii[i*nest->nc+j];
1808: PetscInt ir;
1810: for (ir=rst; ir<ncr+rst; ++ir) {
1811: ii[ir+1] += nii[1]-nii[0];
1812: nii++;
1813: }
1814: }
1815: }
1816: }
1817: for (i=0; i<nr; i++) ii[i+1] += ii[i];
1819: /* construct CSR for the new matrix */
1820: PetscCalloc1(nr,&ci);
1821: for (i=0; i<nest->nr; ++i) {
1822: PetscInt ncr,rst;
1824: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1825: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1826: for (j=0; j<nest->nc; ++j) {
1827: if (aii[i*nest->nc+j]) {
1828: PetscScalar *nvv = avv[i*nest->nc+j];
1829: PetscInt *nii = aii[i*nest->nc+j];
1830: PetscInt *njj = ajj[i*nest->nc+j];
1831: PetscInt ir,cst;
1833: ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1834: for (ir=rst; ir<ncr+rst; ++ir) {
1835: PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1837: for (ij=0;ij<rsize;ij++) {
1838: jj[ist+ij] = *njj+cst;
1839: vv[ist+ij] = *nvv;
1840: njj++;
1841: nvv++;
1842: }
1843: ci[ir] += rsize;
1844: nii++;
1845: }
1846: }
1847: }
1848: }
1849: PetscFree(ci);
1851: /* restore info */
1852: for (i=0; i<nest->nr; ++i) {
1853: for (j=0; j<nest->nc; ++j) {
1854: Mat B = nest->m[i][j];
1855: if (B) {
1856: PetscInt nnr = 0, k = i*nest->nc+j;
1858: B = (trans[k] ? trans[k] : B);
1859: MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1860: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1861: MatSeqAIJRestoreArray(B,&avv[k]);
1862: MatDestroy(&trans[k]);
1863: }
1864: }
1865: }
1866: PetscFree4(aii,ajj,avv,trans);
1868: /* finalize newmat */
1869: if (reuse == MAT_INITIAL_MATRIX) {
1870: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1871: } else if (reuse == MAT_INPLACE_MATRIX) {
1872: Mat B;
1874: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1875: MatHeaderReplace(A,&B);
1876: }
1877: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1878: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1879: {
1880: Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1881: a->free_a = PETSC_TRUE;
1882: a->free_ij = PETSC_TRUE;
1883: }
1884: return(0);
1885: }
1887: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1888: {
1890: Mat_Nest *nest = (Mat_Nest*)A->data;
1891: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1892: PetscInt cstart,cend;
1893: PetscMPIInt size;
1894: Mat C;
1897: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1898: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1899: PetscInt nf;
1900: PetscBool fast;
1902: PetscStrcmp(newtype,MATAIJ,&fast);
1903: if (!fast) {
1904: PetscStrcmp(newtype,MATSEQAIJ,&fast);
1905: }
1906: for (i=0; i<nest->nr && fast; ++i) {
1907: for (j=0; j<nest->nc && fast; ++j) {
1908: Mat B = nest->m[i][j];
1909: if (B) {
1910: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1911: if (!fast) {
1912: PetscBool istrans;
1914: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1915: if (istrans) {
1916: Mat Bt;
1918: MatTransposeGetMat(B,&Bt);
1919: PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1920: }
1921: }
1922: }
1923: }
1924: }
1925: for (i=0, nf=0; i<nest->nr && fast; ++i) {
1926: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1927: if (fast) {
1928: PetscInt f,s;
1930: ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1931: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1932: else {
1933: ISGetSize(nest->isglobal.row[i],&f);
1934: nf += f;
1935: }
1936: }
1937: }
1938: for (i=0, nf=0; i<nest->nc && fast; ++i) {
1939: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1940: if (fast) {
1941: PetscInt f,s;
1943: ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1944: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1945: else {
1946: ISGetSize(nest->isglobal.col[i],&f);
1947: nf += f;
1948: }
1949: }
1950: }
1951: if (fast) {
1952: MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1953: return(0);
1954: }
1955: }
1956: MatGetSize(A,&M,&N);
1957: MatGetLocalSize(A,&m,&n);
1958: MatGetOwnershipRangeColumn(A,&cstart,&cend);
1959: switch (reuse) {
1960: case MAT_INITIAL_MATRIX:
1961: MatCreate(PetscObjectComm((PetscObject)A),&C);
1962: MatSetType(C,newtype);
1963: MatSetSizes(C,m,n,M,N);
1964: *newmat = C;
1965: break;
1966: case MAT_REUSE_MATRIX:
1967: C = *newmat;
1968: break;
1969: default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1970: }
1971: PetscMalloc1(2*m,&dnnz);
1972: onnz = dnnz + m;
1973: for (k=0; k<m; k++) {
1974: dnnz[k] = 0;
1975: onnz[k] = 0;
1976: }
1977: for (j=0; j<nest->nc; ++j) {
1978: IS bNis;
1979: PetscInt bN;
1980: const PetscInt *bNindices;
1981: /* Using global column indices and ISAllGather() is not scalable. */
1982: ISAllGather(nest->isglobal.col[j], &bNis);
1983: ISGetSize(bNis, &bN);
1984: ISGetIndices(bNis,&bNindices);
1985: for (i=0; i<nest->nr; ++i) {
1986: PetscSF bmsf;
1987: PetscSFNode *iremote;
1988: Mat B;
1989: PetscInt bm, *sub_dnnz,*sub_onnz, br;
1990: const PetscInt *bmindices;
1991: B = nest->m[i][j];
1992: if (!B) continue;
1993: ISGetLocalSize(nest->isglobal.row[i],&bm);
1994: ISGetIndices(nest->isglobal.row[i],&bmindices);
1995: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1996: PetscMalloc1(bm,&iremote);
1997: PetscMalloc1(bm,&sub_dnnz);
1998: PetscMalloc1(bm,&sub_onnz);
1999: for (k = 0; k < bm; ++k){
2000: sub_dnnz[k] = 0;
2001: sub_onnz[k] = 0;
2002: }
2003: /*
2004: Locate the owners for all of the locally-owned global row indices for this row block.
2005: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2006: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2007: */
2008: MatGetOwnershipRange(B,&rstart,NULL);
2009: for (br = 0; br < bm; ++br) {
2010: PetscInt row = bmindices[br], brncols, col;
2011: const PetscInt *brcols;
2012: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2013: PetscMPIInt rowowner = 0;
2014: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2015: /* how many roots */
2016: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2017: /* get nonzero pattern */
2018: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2019: for (k=0; k<brncols; k++) {
2020: col = bNindices[brcols[k]];
2021: if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2022: sub_dnnz[br]++;
2023: } else {
2024: sub_onnz[br]++;
2025: }
2026: }
2027: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2028: }
2029: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2030: /* bsf will have to take care of disposing of bedges. */
2031: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2032: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2033: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2034: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2035: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2036: PetscFree(sub_dnnz);
2037: PetscFree(sub_onnz);
2038: PetscSFDestroy(&bmsf);
2039: }
2040: ISRestoreIndices(bNis,&bNindices);
2041: ISDestroy(&bNis);
2042: }
2043: /* Resize preallocation if overestimated */
2044: for (i=0;i<m;i++) {
2045: dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2046: onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2047: }
2048: MatSeqAIJSetPreallocation(C,0,dnnz);
2049: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2050: PetscFree(dnnz);
2052: /* Fill by row */
2053: for (j=0; j<nest->nc; ++j) {
2054: /* Using global column indices and ISAllGather() is not scalable. */
2055: IS bNis;
2056: PetscInt bN;
2057: const PetscInt *bNindices;
2058: ISAllGather(nest->isglobal.col[j], &bNis);
2059: ISGetSize(bNis,&bN);
2060: ISGetIndices(bNis,&bNindices);
2061: for (i=0; i<nest->nr; ++i) {
2062: Mat B;
2063: PetscInt bm, br;
2064: const PetscInt *bmindices;
2065: B = nest->m[i][j];
2066: if (!B) continue;
2067: ISGetLocalSize(nest->isglobal.row[i],&bm);
2068: ISGetIndices(nest->isglobal.row[i],&bmindices);
2069: MatGetOwnershipRange(B,&rstart,NULL);
2070: for (br = 0; br < bm; ++br) {
2071: PetscInt row = bmindices[br], brncols, *cols;
2072: const PetscInt *brcols;
2073: const PetscScalar *brcoldata;
2074: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2075: PetscMalloc1(brncols,&cols);
2076: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2077: /*
2078: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2079: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2080: */
2081: MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2082: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2083: PetscFree(cols);
2084: }
2085: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2086: }
2087: ISRestoreIndices(bNis,&bNindices);
2088: ISDestroy(&bNis);
2089: }
2090: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2091: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2092: return(0);
2093: }
2095: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2096: {
2097: Mat_Nest *bA = (Mat_Nest*)mat->data;
2098: MatOperation opAdd;
2099: PetscInt i,j,nr = bA->nr,nc = bA->nc;
2100: PetscBool flg;
2104: *has = PETSC_FALSE;
2105: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2106: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2107: for (j=0; j<nc; j++) {
2108: for (i=0; i<nr; i++) {
2109: if (!bA->m[i][j]) continue;
2110: MatHasOperation(bA->m[i][j],opAdd,&flg);
2111: if (!flg) return(0);
2112: }
2113: }
2114: }
2115: if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2116: return(0);
2117: }
2119: /*MC
2120: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2122: Level: intermediate
2124: Notes:
2125: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2126: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2127: It is usually used with DMComposite and DMCreateMatrix()
2129: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2130: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2131: than the nest matrix.
2133: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2134: VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2135: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2136: M*/
2137: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2138: {
2139: Mat_Nest *s;
2143: PetscNewLog(A,&s);
2144: A->data = (void*)s;
2146: s->nr = -1;
2147: s->nc = -1;
2148: s->m = NULL;
2149: s->splitassembly = PETSC_FALSE;
2151: PetscMemzero(A->ops,sizeof(*A->ops));
2153: A->ops->mult = MatMult_Nest;
2154: A->ops->multadd = MatMultAdd_Nest;
2155: A->ops->multtranspose = MatMultTranspose_Nest;
2156: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2157: A->ops->transpose = MatTranspose_Nest;
2158: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2159: A->ops->assemblyend = MatAssemblyEnd_Nest;
2160: A->ops->zeroentries = MatZeroEntries_Nest;
2161: A->ops->copy = MatCopy_Nest;
2162: A->ops->axpy = MatAXPY_Nest;
2163: A->ops->duplicate = MatDuplicate_Nest;
2164: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2165: A->ops->destroy = MatDestroy_Nest;
2166: A->ops->view = MatView_Nest;
2167: A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2168: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2169: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2170: A->ops->getdiagonal = MatGetDiagonal_Nest;
2171: A->ops->diagonalscale = MatDiagonalScale_Nest;
2172: A->ops->scale = MatScale_Nest;
2173: A->ops->shift = MatShift_Nest;
2174: A->ops->diagonalset = MatDiagonalSet_Nest;
2175: A->ops->setrandom = MatSetRandom_Nest;
2176: A->ops->hasoperation = MatHasOperation_Nest;
2177: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2179: A->spptr = 0;
2180: A->assembled = PETSC_FALSE;
2182: /* expose Nest api's */
2183: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2184: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2185: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2186: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
2187: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
2188: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2189: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
2190: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2191: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2192: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2193: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2194: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);
2195: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2196: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2197: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);
2199: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2200: return(0);
2201: }