Actual source code: matnest.c
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;
17: *m = *n = *M = *N = 0;
18: for (i=0; i<bA->nr; i++) { /* rows */
19: PetscInt sm,sM;
20: ISGetLocalSize(bA->isglobal.row[i],&sm);
21: ISGetSize(bA->isglobal.row[i],&sM);
22: *m += sm;
23: *M += sM;
24: }
25: for (j=0; j<bA->nc; j++) { /* cols */
26: PetscInt sn,sN;
27: ISGetLocalSize(bA->isglobal.col[j],&sn);
28: ISGetSize(bA->isglobal.col[j],&sN);
29: *n += sn;
30: *N += sN;
31: }
32: return 0;
33: }
35: /* operations */
36: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
37: {
38: Mat_Nest *bA = (Mat_Nest*)A->data;
39: Vec *bx = bA->right,*by = bA->left;
40: PetscInt i,j,nr = bA->nr,nc = bA->nc;
42: for (i=0; i<nr; i++) VecGetSubVector(y,bA->isglobal.row[i],&by[i]);
43: for (i=0; i<nc; i++) VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);
44: for (i=0; i<nr; i++) {
45: VecZeroEntries(by[i]);
46: for (j=0; j<nc; j++) {
47: if (!bA->m[i][j]) continue;
48: /* y[i] <- y[i] + A[i][j] * x[j] */
49: MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);
50: }
51: }
52: for (i=0; i<nr; i++) VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);
53: for (i=0; i<nc; i++) VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);
54: return 0;
55: }
57: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
58: {
59: Mat_Nest *bA = (Mat_Nest*)A->data;
60: Vec *bx = bA->right,*bz = bA->left;
61: PetscInt i,j,nr = bA->nr,nc = bA->nc;
63: for (i=0; i<nr; i++) VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);
64: for (i=0; i<nc; i++) VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);
65: for (i=0; i<nr; i++) {
66: if (y != z) {
67: Vec by;
68: VecGetSubVector(y,bA->isglobal.row[i],&by);
69: VecCopy(by,bz[i]);
70: VecRestoreSubVector(y,bA->isglobal.row[i],&by);
71: }
72: for (j=0; j<nc; j++) {
73: if (!bA->m[i][j]) continue;
74: /* y[i] <- y[i] + A[i][j] * x[j] */
75: MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
76: }
77: }
78: for (i=0; i<nr; i++) VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);
79: for (i=0; i<nc; i++) VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);
80: return 0;
81: }
83: typedef struct {
84: Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
85: PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */
86: PetscInt *dm,*dn,k; /* displacements and number of submatrices */
87: } Nest_Dense;
89: PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
90: {
91: Mat_Nest *bA;
92: Nest_Dense *contents;
93: Mat viewB,viewC,productB,workC;
94: const PetscScalar *barray;
95: PetscScalar *carray;
96: PetscInt i,j,M,N,nr,nc,ldb,ldc;
97: Mat A,B;
99: MatCheckProduct(C,3);
100: A = C->product->A;
101: B = C->product->B;
102: MatGetSize(B,NULL,&N);
103: if (!N) {
104: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
105: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
106: return 0;
107: }
108: contents = (Nest_Dense*)C->product->data;
110: bA = (Mat_Nest*)A->data;
111: nr = bA->nr;
112: nc = bA->nc;
113: MatDenseGetLDA(B,&ldb);
114: MatDenseGetLDA(C,&ldc);
115: MatZeroEntries(C);
116: MatDenseGetArrayRead(B,&barray);
117: MatDenseGetArray(C,&carray);
118: for (i=0; i<nr; i++) {
119: ISGetSize(bA->isglobal.row[i],&M);
120: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
121: MatDenseSetLDA(viewC,ldc);
122: for (j=0; j<nc; j++) {
123: if (!bA->m[i][j]) continue;
124: ISGetSize(bA->isglobal.col[j],&M);
125: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
126: MatDenseSetLDA(viewB,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: MatProductNumeric(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;
154: PetscFree(contents->tarray);
155: for (i=0; i<contents->k; i++) {
156: MatDestroy(contents->workC + i);
157: }
158: PetscFree3(contents->dm,contents->dn,contents->workC);
159: PetscFree(contents);
160: return 0;
161: }
163: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
164: {
165: Mat_Nest *bA;
166: Mat viewB,workC;
167: const PetscScalar *barray;
168: PetscInt i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
169: Nest_Dense *contents=NULL;
170: PetscBool cisdense;
171: Mat A,B;
172: PetscReal fill;
174: MatCheckProduct(C,4);
176: A = C->product->A;
177: B = C->product->B;
178: fill = C->product->fill;
179: bA = (Mat_Nest*)A->data;
180: nr = bA->nr;
181: nc = bA->nc;
182: MatGetLocalSize(C,&m,&n);
183: MatGetSize(C,&M,&N);
184: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
185: MatGetLocalSize(B,NULL,&n);
186: MatGetSize(B,NULL,&N);
187: MatGetLocalSize(A,&m,NULL);
188: MatGetSize(A,&M,NULL);
189: MatSetSizes(C,m,n,M,N);
190: }
191: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"");
192: if (!cisdense) {
193: MatSetType(C,((PetscObject)B)->type_name);
194: }
195: MatSetUp(C);
196: if (!N) {
197: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
198: return 0;
199: }
201: PetscNew(&contents);
202: C->product->data = contents;
203: C->product->destroy = MatNest_DenseDestroy;
204: PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
205: contents->k = nr*nc;
206: for (i=0; i<nr; i++) {
207: ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
208: maxm = PetscMax(maxm,contents->dm[i+1]);
209: contents->dm[i+1] += contents->dm[i];
210: }
211: for (i=0; i<nc; i++) {
212: ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
213: contents->dn[i+1] += contents->dn[i];
214: }
215: PetscMalloc1(maxm*N,&contents->tarray);
216: MatDenseGetLDA(B,&ldb);
217: MatGetSize(B,NULL,&N);
218: MatDenseGetArrayRead(B,&barray);
219: /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
220: for (j=0; j<nc; j++) {
221: ISGetSize(bA->isglobal.col[j],&M);
222: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
223: MatDenseSetLDA(viewB,ldb);
224: for (i=0; i<nr; i++) {
225: if (!bA->m[i][j]) continue;
226: /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
228: MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
229: workC = contents->workC[i*nc + j];
230: MatProductSetType(workC,MATPRODUCT_AB);
231: MatProductSetAlgorithm(workC,"default");
232: MatProductSetFill(workC,fill);
233: MatProductSetFromOptions(workC);
234: MatProductSymbolic(workC);
236: /* since tarray will be shared by all Mat */
237: MatSeqDenseSetPreallocation(workC,contents->tarray);
238: MatMPIDenseSetPreallocation(workC,contents->tarray);
239: }
240: MatDestroy(&viewB);
241: }
242: MatDenseRestoreArrayRead(B,&barray);
244: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
245: return 0;
246: }
248: /* --------------------------------------------------------- */
249: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
250: {
251: C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
252: return 0;
253: }
255: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
256: {
257: Mat_Product *product = C->product;
259: if (product->type == MATPRODUCT_AB) {
260: MatProductSetFromOptions_Nest_Dense_AB(C);
261: }
262: return 0;
263: }
264: /* --------------------------------------------------------- */
266: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
267: {
268: Mat_Nest *bA = (Mat_Nest*)A->data;
269: Vec *bx = bA->left,*by = bA->right;
270: PetscInt i,j,nr = bA->nr,nc = bA->nc;
272: for (i=0; i<nr; i++) VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);
273: for (i=0; i<nc; i++) VecGetSubVector(y,bA->isglobal.col[i],&by[i]);
274: for (j=0; j<nc; j++) {
275: VecZeroEntries(by[j]);
276: for (i=0; i<nr; i++) {
277: if (!bA->m[i][j]) continue;
278: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
279: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
280: }
281: }
282: for (i=0; i<nr; i++) VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);
283: for (i=0; i<nc; i++) VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);
284: return 0;
285: }
287: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
288: {
289: Mat_Nest *bA = (Mat_Nest*)A->data;
290: Vec *bx = bA->left,*bz = bA->right;
291: PetscInt i,j,nr = bA->nr,nc = bA->nc;
293: for (i=0; i<nr; i++) VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);
294: for (i=0; i<nc; i++) VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);
295: for (j=0; j<nc; j++) {
296: if (y != z) {
297: Vec by;
298: VecGetSubVector(y,bA->isglobal.col[j],&by);
299: VecCopy(by,bz[j]);
300: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
301: }
302: for (i=0; i<nr; i++) {
303: if (!bA->m[i][j]) continue;
304: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
305: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
306: }
307: }
308: for (i=0; i<nr; i++) VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);
309: for (i=0; i<nc; i++) VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);
310: return 0;
311: }
313: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
314: {
315: Mat_Nest *bA = (Mat_Nest*)A->data, *bC;
316: Mat C;
317: PetscInt i,j,nr = bA->nr,nc = bA->nc;
321: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
322: Mat *subs;
323: IS *is_row,*is_col;
325: PetscCalloc1(nr * nc,&subs);
326: PetscMalloc2(nr,&is_row,nc,&is_col);
327: MatNestGetISs(A,is_row,is_col);
328: if (reuse == MAT_INPLACE_MATRIX) {
329: for (i=0; i<nr; i++) {
330: for (j=0; j<nc; j++) {
331: subs[i + nr * j] = bA->m[i][j];
332: }
333: }
334: }
336: MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
337: PetscFree(subs);
338: PetscFree2(is_row,is_col);
339: } else {
340: C = *B;
341: }
343: bC = (Mat_Nest*)C->data;
344: for (i=0; i<nr; i++) {
345: for (j=0; j<nc; j++) {
346: if (bA->m[i][j]) {
347: MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
348: } else {
349: bC->m[j][i] = NULL;
350: }
351: }
352: }
354: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
355: *B = C;
356: } else {
357: MatHeaderMerge(A, &C);
358: }
359: return 0;
360: }
362: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
363: {
364: IS *lst = *list;
365: PetscInt i;
367: if (!lst) return 0;
368: for (i=0; i<n; i++) if (lst[i]) ISDestroy(&lst[i]);
369: PetscFree(lst);
370: *list = NULL;
371: return 0;
372: }
374: static PetscErrorCode MatReset_Nest(Mat A)
375: {
376: Mat_Nest *vs = (Mat_Nest*)A->data;
377: PetscInt i,j;
379: /* release the matrices and the place holders */
380: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
381: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
382: MatNestDestroyISList(vs->nr,&vs->islocal.row);
383: MatNestDestroyISList(vs->nc,&vs->islocal.col);
385: PetscFree(vs->row_len);
386: PetscFree(vs->col_len);
387: PetscFree(vs->nnzstate);
389: PetscFree2(vs->left,vs->right);
391: /* release the matrices and the place holders */
392: if (vs->m) {
393: for (i=0; i<vs->nr; i++) {
394: for (j=0; j<vs->nc; j++) {
395: MatDestroy(&vs->m[i][j]);
396: }
397: PetscFree(vs->m[i]);
398: }
399: PetscFree(vs->m);
400: }
402: /* restore defaults */
403: vs->nr = 0;
404: vs->nc = 0;
405: vs->splitassembly = PETSC_FALSE;
406: return 0;
407: }
409: static PetscErrorCode MatDestroy_Nest(Mat A)
410: {
411: MatReset_Nest(A);
412: PetscFree(A->data);
413: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",NULL);
414: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",NULL);
415: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",NULL);
416: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",NULL);
417: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",NULL);
418: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",NULL);
419: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",NULL);
420: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",NULL);
421: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",NULL);
422: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",NULL);
423: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",NULL);
424: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",NULL);
425: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",NULL);
426: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",NULL);
427: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
428: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
429: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
430: return 0;
431: }
433: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
434: {
435: Mat_Nest *vs = (Mat_Nest*)mat->data;
436: PetscInt i;
438: if (dd) *dd = 0;
439: if (!vs->nr) {
440: *missing = PETSC_TRUE;
441: return 0;
442: }
443: *missing = PETSC_FALSE;
444: for (i = 0; i < vs->nr && !(*missing); i++) {
445: *missing = PETSC_TRUE;
446: if (vs->m[i][i]) {
447: MatMissingDiagonal(vs->m[i][i],missing,NULL);
449: }
450: }
451: return 0;
452: }
454: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
455: {
456: Mat_Nest *vs = (Mat_Nest*)A->data;
457: PetscInt i,j;
458: PetscBool nnzstate = PETSC_FALSE;
460: for (i=0; i<vs->nr; i++) {
461: for (j=0; j<vs->nc; j++) {
462: PetscObjectState subnnzstate = 0;
463: if (vs->m[i][j]) {
464: MatAssemblyBegin(vs->m[i][j],type);
465: if (!vs->splitassembly) {
466: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
467: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
468: * already performing an assembly, but the result would by more complicated and appears to offer less
469: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
470: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
471: */
472: MatAssemblyEnd(vs->m[i][j],type);
473: MatGetNonzeroState(vs->m[i][j],&subnnzstate);
474: }
475: }
476: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
477: vs->nnzstate[i*vs->nc+j] = subnnzstate;
478: }
479: }
480: if (nnzstate) A->nonzerostate++;
481: return 0;
482: }
484: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
485: {
486: Mat_Nest *vs = (Mat_Nest*)A->data;
487: PetscInt i,j;
489: for (i=0; i<vs->nr; i++) {
490: for (j=0; j<vs->nc; j++) {
491: if (vs->m[i][j]) {
492: if (vs->splitassembly) {
493: MatAssemblyEnd(vs->m[i][j],type);
494: }
495: }
496: }
497: }
498: return 0;
499: }
501: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
502: {
503: Mat_Nest *vs = (Mat_Nest*)A->data;
504: PetscInt j;
505: Mat sub;
507: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
508: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
509: if (sub) MatSetUp(sub); /* Ensure that the sizes are available */
510: *B = sub;
511: return 0;
512: }
514: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
515: {
516: Mat_Nest *vs = (Mat_Nest*)A->data;
517: PetscInt i;
518: Mat sub;
520: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
521: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
522: if (sub) MatSetUp(sub); /* Ensure that the sizes are available */
523: *B = sub;
524: return 0;
525: }
527: static PetscErrorCode MatNestFindISRange(Mat A,PetscInt n,const IS list[],IS is,PetscInt *begin,PetscInt *end)
528: {
529: PetscInt i,j,size,m;
530: PetscBool flg;
531: IS out,concatenate[2];
535: if (begin) {
537: *begin = -1;
538: }
539: if (end) {
541: *end = -1;
542: }
543: for (i=0; i<n; i++) {
544: if (!list[i]) continue;
545: ISEqualUnsorted(list[i],is,&flg);
546: if (flg) {
547: if (begin) *begin = i;
548: if (end) *end = i+1;
549: return 0;
550: }
551: }
552: ISGetSize(is,&size);
553: for (i=0; i<n-1; i++) {
554: if (!list[i]) continue;
555: m = 0;
556: ISConcatenate(PetscObjectComm((PetscObject)A),2,list+i,&out);
557: ISGetSize(out,&m);
558: for (j=i+2; j<n && m<size; j++) {
559: if (list[j]) {
560: concatenate[0] = out;
561: concatenate[1] = list[j];
562: ISConcatenate(PetscObjectComm((PetscObject)A),2,concatenate,&out);
563: ISDestroy(concatenate);
564: ISGetSize(out,&m);
565: }
566: }
567: if (m == size) {
568: ISEqualUnsorted(out,is,&flg);
569: if (flg) {
570: if (begin) *begin = i;
571: if (end) *end = j;
572: ISDestroy(&out);
573: return 0;
574: }
575: }
576: ISDestroy(&out);
577: }
578: return 0;
579: }
581: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A,PetscInt i,PetscInt j,Mat *B)
582: {
583: Mat_Nest *vs = (Mat_Nest*)A->data;
584: PetscInt lr,lc;
586: MatCreate(PetscObjectComm((PetscObject)A),B);
587: ISGetLocalSize(vs->isglobal.row[i],&lr);
588: ISGetLocalSize(vs->isglobal.col[j],&lc);
589: MatSetSizes(*B,lr,lc,PETSC_DECIDE,PETSC_DECIDE);
590: MatSetType(*B,MATAIJ);
591: MatSeqAIJSetPreallocation(*B,0,NULL);
592: MatMPIAIJSetPreallocation(*B,0,NULL,0,NULL);
593: MatSetUp(*B);
594: MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
595: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
596: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
597: return 0;
598: }
600: static PetscErrorCode MatNestGetBlock_Private(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *B)
601: {
602: Mat_Nest *vs = (Mat_Nest*)A->data;
603: Mat *a;
604: PetscInt i,j,k,l,nr=rend-rbegin,nc=cend-cbegin;
605: char keyname[256];
606: PetscBool *b;
607: PetscBool flg;
609: *B = NULL;
610: PetscSNPrintf(keyname,sizeof(keyname),"NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT,rbegin,rend,cbegin,cend);
611: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
612: if (*B) return 0;
614: PetscMalloc2(nr*nc,&a,nr*nc,&b);
615: for (i=0; i<nr; i++) {
616: for (j=0; j<nc; j++) {
617: a[i*nc + j] = vs->m[rbegin+i][cbegin+j];
618: b[i*nc + j] = PETSC_FALSE;
619: }
620: }
621: if (nc!=vs->nc&&nr!=vs->nr) {
622: for (i=0; i<nr; i++) {
623: for (j=0; j<nc; j++) {
624: flg = PETSC_FALSE;
625: for (k=0; (k<nr&&!flg); k++) {
626: if (a[j + k*nc]) flg = PETSC_TRUE;
627: }
628: if (flg) {
629: flg = PETSC_FALSE;
630: for (l=0; (l<nc&&!flg); l++) {
631: if (a[i*nc + l]) flg = PETSC_TRUE;
632: }
633: }
634: if (!flg) {
635: b[i*nc + j] = PETSC_TRUE;
636: MatNestFillEmptyMat_Private(A,rbegin+i,cbegin+j,a + i*nc + j);
637: }
638: }
639: }
640: }
641: MatCreateNest(PetscObjectComm((PetscObject)A),nr,nr!=vs->nr?NULL:vs->isglobal.row,nc,nc!=vs->nc?NULL:vs->isglobal.col,a,B);
642: for (i=0; i<nr; i++) {
643: for (j=0; j<nc; j++) {
644: if (b[i*nc + j]) {
645: MatDestroy(a + i*nc + j);
646: }
647: }
648: }
649: PetscFree2(a,b);
650: (*B)->assembled = A->assembled;
651: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
652: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
653: return 0;
654: }
656: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
657: {
658: Mat_Nest *vs = (Mat_Nest*)A->data;
659: PetscInt rbegin,rend,cbegin,cend;
661: MatNestFindISRange(A,vs->nr,is->row,isrow,&rbegin,&rend);
662: MatNestFindISRange(A,vs->nc,is->col,iscol,&cbegin,&cend);
663: if (rend == rbegin + 1 && cend == cbegin + 1) {
664: if (!vs->m[rbegin][cbegin]) {
665: MatNestFillEmptyMat_Private(A,rbegin,cbegin,vs->m[rbegin] + cbegin);
666: }
667: *B = vs->m[rbegin][cbegin];
668: } else if (rbegin != -1 && cbegin != -1) {
669: MatNestGetBlock_Private(A,rbegin,rend,cbegin,cend,B);
670: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
671: return 0;
672: }
674: /*
675: TODO: This does not actually returns a submatrix we can modify
676: */
677: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
678: {
679: Mat_Nest *vs = (Mat_Nest*)A->data;
680: Mat sub;
682: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
683: switch (reuse) {
684: case MAT_INITIAL_MATRIX:
685: if (sub) PetscObjectReference((PetscObject)sub);
686: *B = sub;
687: break;
688: case MAT_REUSE_MATRIX:
690: break;
691: case MAT_IGNORE_MATRIX: /* Nothing to do */
692: break;
693: case MAT_INPLACE_MATRIX: /* Nothing to do */
694: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
695: }
696: return 0;
697: }
699: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
700: {
701: Mat_Nest *vs = (Mat_Nest*)A->data;
702: Mat sub;
704: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
705: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
706: if (sub) PetscObjectReference((PetscObject)sub);
707: *B = sub;
708: return 0;
709: }
711: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
712: {
713: Mat_Nest *vs = (Mat_Nest*)A->data;
714: Mat sub;
716: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
718: if (sub) {
720: MatDestroy(B);
721: }
722: return 0;
723: }
725: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
726: {
727: Mat_Nest *bA = (Mat_Nest*)A->data;
728: PetscInt i;
730: for (i=0; i<bA->nr; i++) {
731: Vec bv;
732: VecGetSubVector(v,bA->isglobal.row[i],&bv);
733: if (bA->m[i][i]) {
734: MatGetDiagonal(bA->m[i][i],bv);
735: } else {
736: VecSet(bv,0.0);
737: }
738: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
739: }
740: return 0;
741: }
743: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
744: {
745: Mat_Nest *bA = (Mat_Nest*)A->data;
746: Vec bl,*br;
747: PetscInt i,j;
749: PetscCalloc1(bA->nc,&br);
750: if (r) {
751: for (j=0; j<bA->nc; j++) VecGetSubVector(r,bA->isglobal.col[j],&br[j]);
752: }
753: bl = NULL;
754: for (i=0; i<bA->nr; i++) {
755: if (l) {
756: VecGetSubVector(l,bA->isglobal.row[i],&bl);
757: }
758: for (j=0; j<bA->nc; j++) {
759: if (bA->m[i][j]) {
760: MatDiagonalScale(bA->m[i][j],bl,br[j]);
761: }
762: }
763: if (l) {
764: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
765: }
766: }
767: if (r) {
768: for (j=0; j<bA->nc; j++) VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);
769: }
770: PetscFree(br);
771: return 0;
772: }
774: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
775: {
776: Mat_Nest *bA = (Mat_Nest*)A->data;
777: PetscInt i,j;
779: for (i=0; i<bA->nr; i++) {
780: for (j=0; j<bA->nc; j++) {
781: if (bA->m[i][j]) {
782: MatScale(bA->m[i][j],a);
783: }
784: }
785: }
786: return 0;
787: }
789: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
790: {
791: Mat_Nest *bA = (Mat_Nest*)A->data;
792: PetscInt i;
793: PetscBool nnzstate = PETSC_FALSE;
795: for (i=0; i<bA->nr; i++) {
796: PetscObjectState subnnzstate = 0;
798: MatShift(bA->m[i][i],a);
799: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
800: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
801: bA->nnzstate[i*bA->nc+i] = subnnzstate;
802: }
803: if (nnzstate) A->nonzerostate++;
804: return 0;
805: }
807: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
808: {
809: Mat_Nest *bA = (Mat_Nest*)A->data;
810: PetscInt i;
811: 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;
834: for (i=0; i<bA->nr; i++) {
835: for (j=0; j<bA->nc; j++) {
836: if (bA->m[i][j]) {
837: MatSetRandom(bA->m[i][j],rctx);
838: }
839: }
840: }
841: return 0;
842: }
844: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
845: {
846: Mat_Nest *bA = (Mat_Nest*)A->data;
847: Vec *L,*R;
848: MPI_Comm comm;
849: PetscInt i,j;
851: PetscObjectGetComm((PetscObject)A,&comm);
852: if (right) {
853: /* allocate R */
854: PetscMalloc1(bA->nc, &R);
855: /* Create the right vectors */
856: for (j=0; j<bA->nc; j++) {
857: for (i=0; i<bA->nr; i++) {
858: if (bA->m[i][j]) {
859: MatCreateVecs(bA->m[i][j],&R[j],NULL);
860: break;
861: }
862: }
864: }
865: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
866: /* hand back control to the nest vector */
867: for (j=0; j<bA->nc; j++) {
868: VecDestroy(&R[j]);
869: }
870: PetscFree(R);
871: }
873: if (left) {
874: /* allocate L */
875: PetscMalloc1(bA->nr, &L);
876: /* Create the left vectors */
877: for (i=0; i<bA->nr; i++) {
878: for (j=0; j<bA->nc; j++) {
879: if (bA->m[i][j]) {
880: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
881: break;
882: }
883: }
885: }
887: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
888: for (i=0; i<bA->nr; i++) {
889: VecDestroy(&L[i]);
890: }
892: PetscFree(L);
893: }
894: return 0;
895: }
897: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
898: {
899: Mat_Nest *bA = (Mat_Nest*)A->data;
900: PetscBool isascii,viewSub = PETSC_FALSE;
901: PetscInt i,j;
903: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
904: if (isascii) {
906: PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
907: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
908: PetscViewerASCIIPushTab(viewer);
909: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n",bA->nr,bA->nc);
911: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
912: for (i=0; i<bA->nr; i++) {
913: for (j=0; j<bA->nc; j++) {
914: MatType type;
915: char name[256] = "",prefix[256] = "";
916: PetscInt NR,NC;
917: PetscBool isNest = PETSC_FALSE;
919: if (!bA->m[i][j]) {
920: PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL \n",i,j);
921: continue;
922: }
923: MatGetSize(bA->m[i][j],&NR,&NC);
924: MatGetType(bA->m[i][j], &type);
925: if (((PetscObject)bA->m[i][j])->name) PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);
926: if (((PetscObject)bA->m[i][j])->prefix) PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);
927: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
929: PetscViewerASCIIPrintf(viewer,"(%" PetscInt_FMT ",%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n",i,j,name,prefix,type,NR,NC);
931: if (isNest || viewSub) {
932: PetscViewerASCIIPushTab(viewer); /* push1 */
933: MatView(bA->m[i][j],viewer);
934: PetscViewerASCIIPopTab(viewer); /* pop1 */
935: }
936: }
937: }
938: PetscViewerASCIIPopTab(viewer); /* pop0 */
939: }
940: return 0;
941: }
943: static PetscErrorCode MatZeroEntries_Nest(Mat A)
944: {
945: Mat_Nest *bA = (Mat_Nest*)A->data;
946: PetscInt i,j;
948: for (i=0; i<bA->nr; i++) {
949: for (j=0; j<bA->nc; j++) {
950: if (!bA->m[i][j]) continue;
951: MatZeroEntries(bA->m[i][j]);
952: }
953: }
954: return 0;
955: }
957: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
958: {
959: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
960: PetscInt i,j,nr = bA->nr,nc = bA->nc;
961: PetscBool nnzstate = PETSC_FALSE;
964: for (i=0; i<nr; i++) {
965: for (j=0; j<nc; j++) {
966: PetscObjectState subnnzstate = 0;
967: if (bA->m[i][j] && bB->m[i][j]) {
968: MatCopy(bA->m[i][j],bB->m[i][j],str);
970: MatGetNonzeroState(bB->m[i][j],&subnnzstate);
971: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
972: bB->nnzstate[i*nc+j] = subnnzstate;
973: }
974: }
975: if (nnzstate) B->nonzerostate++;
976: return 0;
977: }
979: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
980: {
981: Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
982: PetscInt i,j,nr = bY->nr,nc = bY->nc;
983: PetscBool nnzstate = PETSC_FALSE;
986: for (i=0; i<nr; i++) {
987: for (j=0; j<nc; j++) {
988: PetscObjectState subnnzstate = 0;
989: if (bY->m[i][j] && bX->m[i][j]) {
990: MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
991: } else if (bX->m[i][j]) {
992: Mat M;
995: MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
996: MatNestSetSubMat(Y,i,j,M);
997: MatDestroy(&M);
998: }
999: if (bY->m[i][j]) MatGetNonzeroState(bY->m[i][j],&subnnzstate);
1000: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1001: bY->nnzstate[i*nc+j] = subnnzstate;
1002: }
1003: }
1004: if (nnzstate) Y->nonzerostate++;
1005: return 0;
1006: }
1008: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1009: {
1010: Mat_Nest *bA = (Mat_Nest*)A->data;
1011: Mat *b;
1012: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1014: PetscMalloc1(nr*nc,&b);
1015: for (i=0; i<nr; i++) {
1016: for (j=0; j<nc; j++) {
1017: if (bA->m[i][j]) {
1018: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1019: } else {
1020: b[i*nc+j] = NULL;
1021: }
1022: }
1023: }
1024: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1025: /* Give the new MatNest exclusive ownership */
1026: for (i=0; i<nr*nc; i++) {
1027: MatDestroy(&b[i]);
1028: }
1029: PetscFree(b);
1031: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1032: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1033: return 0;
1034: }
1036: /* nest api */
1037: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1038: {
1039: Mat_Nest *bA = (Mat_Nest*)A->data;
1043: *mat = bA->m[idxm][jdxm];
1044: return 0;
1045: }
1047: /*@
1048: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
1050: Not collective
1052: Input Parameters:
1053: + A - nest matrix
1054: . idxm - index of the matrix within the nest matrix
1055: - jdxm - index of the matrix within the nest matrix
1057: Output Parameter:
1058: . sub - matrix at index idxm,jdxm within the nest matrix
1060: Level: developer
1062: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1063: MatNestGetLocalISs(), MatNestGetISs()
1064: @*/
1065: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1066: {
1067: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1068: return 0;
1069: }
1071: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1072: {
1073: Mat_Nest *bA = (Mat_Nest*)A->data;
1074: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1078: MatGetLocalSize(mat,&m,&n);
1079: MatGetSize(mat,&M,&N);
1080: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1081: ISGetSize(bA->isglobal.row[idxm],&Mi);
1082: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1083: ISGetSize(bA->isglobal.col[jdxm],&Ni);
1087: /* do not increase object state */
1088: if (mat == bA->m[idxm][jdxm]) return 0;
1090: PetscObjectReference((PetscObject)mat);
1091: MatDestroy(&bA->m[idxm][jdxm]);
1092: bA->m[idxm][jdxm] = mat;
1093: PetscObjectStateIncrease((PetscObject)A);
1094: MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1095: A->nonzerostate++;
1096: return 0;
1097: }
1099: /*@
1100: MatNestSetSubMat - Set a single submatrix in the nest matrix.
1102: Logically collective on the submatrix communicator
1104: Input Parameters:
1105: + A - nest matrix
1106: . idxm - index of the matrix within the nest matrix
1107: . jdxm - index of the matrix within the nest matrix
1108: - sub - matrix at index idxm,jdxm within the nest matrix
1110: Notes:
1111: The new submatrix must have the same size and communicator as that block of the nest.
1113: This increments the reference count of the submatrix.
1115: Level: developer
1117: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1118: MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1119: @*/
1120: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1121: {
1122: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1123: return 0;
1124: }
1126: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1127: {
1128: Mat_Nest *bA = (Mat_Nest*)A->data;
1130: if (M) *M = bA->nr;
1131: if (N) *N = bA->nc;
1132: if (mat) *mat = bA->m;
1133: return 0;
1134: }
1136: /*@C
1137: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
1139: Not collective
1141: Input Parameter:
1142: . A - nest matrix
1144: Output Parameters:
1145: + M - number of rows in the nest matrix
1146: . N - number of cols in the nest matrix
1147: - mat - 2d array of matrices
1149: Notes:
1151: The user should not free the array mat.
1153: In Fortran, this routine has a calling sequence
1154: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1155: where the space allocated for the optional argument mat is assumed large enough (if provided).
1157: Level: developer
1159: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1160: MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1161: @*/
1162: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1163: {
1164: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1165: return 0;
1166: }
1168: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1169: {
1170: Mat_Nest *bA = (Mat_Nest*)A->data;
1172: if (M) *M = bA->nr;
1173: if (N) *N = bA->nc;
1174: return 0;
1175: }
1177: /*@
1178: MatNestGetSize - Returns the size of the nest matrix.
1180: Not collective
1182: Input Parameter:
1183: . A - nest matrix
1185: Output Parameters:
1186: + M - number of rows in the nested mat
1187: - N - number of cols in the nested mat
1189: Notes:
1191: Level: developer
1193: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1194: MatNestGetISs()
1195: @*/
1196: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1197: {
1198: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1199: return 0;
1200: }
1202: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1203: {
1204: Mat_Nest *vs = (Mat_Nest*)A->data;
1205: PetscInt i;
1207: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1208: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1209: return 0;
1210: }
1212: /*@C
1213: MatNestGetISs - Returns the index sets partitioning the row and column spaces
1215: Not collective
1217: Input Parameter:
1218: . A - nest matrix
1220: Output Parameters:
1221: + rows - array of row index sets
1222: - cols - array of column index sets
1224: Level: advanced
1226: Notes:
1227: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1229: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1230: MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1231: @*/
1232: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1233: {
1235: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1236: return 0;
1237: }
1239: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1240: {
1241: Mat_Nest *vs = (Mat_Nest*)A->data;
1242: PetscInt i;
1244: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1245: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1246: return 0;
1247: }
1249: /*@C
1250: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1252: Not collective
1254: Input Parameter:
1255: . A - nest matrix
1257: Output Parameters:
1258: + rows - array of row index sets (or NULL to ignore)
1259: - cols - array of column index sets (or NULL to ignore)
1261: Level: advanced
1263: Notes:
1264: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1266: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1267: MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1268: @*/
1269: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1270: {
1272: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1273: return 0;
1274: }
1276: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1277: {
1278: PetscBool flg;
1280: PetscStrcmp(vtype,VECNEST,&flg);
1281: /* In reality, this only distinguishes VECNEST and "other" */
1282: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1283: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1284: return 0;
1285: }
1287: /*@C
1288: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1290: Not collective
1292: Input Parameters:
1293: + A - nest matrix
1294: - vtype - type to use for creating vectors
1296: Notes:
1298: Level: developer
1300: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1301: @*/
1302: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1303: {
1304: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1305: return 0;
1306: }
1308: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1309: {
1310: Mat_Nest *s = (Mat_Nest*)A->data;
1311: PetscInt i,j,m,n,M,N;
1312: PetscBool cong,isstd,sametype=PETSC_FALSE;
1313: VecType vtype,type;
1315: MatReset_Nest(A);
1317: s->nr = nr;
1318: s->nc = nc;
1320: /* Create space for submatrices */
1321: PetscMalloc1(nr,&s->m);
1322: for (i=0; i<nr; i++) {
1323: PetscMalloc1(nc,&s->m[i]);
1324: }
1325: for (i=0; i<nr; i++) {
1326: for (j=0; j<nc; j++) {
1327: s->m[i][j] = a[i*nc+j];
1328: if (a[i*nc+j]) {
1329: PetscObjectReference((PetscObject)a[i*nc+j]);
1330: }
1331: }
1332: }
1333: MatGetVecType(A,&vtype);
1334: PetscStrcmp(vtype,VECSTANDARD,&isstd);
1335: if (isstd) {
1336: /* check if all blocks have the same vectype */
1337: vtype = NULL;
1338: for (i=0; i<nr; i++) {
1339: for (j=0; j<nc; j++) {
1340: if (a[i*nc+j]) {
1341: if (!vtype) { /* first visited block */
1342: MatGetVecType(a[i*nc+j],&vtype);
1343: sametype = PETSC_TRUE;
1344: } else if (sametype) {
1345: MatGetVecType(a[i*nc+j],&type);
1346: PetscStrcmp(vtype,type,&sametype);
1347: }
1348: }
1349: }
1350: }
1351: if (sametype) { /* propagate vectype */
1352: MatSetVecType(A,vtype);
1353: }
1354: }
1356: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1358: PetscMalloc1(nr,&s->row_len);
1359: PetscMalloc1(nc,&s->col_len);
1360: for (i=0; i<nr; i++) s->row_len[i]=-1;
1361: for (j=0; j<nc; j++) s->col_len[j]=-1;
1363: PetscCalloc1(nr*nc,&s->nnzstate);
1364: for (i=0; i<nr; i++) {
1365: for (j=0; j<nc; j++) {
1366: if (s->m[i][j]) {
1367: MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1368: }
1369: }
1370: }
1372: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1374: PetscLayoutSetSize(A->rmap,M);
1375: PetscLayoutSetLocalSize(A->rmap,m);
1376: PetscLayoutSetSize(A->cmap,N);
1377: PetscLayoutSetLocalSize(A->cmap,n);
1379: PetscLayoutSetUp(A->rmap);
1380: PetscLayoutSetUp(A->cmap);
1382: /* disable operations that are not supported for non-square matrices,
1383: or matrices for which is_row != is_col */
1384: MatHasCongruentLayouts(A,&cong);
1385: if (cong && nr != nc) cong = PETSC_FALSE;
1386: if (cong) {
1387: for (i = 0; cong && i < nr; i++) {
1388: ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1389: }
1390: }
1391: if (!cong) {
1392: A->ops->missingdiagonal = NULL;
1393: A->ops->getdiagonal = NULL;
1394: A->ops->shift = NULL;
1395: A->ops->diagonalset = NULL;
1396: }
1398: PetscCalloc2(nr,&s->left,nc,&s->right);
1399: PetscObjectStateIncrease((PetscObject)A);
1400: A->nonzerostate++;
1401: return 0;
1402: }
1404: /*@
1405: MatNestSetSubMats - Sets the nested submatrices
1407: Collective on Mat
1409: Input Parameters:
1410: + A - nested matrix
1411: . nr - number of nested row blocks
1412: . is_row - index sets for each nested row block, or NULL to make contiguous
1413: . nc - number of nested column blocks
1414: . is_col - index sets for each nested column block, or NULL to make contiguous
1415: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1417: Notes: this always resets any submatrix information previously set
1419: Level: advanced
1421: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1422: @*/
1423: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1424: {
1425: PetscInt i;
1429: if (nr && is_row) {
1432: }
1434: if (nc && is_col) {
1437: }
1439: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1440: return 0;
1441: }
1443: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1444: {
1445: PetscBool flg;
1446: PetscInt i,j,m,mi,*ix;
1448: *ltog = NULL;
1449: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1450: if (islocal[i]) {
1451: ISGetLocalSize(islocal[i],&mi);
1452: flg = PETSC_TRUE; /* We found a non-trivial entry */
1453: } else {
1454: ISGetLocalSize(isglobal[i],&mi);
1455: }
1456: m += mi;
1457: }
1458: if (!flg) return 0;
1460: PetscMalloc1(m,&ix);
1461: for (i=0,m=0; i<n; i++) {
1462: ISLocalToGlobalMapping smap = NULL;
1463: Mat sub = NULL;
1464: PetscSF sf;
1465: PetscLayout map;
1466: const PetscInt *ix2;
1468: if (!colflg) {
1469: MatNestFindNonzeroSubMatRow(A,i,&sub);
1470: } else {
1471: MatNestFindNonzeroSubMatCol(A,i,&sub);
1472: }
1473: if (sub) {
1474: if (!colflg) {
1475: MatGetLocalToGlobalMapping(sub,&smap,NULL);
1476: } else {
1477: MatGetLocalToGlobalMapping(sub,NULL,&smap);
1478: }
1479: }
1480: /*
1481: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1482: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1483: */
1484: ISGetIndices(isglobal[i],&ix2);
1485: if (islocal[i]) {
1486: PetscInt *ilocal,*iremote;
1487: PetscInt mil,nleaves;
1489: ISGetLocalSize(islocal[i],&mi);
1491: for (j=0; j<mi; j++) ix[m+j] = j;
1492: ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);
1494: /* PetscSFSetGraphLayout does not like negative indices */
1495: PetscMalloc2(mi,&ilocal,mi,&iremote);
1496: for (j=0, nleaves = 0; j<mi; j++) {
1497: if (ix[m+j] < 0) continue;
1498: ilocal[nleaves] = j;
1499: iremote[nleaves] = ix[m+j];
1500: nleaves++;
1501: }
1502: ISGetLocalSize(isglobal[i],&mil);
1503: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1504: PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1505: PetscLayoutSetLocalSize(map,mil);
1506: PetscLayoutSetUp(map);
1507: PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1508: PetscLayoutDestroy(&map);
1509: PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1510: PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1511: PetscSFDestroy(&sf);
1512: PetscFree2(ilocal,iremote);
1513: } else {
1514: ISGetLocalSize(isglobal[i],&mi);
1515: for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1516: }
1517: ISRestoreIndices(isglobal[i],&ix2);
1518: m += mi;
1519: }
1520: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1521: return 0;
1522: }
1524: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1525: /*
1526: nprocessors = NP
1527: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1528: proc 0: => (g_0,h_0,)
1529: proc 1: => (g_1,h_1,)
1530: ...
1531: proc nprocs-1: => (g_NP-1,h_NP-1,)
1533: proc 0: proc 1: proc nprocs-1:
1534: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1536: proc 0:
1537: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1538: proc 1:
1539: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1541: proc NP-1:
1542: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1543: */
1544: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1545: {
1546: Mat_Nest *vs = (Mat_Nest*)A->data;
1547: PetscInt i,j,offset,n,nsum,bs;
1548: Mat sub = NULL;
1550: PetscMalloc1(nr,&vs->isglobal.row);
1551: PetscMalloc1(nc,&vs->isglobal.col);
1552: if (is_row) { /* valid IS is passed in */
1553: /* refs on is[] are incremented */
1554: for (i=0; i<vs->nr; i++) {
1555: PetscObjectReference((PetscObject)is_row[i]);
1557: vs->isglobal.row[i] = is_row[i];
1558: }
1559: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1560: nsum = 0;
1561: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1562: MatNestFindNonzeroSubMatRow(A,i,&sub);
1564: MatGetLocalSize(sub,&n,NULL);
1566: nsum += n;
1567: }
1568: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1569: offset -= nsum;
1570: for (i=0; i<vs->nr; i++) {
1571: MatNestFindNonzeroSubMatRow(A,i,&sub);
1572: MatGetLocalSize(sub,&n,NULL);
1573: MatGetBlockSizes(sub,&bs,NULL);
1574: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1575: ISSetBlockSize(vs->isglobal.row[i],bs);
1576: offset += n;
1577: }
1578: }
1580: if (is_col) { /* valid IS is passed in */
1581: /* refs on is[] are incremented */
1582: for (j=0; j<vs->nc; j++) {
1583: PetscObjectReference((PetscObject)is_col[j]);
1585: vs->isglobal.col[j] = is_col[j];
1586: }
1587: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1588: offset = A->cmap->rstart;
1589: nsum = 0;
1590: for (j=0; j<vs->nc; j++) {
1591: MatNestFindNonzeroSubMatCol(A,j,&sub);
1593: MatGetLocalSize(sub,NULL,&n);
1595: nsum += n;
1596: }
1597: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1598: offset -= nsum;
1599: for (j=0; j<vs->nc; j++) {
1600: MatNestFindNonzeroSubMatCol(A,j,&sub);
1601: MatGetLocalSize(sub,NULL,&n);
1602: MatGetBlockSizes(sub,NULL,&bs);
1603: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1604: ISSetBlockSize(vs->isglobal.col[j],bs);
1605: offset += n;
1606: }
1607: }
1609: /* Set up the local ISs */
1610: PetscMalloc1(vs->nr,&vs->islocal.row);
1611: PetscMalloc1(vs->nc,&vs->islocal.col);
1612: for (i=0,offset=0; i<vs->nr; i++) {
1613: IS isloc;
1614: ISLocalToGlobalMapping rmap = NULL;
1615: PetscInt nlocal,bs;
1616: MatNestFindNonzeroSubMatRow(A,i,&sub);
1617: if (sub) MatGetLocalToGlobalMapping(sub,&rmap,NULL);
1618: if (rmap) {
1619: MatGetBlockSizes(sub,&bs,NULL);
1620: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1621: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1622: ISSetBlockSize(isloc,bs);
1623: } else {
1624: nlocal = 0;
1625: isloc = NULL;
1626: }
1627: vs->islocal.row[i] = isloc;
1628: offset += nlocal;
1629: }
1630: for (i=0,offset=0; i<vs->nc; i++) {
1631: IS isloc;
1632: ISLocalToGlobalMapping cmap = NULL;
1633: PetscInt nlocal,bs;
1634: MatNestFindNonzeroSubMatCol(A,i,&sub);
1635: if (sub) MatGetLocalToGlobalMapping(sub,NULL,&cmap);
1636: if (cmap) {
1637: MatGetBlockSizes(sub,NULL,&bs);
1638: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1639: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1640: ISSetBlockSize(isloc,bs);
1641: } else {
1642: nlocal = 0;
1643: isloc = NULL;
1644: }
1645: vs->islocal.col[i] = isloc;
1646: offset += nlocal;
1647: }
1649: /* Set up the aggregate ISLocalToGlobalMapping */
1650: {
1651: ISLocalToGlobalMapping rmap,cmap;
1652: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1653: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1654: if (rmap && cmap) MatSetLocalToGlobalMapping(A,rmap,cmap);
1655: ISLocalToGlobalMappingDestroy(&rmap);
1656: ISLocalToGlobalMappingDestroy(&cmap);
1657: }
1659: if (PetscDefined(USE_DEBUG)) {
1660: for (i=0; i<vs->nr; i++) {
1661: for (j=0; j<vs->nc; j++) {
1662: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1663: Mat B = vs->m[i][j];
1664: if (!B) continue;
1665: MatGetSize(B,&M,&N);
1666: MatGetLocalSize(B,&m,&n);
1667: ISGetSize(vs->isglobal.row[i],&Mi);
1668: ISGetSize(vs->isglobal.col[j],&Ni);
1669: ISGetLocalSize(vs->isglobal.row[i],&mi);
1670: ISGetLocalSize(vs->isglobal.col[j],&ni);
1673: }
1674: }
1675: }
1677: /* Set A->assembled if all non-null blocks are currently assembled */
1678: for (i=0; i<vs->nr; i++) {
1679: for (j=0; j<vs->nc; j++) {
1680: if (vs->m[i][j] && !vs->m[i][j]->assembled) return 0;
1681: }
1682: }
1683: A->assembled = PETSC_TRUE;
1684: return 0;
1685: }
1687: /*@C
1688: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1690: Collective on Mat
1692: Input Parameters:
1693: + comm - Communicator for the new Mat
1694: . nr - number of nested row blocks
1695: . is_row - index sets for each nested row block, or NULL to make contiguous
1696: . nc - number of nested column blocks
1697: . is_col - index sets for each nested column block, or NULL to make contiguous
1698: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1700: Output Parameter:
1701: . B - new matrix
1703: Level: advanced
1705: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1706: MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1707: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1708: @*/
1709: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1710: {
1711: Mat A;
1713: *B = NULL;
1714: MatCreate(comm,&A);
1715: MatSetType(A,MATNEST);
1716: A->preallocated = PETSC_TRUE;
1717: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1718: *B = A;
1719: return 0;
1720: }
1722: PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1723: {
1724: Mat_Nest *nest = (Mat_Nest*)A->data;
1725: Mat *trans;
1726: PetscScalar **avv;
1727: PetscScalar *vv;
1728: PetscInt **aii,**ajj;
1729: PetscInt *ii,*jj,*ci;
1730: PetscInt nr,nc,nnz,i,j;
1731: PetscBool done;
1733: MatGetSize(A,&nr,&nc);
1734: if (reuse == MAT_REUSE_MATRIX) {
1735: PetscInt rnr;
1737: MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1740: MatSeqAIJGetArray(*newmat,&vv);
1741: }
1742: /* extract CSR for nested SeqAIJ matrices */
1743: nnz = 0;
1744: PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1745: for (i=0; i<nest->nr; ++i) {
1746: for (j=0; j<nest->nc; ++j) {
1747: Mat B = nest->m[i][j];
1748: if (B) {
1749: PetscScalar *naa;
1750: PetscInt *nii,*njj,nnr;
1751: PetscBool istrans;
1753: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1754: if (istrans) {
1755: Mat Bt;
1757: MatTransposeGetMat(B,&Bt);
1758: MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1759: B = trans[i*nest->nc+j];
1760: }
1761: MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1763: MatSeqAIJGetArray(B,&naa);
1764: nnz += nii[nnr];
1766: aii[i*nest->nc+j] = nii;
1767: ajj[i*nest->nc+j] = njj;
1768: avv[i*nest->nc+j] = naa;
1769: }
1770: }
1771: }
1772: if (reuse != MAT_REUSE_MATRIX) {
1773: PetscMalloc1(nr+1,&ii);
1774: PetscMalloc1(nnz,&jj);
1775: PetscMalloc1(nnz,&vv);
1776: } else {
1778: }
1780: /* new row pointer */
1781: PetscArrayzero(ii,nr+1);
1782: for (i=0; i<nest->nr; ++i) {
1783: PetscInt ncr,rst;
1785: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1786: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1787: for (j=0; j<nest->nc; ++j) {
1788: if (aii[i*nest->nc+j]) {
1789: PetscInt *nii = aii[i*nest->nc+j];
1790: PetscInt ir;
1792: for (ir=rst; ir<ncr+rst; ++ir) {
1793: ii[ir+1] += nii[1]-nii[0];
1794: nii++;
1795: }
1796: }
1797: }
1798: }
1799: for (i=0; i<nr; i++) ii[i+1] += ii[i];
1801: /* construct CSR for the new matrix */
1802: PetscCalloc1(nr,&ci);
1803: for (i=0; i<nest->nr; ++i) {
1804: PetscInt ncr,rst;
1806: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1807: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1808: for (j=0; j<nest->nc; ++j) {
1809: if (aii[i*nest->nc+j]) {
1810: PetscScalar *nvv = avv[i*nest->nc+j];
1811: PetscInt *nii = aii[i*nest->nc+j];
1812: PetscInt *njj = ajj[i*nest->nc+j];
1813: PetscInt ir,cst;
1815: ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1816: for (ir=rst; ir<ncr+rst; ++ir) {
1817: PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1819: for (ij=0;ij<rsize;ij++) {
1820: jj[ist+ij] = *njj+cst;
1821: vv[ist+ij] = *nvv;
1822: njj++;
1823: nvv++;
1824: }
1825: ci[ir] += rsize;
1826: nii++;
1827: }
1828: }
1829: }
1830: }
1831: PetscFree(ci);
1833: /* restore info */
1834: for (i=0; i<nest->nr; ++i) {
1835: for (j=0; j<nest->nc; ++j) {
1836: Mat B = nest->m[i][j];
1837: if (B) {
1838: PetscInt nnr = 0, k = i*nest->nc+j;
1840: B = (trans[k] ? trans[k] : B);
1841: MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1843: MatSeqAIJRestoreArray(B,&avv[k]);
1844: MatDestroy(&trans[k]);
1845: }
1846: }
1847: }
1848: PetscFree4(aii,ajj,avv,trans);
1850: /* finalize newmat */
1851: if (reuse == MAT_INITIAL_MATRIX) {
1852: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1853: } else if (reuse == MAT_INPLACE_MATRIX) {
1854: Mat B;
1856: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1857: MatHeaderReplace(A,&B);
1858: }
1859: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1860: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1861: {
1862: Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1863: a->free_a = PETSC_TRUE;
1864: a->free_ij = PETSC_TRUE;
1865: }
1866: return 0;
1867: }
1869: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y,PetscScalar a,Mat X)
1870: {
1871: Mat_Nest *nest = (Mat_Nest*)X->data;
1872: PetscInt i,j,k,rstart;
1873: PetscBool flg;
1875: /* Fill by row */
1876: for (j=0; j<nest->nc; ++j) {
1877: /* Using global column indices and ISAllGather() is not scalable. */
1878: IS bNis;
1879: PetscInt bN;
1880: const PetscInt *bNindices;
1881: ISAllGather(nest->isglobal.col[j], &bNis);
1882: ISGetSize(bNis,&bN);
1883: ISGetIndices(bNis,&bNindices);
1884: for (i=0; i<nest->nr; ++i) {
1885: Mat B,D=NULL;
1886: PetscInt bm, br;
1887: const PetscInt *bmindices;
1888: B = nest->m[i][j];
1889: if (!B) continue;
1890: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&flg);
1891: if (flg) {
1892: PetscTryMethod(B,"MatTransposeGetMat_C",(Mat,Mat*),(B,&D));
1893: PetscTryMethod(B,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(B,&D));
1894: MatConvert(B,((PetscObject)D)->type_name,MAT_INITIAL_MATRIX,&D);
1895: B = D;
1896: }
1897: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQSBAIJ,MATMPISBAIJ,"");
1898: if (flg) {
1899: if (D) {
1900: MatConvert(D,MATBAIJ,MAT_INPLACE_MATRIX,&D);
1901: } else {
1902: MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&D);
1903: }
1904: B = D;
1905: }
1906: ISGetLocalSize(nest->isglobal.row[i],&bm);
1907: ISGetIndices(nest->isglobal.row[i],&bmindices);
1908: MatGetOwnershipRange(B,&rstart,NULL);
1909: for (br = 0; br < bm; ++br) {
1910: PetscInt row = bmindices[br], brncols, *cols;
1911: const PetscInt *brcols;
1912: const PetscScalar *brcoldata;
1913: PetscScalar *vals = NULL;
1914: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1915: PetscMalloc1(brncols,&cols);
1916: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1917: /*
1918: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1919: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1920: */
1921: if (a != 1.0) {
1922: PetscMalloc1(brncols,&vals);
1923: for (k=0; k<brncols; k++) vals[k] = a * brcoldata[k];
1924: MatSetValues(Y,1,&row,brncols,cols,vals,ADD_VALUES);
1925: PetscFree(vals);
1926: } else {
1927: MatSetValues(Y,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1928: }
1929: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1930: PetscFree(cols);
1931: }
1932: if (D) {
1933: MatDestroy(&D);
1934: }
1935: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1936: }
1937: ISRestoreIndices(bNis,&bNindices);
1938: ISDestroy(&bNis);
1939: }
1940: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
1941: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
1942: return 0;
1943: }
1945: PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1946: {
1947: Mat_Nest *nest = (Mat_Nest*)A->data;
1948: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart,cstart,cend;
1949: PetscMPIInt size;
1950: Mat C;
1952: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1953: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1954: PetscInt nf;
1955: PetscBool fast;
1957: PetscStrcmp(newtype,MATAIJ,&fast);
1958: if (!fast) {
1959: PetscStrcmp(newtype,MATSEQAIJ,&fast);
1960: }
1961: for (i=0; i<nest->nr && fast; ++i) {
1962: for (j=0; j<nest->nc && fast; ++j) {
1963: Mat B = nest->m[i][j];
1964: if (B) {
1965: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1966: if (!fast) {
1967: PetscBool istrans;
1969: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1970: if (istrans) {
1971: Mat Bt;
1973: MatTransposeGetMat(B,&Bt);
1974: PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1975: }
1976: }
1977: }
1978: }
1979: }
1980: for (i=0, nf=0; i<nest->nr && fast; ++i) {
1981: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1982: if (fast) {
1983: PetscInt f,s;
1985: ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1986: if (f != nf || s != 1) { fast = PETSC_FALSE; }
1987: else {
1988: ISGetSize(nest->isglobal.row[i],&f);
1989: nf += f;
1990: }
1991: }
1992: }
1993: for (i=0, nf=0; i<nest->nc && fast; ++i) {
1994: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1995: if (fast) {
1996: PetscInt f,s;
1998: ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1999: if (f != nf || s != 1) { fast = PETSC_FALSE; }
2000: else {
2001: ISGetSize(nest->isglobal.col[i],&f);
2002: nf += f;
2003: }
2004: }
2005: }
2006: if (fast) {
2007: MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
2008: return 0;
2009: }
2010: }
2011: MatGetSize(A,&M,&N);
2012: MatGetLocalSize(A,&m,&n);
2013: MatGetOwnershipRangeColumn(A,&cstart,&cend);
2014: if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2015: else {
2016: MatCreate(PetscObjectComm((PetscObject)A),&C);
2017: MatSetType(C,newtype);
2018: MatSetSizes(C,m,n,M,N);
2019: }
2020: PetscMalloc1(2*m,&dnnz);
2021: onnz = dnnz + m;
2022: for (k=0; k<m; k++) {
2023: dnnz[k] = 0;
2024: onnz[k] = 0;
2025: }
2026: for (j=0; j<nest->nc; ++j) {
2027: IS bNis;
2028: PetscInt bN;
2029: const PetscInt *bNindices;
2030: /* Using global column indices and ISAllGather() is not scalable. */
2031: ISAllGather(nest->isglobal.col[j], &bNis);
2032: ISGetSize(bNis, &bN);
2033: ISGetIndices(bNis,&bNindices);
2034: for (i=0; i<nest->nr; ++i) {
2035: PetscSF bmsf;
2036: PetscSFNode *iremote;
2037: Mat B;
2038: PetscInt bm, *sub_dnnz,*sub_onnz, br;
2039: const PetscInt *bmindices;
2040: B = nest->m[i][j];
2041: if (!B) continue;
2042: ISGetLocalSize(nest->isglobal.row[i],&bm);
2043: ISGetIndices(nest->isglobal.row[i],&bmindices);
2044: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2045: PetscMalloc1(bm,&iremote);
2046: PetscMalloc1(bm,&sub_dnnz);
2047: PetscMalloc1(bm,&sub_onnz);
2048: for (k = 0; k < bm; ++k) {
2049: sub_dnnz[k] = 0;
2050: sub_onnz[k] = 0;
2051: }
2052: /*
2053: Locate the owners for all of the locally-owned global row indices for this row block.
2054: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2055: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2056: */
2057: MatGetOwnershipRange(B,&rstart,NULL);
2058: for (br = 0; br < bm; ++br) {
2059: PetscInt row = bmindices[br], brncols, col;
2060: const PetscInt *brcols;
2061: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2062: PetscMPIInt rowowner = 0;
2063: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2064: /* how many roots */
2065: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2066: /* get nonzero pattern */
2067: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2068: for (k=0; k<brncols; k++) {
2069: col = bNindices[brcols[k]];
2070: if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2071: sub_dnnz[br]++;
2072: } else {
2073: sub_onnz[br]++;
2074: }
2075: }
2076: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2077: }
2078: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2079: /* bsf will have to take care of disposing of bedges. */
2080: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2081: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2082: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2083: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2084: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2085: PetscFree(sub_dnnz);
2086: PetscFree(sub_onnz);
2087: PetscSFDestroy(&bmsf);
2088: }
2089: ISRestoreIndices(bNis,&bNindices);
2090: ISDestroy(&bNis);
2091: }
2092: /* Resize preallocation if overestimated */
2093: for (i=0;i<m;i++) {
2094: dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2095: onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2096: }
2097: MatSeqAIJSetPreallocation(C,0,dnnz);
2098: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2099: PetscFree(dnnz);
2100: MatAXPY_Dense_Nest(C,1.0,A);
2101: if (reuse == MAT_INPLACE_MATRIX) {
2102: MatHeaderReplace(A,&C);
2103: } else *newmat = C;
2104: return 0;
2105: }
2107: PetscErrorCode MatConvert_Nest_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2108: {
2109: Mat B;
2110: PetscInt m,n,M,N;
2112: MatGetSize(A,&M,&N);
2113: MatGetLocalSize(A,&m,&n);
2114: if (reuse == MAT_REUSE_MATRIX) {
2115: B = *newmat;
2116: MatZeroEntries(B);
2117: } else {
2118: MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,&B);
2119: }
2120: MatAXPY_Dense_Nest(B,1.0,A);
2121: if (reuse == MAT_INPLACE_MATRIX) {
2122: MatHeaderReplace(A,&B);
2123: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2124: return 0;
2125: }
2127: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2128: {
2129: Mat_Nest *bA = (Mat_Nest*)mat->data;
2130: MatOperation opAdd;
2131: PetscInt i,j,nr = bA->nr,nc = bA->nc;
2132: PetscBool flg;
2134: *has = PETSC_FALSE;
2135: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2136: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2137: for (j=0; j<nc; j++) {
2138: for (i=0; i<nr; i++) {
2139: if (!bA->m[i][j]) continue;
2140: MatHasOperation(bA->m[i][j],opAdd,&flg);
2141: if (!flg) return 0;
2142: }
2143: }
2144: }
2145: if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2146: return 0;
2147: }
2149: /*MC
2150: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2152: Level: intermediate
2154: Notes:
2155: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2156: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2157: It is usually used with DMComposite and DMCreateMatrix()
2159: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2160: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2161: than the nest matrix.
2163: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2164: VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2165: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2166: M*/
2167: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2168: {
2169: Mat_Nest *s;
2171: PetscNewLog(A,&s);
2172: A->data = (void*)s;
2174: s->nr = -1;
2175: s->nc = -1;
2176: s->m = NULL;
2177: s->splitassembly = PETSC_FALSE;
2179: PetscMemzero(A->ops,sizeof(*A->ops));
2181: A->ops->mult = MatMult_Nest;
2182: A->ops->multadd = MatMultAdd_Nest;
2183: A->ops->multtranspose = MatMultTranspose_Nest;
2184: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2185: A->ops->transpose = MatTranspose_Nest;
2186: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2187: A->ops->assemblyend = MatAssemblyEnd_Nest;
2188: A->ops->zeroentries = MatZeroEntries_Nest;
2189: A->ops->copy = MatCopy_Nest;
2190: A->ops->axpy = MatAXPY_Nest;
2191: A->ops->duplicate = MatDuplicate_Nest;
2192: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2193: A->ops->destroy = MatDestroy_Nest;
2194: A->ops->view = MatView_Nest;
2195: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2196: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2197: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2198: A->ops->getdiagonal = MatGetDiagonal_Nest;
2199: A->ops->diagonalscale = MatDiagonalScale_Nest;
2200: A->ops->scale = MatScale_Nest;
2201: A->ops->shift = MatShift_Nest;
2202: A->ops->diagonalset = MatDiagonalSet_Nest;
2203: A->ops->setrandom = MatSetRandom_Nest;
2204: A->ops->hasoperation = MatHasOperation_Nest;
2205: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2207: A->spptr = NULL;
2208: A->assembled = PETSC_FALSE;
2210: /* expose Nest api's */
2211: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2212: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2213: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2214: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
2215: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
2216: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2217: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
2218: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2219: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2220: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2221: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2222: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);
2223: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",MatConvert_Nest_Dense);
2224: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",MatConvert_Nest_Dense);
2225: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2226: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2227: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);
2229: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2230: return 0;
2231: }