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;
19: *m = *n = *M = *N = 0;
20: for (i=0; i<bA->nr; i++) { /* rows */
21: PetscInt sm,sM;
22: ISGetLocalSize(bA->isglobal.row[i],&sm);
23: ISGetSize(bA->isglobal.row[i],&sM);
24: *m += sm;
25: *M += sM;
26: }
27: for (j=0; j<bA->nc; j++) { /* cols */
28: PetscInt sn,sN;
29: ISGetLocalSize(bA->isglobal.col[j],&sn);
30: ISGetSize(bA->isglobal.col[j],&sN);
31: *n += sn;
32: *N += sN;
33: }
34: return(0);
35: }
37: /* operations */
38: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
39: {
40: Mat_Nest *bA = (Mat_Nest*)A->data;
41: Vec *bx = bA->right,*by = bA->left;
42: PetscInt i,j,nr = bA->nr,nc = bA->nc;
46: for (i=0; i<nr; i++) {VecGetSubVector(y,bA->isglobal.row[i],&by[i]);}
47: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
48: for (i=0; i<nr; i++) {
49: VecZeroEntries(by[i]);
50: for (j=0; j<nc; j++) {
51: if (!bA->m[i][j]) continue;
52: /* y[i] <- y[i] + A[i][j] * x[j] */
53: MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);
54: }
55: }
56: for (i=0; i<nr; i++) {VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);}
57: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
58: return(0);
59: }
61: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
62: {
63: Mat_Nest *bA = (Mat_Nest*)A->data;
64: Vec *bx = bA->right,*bz = bA->left;
65: PetscInt i,j,nr = bA->nr,nc = bA->nc;
69: for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
70: for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
71: for (i=0; i<nr; i++) {
72: if (y != z) {
73: Vec by;
74: VecGetSubVector(y,bA->isglobal.row[i],&by);
75: VecCopy(by,bz[i]);
76: VecRestoreSubVector(y,bA->isglobal.row[i],&by);
77: }
78: for (j=0; j<nc; j++) {
79: if (!bA->m[i][j]) continue;
80: /* y[i] <- y[i] + A[i][j] * x[j] */
81: MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
82: }
83: }
84: for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
85: for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
86: return(0);
87: }
89: typedef struct {
90: Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
91: PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */
92: PetscInt *dm,*dn,k; /* displacements and number of submatrices */
93: } Nest_Dense;
95: PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
96: {
97: Mat_Nest *bA;
98: Nest_Dense *contents;
99: Mat viewB,viewC,productB,workC;
100: const PetscScalar *barray;
101: PetscScalar *carray;
102: PetscInt i,j,M,N,nr,nc,ldb,ldc;
103: PetscErrorCode ierr;
104: Mat A,B;
107: MatCheckProduct(C,3);
108: A = C->product->A;
109: B = C->product->B;
110: MatGetSize(B,NULL,&N);
111: if (!N) {
112: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
114: return(0);
115: }
116: contents = (Nest_Dense*)C->product->data;
117: if (!contents) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
118: bA = (Mat_Nest*)A->data;
119: nr = bA->nr;
120: nc = bA->nc;
121: MatDenseGetLDA(B,&ldb);
122: MatDenseGetLDA(C,&ldc);
123: MatZeroEntries(C);
124: MatDenseGetArrayRead(B,&barray);
125: MatDenseGetArray(C,&carray);
126: for (i=0; i<nr; i++) {
127: ISGetSize(bA->isglobal.row[i],&M);
128: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
129: MatDenseSetLDA(viewC,ldc);
130: for (j=0; j<nc; j++) {
131: if (!bA->m[i][j]) continue;
132: ISGetSize(bA->isglobal.col[j],&M);
133: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
134: MatDenseSetLDA(viewB,ldb);
136: /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
137: workC = contents->workC[i*nc + j];
138: productB = workC->product->B;
139: workC->product->B = viewB; /* use newly created dense matrix viewB */
140: MatProductNumeric(workC);
141: MatDestroy(&viewB);
142: workC->product->B = productB; /* resume original B */
144: /* C[i] <- workC + C[i] */
145: MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
146: }
147: MatDestroy(&viewC);
148: }
149: MatDenseRestoreArray(C,&carray);
150: MatDenseRestoreArrayRead(B,&barray);
152: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
154: return(0);
155: }
157: PetscErrorCode MatNest_DenseDestroy(void *ctx)
158: {
159: Nest_Dense *contents = (Nest_Dense*)ctx;
160: PetscInt i;
164: PetscFree(contents->tarray);
165: for (i=0; i<contents->k; i++) {
166: MatDestroy(contents->workC + i);
167: }
168: PetscFree3(contents->dm,contents->dn,contents->workC);
169: PetscFree(contents);
170: return(0);
171: }
173: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
174: {
175: Mat_Nest *bA;
176: Mat viewB,workC;
177: const PetscScalar *barray;
178: PetscInt i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
179: Nest_Dense *contents=NULL;
180: PetscBool cisdense;
181: PetscErrorCode ierr;
182: Mat A,B;
183: PetscReal fill;
186: MatCheckProduct(C,4);
187: if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
188: A = C->product->A;
189: B = C->product->B;
190: fill = C->product->fill;
191: bA = (Mat_Nest*)A->data;
192: nr = bA->nr;
193: nc = bA->nc;
194: MatGetLocalSize(C,&m,&n);
195: MatGetSize(C,&M,&N);
196: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
197: MatGetLocalSize(B,NULL,&n);
198: MatGetSize(B,NULL,&N);
199: MatGetLocalSize(A,&m,NULL);
200: MatGetSize(A,&M,NULL);
201: MatSetSizes(C,m,n,M,N);
202: }
203: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"");
204: if (!cisdense) {
205: MatSetType(C,((PetscObject)B)->type_name);
206: }
207: MatSetUp(C);
208: if (!N) {
209: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
210: return(0);
211: }
213: PetscNew(&contents);
214: C->product->data = contents;
215: C->product->destroy = MatNest_DenseDestroy;
216: PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
217: contents->k = nr*nc;
218: for (i=0; i<nr; i++) {
219: ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
220: maxm = PetscMax(maxm,contents->dm[i+1]);
221: contents->dm[i+1] += contents->dm[i];
222: }
223: for (i=0; i<nc; i++) {
224: ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
225: contents->dn[i+1] += contents->dn[i];
226: }
227: PetscMalloc1(maxm*N,&contents->tarray);
228: MatDenseGetLDA(B,&ldb);
229: MatGetSize(B,NULL,&N);
230: MatDenseGetArrayRead(B,&barray);
231: /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
232: for (j=0; j<nc; j++) {
233: ISGetSize(bA->isglobal.col[j],&M);
234: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
235: MatDenseSetLDA(viewB,ldb);
236: for (i=0; i<nr; i++) {
237: if (!bA->m[i][j]) continue;
238: /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
240: MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
241: workC = contents->workC[i*nc + j];
242: MatProductSetType(workC,MATPRODUCT_AB);
243: MatProductSetAlgorithm(workC,"default");
244: MatProductSetFill(workC,fill);
245: MatProductSetFromOptions(workC);
246: MatProductSymbolic(workC);
248: /* since tarray will be shared by all Mat */
249: MatSeqDenseSetPreallocation(workC,contents->tarray);
250: MatMPIDenseSetPreallocation(workC,contents->tarray);
251: }
252: MatDestroy(&viewB);
253: }
254: MatDenseRestoreArrayRead(B,&barray);
256: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
257: return(0);
258: }
260: /* --------------------------------------------------------- */
261: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
262: {
264: C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
265: return(0);
266: }
268: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
269: {
271: Mat_Product *product = C->product;
274: if (product->type == MATPRODUCT_AB) {
275: MatProductSetFromOptions_Nest_Dense_AB(C);
276: }
277: return(0);
278: }
279: /* --------------------------------------------------------- */
281: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
282: {
283: Mat_Nest *bA = (Mat_Nest*)A->data;
284: Vec *bx = bA->left,*by = bA->right;
285: PetscInt i,j,nr = bA->nr,nc = bA->nc;
289: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
290: for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
291: for (j=0; j<nc; j++) {
292: VecZeroEntries(by[j]);
293: for (i=0; i<nr; i++) {
294: if (!bA->m[i][j]) continue;
295: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
296: MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
297: }
298: }
299: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
300: for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
301: return(0);
302: }
304: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
305: {
306: Mat_Nest *bA = (Mat_Nest*)A->data;
307: Vec *bx = bA->left,*bz = bA->right;
308: PetscInt i,j,nr = bA->nr,nc = bA->nc;
312: for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
313: for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
314: for (j=0; j<nc; j++) {
315: if (y != z) {
316: Vec by;
317: VecGetSubVector(y,bA->isglobal.col[j],&by);
318: VecCopy(by,bz[j]);
319: VecRestoreSubVector(y,bA->isglobal.col[j],&by);
320: }
321: for (i=0; i<nr; i++) {
322: if (!bA->m[i][j]) continue;
323: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
324: MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
325: }
326: }
327: for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
328: for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
329: return(0);
330: }
332: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
333: {
334: Mat_Nest *bA = (Mat_Nest*)A->data, *bC;
335: Mat C;
336: PetscInt i,j,nr = bA->nr,nc = bA->nc;
340: if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");
342: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
343: Mat *subs;
344: IS *is_row,*is_col;
346: PetscCalloc1(nr * nc,&subs);
347: PetscMalloc2(nr,&is_row,nc,&is_col);
348: MatNestGetISs(A,is_row,is_col);
349: if (reuse == MAT_INPLACE_MATRIX) {
350: for (i=0; i<nr; i++) {
351: for (j=0; j<nc; j++) {
352: subs[i + nr * j] = bA->m[i][j];
353: }
354: }
355: }
357: MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
358: PetscFree(subs);
359: PetscFree2(is_row,is_col);
360: } else {
361: C = *B;
362: }
364: bC = (Mat_Nest*)C->data;
365: for (i=0; i<nr; i++) {
366: for (j=0; j<nc; j++) {
367: if (bA->m[i][j]) {
368: MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
369: } else {
370: bC->m[j][i] = NULL;
371: }
372: }
373: }
375: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
376: *B = C;
377: } else {
378: MatHeaderMerge(A, &C);
379: }
380: return(0);
381: }
383: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
384: {
386: IS *lst = *list;
387: PetscInt i;
390: if (!lst) return(0);
391: for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
392: PetscFree(lst);
393: *list = NULL;
394: return(0);
395: }
397: static PetscErrorCode MatReset_Nest(Mat A)
398: {
399: Mat_Nest *vs = (Mat_Nest*)A->data;
400: PetscInt i,j;
404: /* release the matrices and the place holders */
405: MatNestDestroyISList(vs->nr,&vs->isglobal.row);
406: MatNestDestroyISList(vs->nc,&vs->isglobal.col);
407: MatNestDestroyISList(vs->nr,&vs->islocal.row);
408: MatNestDestroyISList(vs->nc,&vs->islocal.col);
410: PetscFree(vs->row_len);
411: PetscFree(vs->col_len);
412: PetscFree(vs->nnzstate);
414: PetscFree2(vs->left,vs->right);
416: /* release the matrices and the place holders */
417: if (vs->m) {
418: for (i=0; i<vs->nr; i++) {
419: for (j=0; j<vs->nc; j++) {
420: MatDestroy(&vs->m[i][j]);
421: }
422: PetscFree(vs->m[i]);
423: }
424: PetscFree(vs->m);
425: }
427: /* restore defaults */
428: vs->nr = 0;
429: vs->nc = 0;
430: vs->splitassembly = PETSC_FALSE;
431: return(0);
432: }
434: static PetscErrorCode MatDestroy_Nest(Mat A)
435: {
439: MatReset_Nest(A);
440: PetscFree(A->data);
441: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",NULL);
442: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",NULL);
443: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",NULL);
444: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",NULL);
445: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",NULL);
446: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",NULL);
447: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",NULL);
448: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",NULL);
449: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",NULL);
450: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",NULL);
451: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",NULL);
452: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",NULL);
453: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",NULL);
454: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",NULL);
455: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
456: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
457: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
458: return(0);
459: }
461: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
462: {
463: Mat_Nest *vs = (Mat_Nest*)mat->data;
464: PetscInt i;
468: if (dd) *dd = 0;
469: if (!vs->nr) {
470: *missing = PETSC_TRUE;
471: return(0);
472: }
473: *missing = PETSC_FALSE;
474: for (i = 0; i < vs->nr && !(*missing); i++) {
475: *missing = PETSC_TRUE;
476: if (vs->m[i][i]) {
477: MatMissingDiagonal(vs->m[i][i],missing,NULL);
478: if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
479: }
480: }
481: return(0);
482: }
484: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
485: {
486: Mat_Nest *vs = (Mat_Nest*)A->data;
487: PetscInt i,j;
489: PetscBool nnzstate = PETSC_FALSE;
492: for (i=0; i<vs->nr; i++) {
493: for (j=0; j<vs->nc; j++) {
494: PetscObjectState subnnzstate = 0;
495: if (vs->m[i][j]) {
496: MatAssemblyBegin(vs->m[i][j],type);
497: if (!vs->splitassembly) {
498: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
499: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
500: * already performing an assembly, but the result would by more complicated and appears to offer less
501: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
502: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
503: */
504: MatAssemblyEnd(vs->m[i][j],type);
505: MatGetNonzeroState(vs->m[i][j],&subnnzstate);
506: }
507: }
508: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
509: vs->nnzstate[i*vs->nc+j] = subnnzstate;
510: }
511: }
512: if (nnzstate) A->nonzerostate++;
513: return(0);
514: }
516: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
517: {
518: Mat_Nest *vs = (Mat_Nest*)A->data;
519: PetscInt i,j;
523: for (i=0; i<vs->nr; i++) {
524: for (j=0; j<vs->nc; j++) {
525: if (vs->m[i][j]) {
526: if (vs->splitassembly) {
527: MatAssemblyEnd(vs->m[i][j],type);
528: }
529: }
530: }
531: }
532: return(0);
533: }
535: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
536: {
538: Mat_Nest *vs = (Mat_Nest*)A->data;
539: PetscInt j;
540: Mat sub;
543: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
544: for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
545: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
546: *B = sub;
547: return(0);
548: }
550: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
551: {
553: Mat_Nest *vs = (Mat_Nest*)A->data;
554: PetscInt i;
555: Mat sub;
558: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
559: for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
560: if (sub) {MatSetUp(sub);} /* Ensure that the sizes are available */
561: *B = sub;
562: return(0);
563: }
565: static PetscErrorCode MatNestFindISRange(Mat A,PetscInt n,const IS list[],IS is,PetscInt *begin,PetscInt *end)
566: {
567: PetscInt i,j,size,m;
568: PetscBool flg;
569: IS out,concatenate[2];
575: if (begin) {
577: *begin = -1;
578: }
579: if (end) {
581: *end = -1;
582: }
583: for (i=0; i<n; i++) {
584: if (!list[i]) continue;
585: ISEqualUnsorted(list[i],is,&flg);
586: if (flg) {
587: if (begin) *begin = i;
588: if (end) *end = i+1;
589: return(0);
590: }
591: }
592: ISGetSize(is,&size);
593: for (i=0; i<n-1; i++) {
594: if (!list[i]) continue;
595: m = 0;
596: ISConcatenate(PetscObjectComm((PetscObject)A),2,list+i,&out);
597: ISGetSize(out,&m);
598: for (j=i+2; j<n && m<size; j++) {
599: if (list[j]) {
600: concatenate[0] = out;
601: concatenate[1] = list[j];
602: ISConcatenate(PetscObjectComm((PetscObject)A),2,concatenate,&out);
603: ISDestroy(concatenate);
604: ISGetSize(out,&m);
605: }
606: }
607: if (m == size) {
608: ISEqualUnsorted(out,is,&flg);
609: if (flg) {
610: if (begin) *begin = i;
611: if (end) *end = j;
612: ISDestroy(&out);
613: return(0);
614: }
615: }
616: ISDestroy(&out);
617: }
618: return(0);
619: }
621: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A,PetscInt i,PetscInt j,Mat *B)
622: {
623: Mat_Nest *vs = (Mat_Nest*)A->data;
624: PetscInt lr,lc;
628: MatCreate(PetscObjectComm((PetscObject)A),B);
629: ISGetLocalSize(vs->isglobal.row[i],&lr);
630: ISGetLocalSize(vs->isglobal.col[j],&lc);
631: MatSetSizes(*B,lr,lc,PETSC_DECIDE,PETSC_DECIDE);
632: MatSetType(*B,MATAIJ);
633: MatSeqAIJSetPreallocation(*B,0,NULL);
634: MatMPIAIJSetPreallocation(*B,0,NULL,0,NULL);
635: MatSetUp(*B);
636: MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
637: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
638: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
639: return(0);
640: }
642: static PetscErrorCode MatNestGetBlock_Private(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *B)
643: {
644: Mat_Nest *vs = (Mat_Nest*)A->data;
645: Mat *a;
646: PetscInt i,j,k,l,nr=rend-rbegin,nc=cend-cbegin;
647: char keyname[256];
648: PetscBool *b;
649: PetscBool flg;
653: *B = NULL;
654: PetscSNPrintf(keyname,sizeof(keyname),"NestBlock_%D-%Dx%D-%D",rbegin,rend,cbegin,cend);
655: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
656: if (*B) return(0);
658: PetscMalloc2(nr*nc,&a,nr*nc,&b);
659: for (i=0; i<nr; i++) {
660: for (j=0; j<nc; j++) {
661: a[i*nc + j] = vs->m[rbegin+i][cbegin+j];
662: b[i*nc + j] = PETSC_FALSE;
663: }
664: }
665: if (nc!=vs->nc&&nr!=vs->nr) {
666: for (i=0; i<nr; i++) {
667: for (j=0; j<nc; j++) {
668: flg = PETSC_FALSE;
669: for (k=0; (k<nr&&!flg); k++) {
670: if (a[j + k*nc]) flg = PETSC_TRUE;
671: }
672: if (flg) {
673: flg = PETSC_FALSE;
674: for (l=0; (l<nc&&!flg); l++) {
675: if (a[i*nc + l]) flg = PETSC_TRUE;
676: }
677: }
678: if (!flg) {
679: b[i*nc + j] = PETSC_TRUE;
680: MatNestFillEmptyMat_Private(A,rbegin+i,cbegin+j,a + i*nc + j);
681: }
682: }
683: }
684: }
685: MatCreateNest(PetscObjectComm((PetscObject)A),nr,nr!=vs->nr?NULL:vs->isglobal.row,nc,nc!=vs->nc?NULL:vs->isglobal.col,a,B);
686: for (i=0; i<nr; i++) {
687: for (j=0; j<nc; j++) {
688: if (b[i*nc + j]) {
689: MatDestroy(a + i*nc + j);
690: }
691: }
692: }
693: PetscFree2(a,b);
694: (*B)->assembled = A->assembled;
695: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
696: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
697: return(0);
698: }
700: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
701: {
702: Mat_Nest *vs = (Mat_Nest*)A->data;
703: PetscInt rbegin,rend,cbegin,cend;
707: MatNestFindISRange(A,vs->nr,is->row,isrow,&rbegin,&rend);
708: MatNestFindISRange(A,vs->nc,is->col,iscol,&cbegin,&cend);
709: if (rend == rbegin + 1 && cend == cbegin + 1) {
710: if (!vs->m[rbegin][cbegin]) {
711: MatNestFillEmptyMat_Private(A,rbegin,cbegin,vs->m[rbegin] + cbegin);
712: }
713: *B = vs->m[rbegin][cbegin];
714: } else if (rbegin != -1 && cbegin != -1) {
715: MatNestGetBlock_Private(A,rbegin,rend,cbegin,cend,B);
716: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
717: return(0);
718: }
720: /*
721: TODO: This does not actually returns a submatrix we can modify
722: */
723: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
724: {
726: Mat_Nest *vs = (Mat_Nest*)A->data;
727: Mat sub;
730: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
731: switch (reuse) {
732: case MAT_INITIAL_MATRIX:
733: if (sub) { PetscObjectReference((PetscObject)sub); }
734: *B = sub;
735: break;
736: case MAT_REUSE_MATRIX:
737: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
738: break;
739: case MAT_IGNORE_MATRIX: /* Nothing to do */
740: break;
741: case MAT_INPLACE_MATRIX: /* Nothing to do */
742: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
743: }
744: return(0);
745: }
747: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
748: {
750: Mat_Nest *vs = (Mat_Nest*)A->data;
751: Mat sub;
754: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
755: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
756: if (sub) {PetscObjectReference((PetscObject)sub);}
757: *B = sub;
758: return(0);
759: }
761: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
762: {
764: Mat_Nest *vs = (Mat_Nest*)A->data;
765: Mat sub;
768: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
769: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
770: if (sub) {
771: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
772: MatDestroy(B);
773: }
774: return(0);
775: }
777: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
778: {
779: Mat_Nest *bA = (Mat_Nest*)A->data;
780: PetscInt i;
784: for (i=0; i<bA->nr; i++) {
785: Vec bv;
786: VecGetSubVector(v,bA->isglobal.row[i],&bv);
787: if (bA->m[i][i]) {
788: MatGetDiagonal(bA->m[i][i],bv);
789: } else {
790: VecSet(bv,0.0);
791: }
792: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
793: }
794: return(0);
795: }
797: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
798: {
799: Mat_Nest *bA = (Mat_Nest*)A->data;
800: Vec bl,*br;
801: PetscInt i,j;
805: PetscCalloc1(bA->nc,&br);
806: if (r) {
807: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
808: }
809: bl = NULL;
810: for (i=0; i<bA->nr; i++) {
811: if (l) {
812: VecGetSubVector(l,bA->isglobal.row[i],&bl);
813: }
814: for (j=0; j<bA->nc; j++) {
815: if (bA->m[i][j]) {
816: MatDiagonalScale(bA->m[i][j],bl,br[j]);
817: }
818: }
819: if (l) {
820: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
821: }
822: }
823: if (r) {
824: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
825: }
826: PetscFree(br);
827: return(0);
828: }
830: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
831: {
832: Mat_Nest *bA = (Mat_Nest*)A->data;
833: PetscInt i,j;
837: for (i=0; i<bA->nr; i++) {
838: for (j=0; j<bA->nc; j++) {
839: if (bA->m[i][j]) {
840: MatScale(bA->m[i][j],a);
841: }
842: }
843: }
844: return(0);
845: }
847: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
848: {
849: Mat_Nest *bA = (Mat_Nest*)A->data;
850: PetscInt i;
852: PetscBool nnzstate = PETSC_FALSE;
855: for (i=0; i<bA->nr; i++) {
856: PetscObjectState subnnzstate = 0;
857: 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);
858: MatShift(bA->m[i][i],a);
859: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
860: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
861: bA->nnzstate[i*bA->nc+i] = subnnzstate;
862: }
863: if (nnzstate) A->nonzerostate++;
864: return(0);
865: }
867: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
868: {
869: Mat_Nest *bA = (Mat_Nest*)A->data;
870: PetscInt i;
872: PetscBool nnzstate = PETSC_FALSE;
875: for (i=0; i<bA->nr; i++) {
876: PetscObjectState subnnzstate = 0;
877: Vec bv;
878: VecGetSubVector(D,bA->isglobal.row[i],&bv);
879: if (bA->m[i][i]) {
880: MatDiagonalSet(bA->m[i][i],bv,is);
881: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
882: }
883: VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
884: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
885: bA->nnzstate[i*bA->nc+i] = subnnzstate;
886: }
887: if (nnzstate) A->nonzerostate++;
888: return(0);
889: }
891: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
892: {
893: Mat_Nest *bA = (Mat_Nest*)A->data;
894: PetscInt i,j;
898: for (i=0; i<bA->nr; i++) {
899: for (j=0; j<bA->nc; j++) {
900: if (bA->m[i][j]) {
901: MatSetRandom(bA->m[i][j],rctx);
902: }
903: }
904: }
905: return(0);
906: }
908: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
909: {
910: Mat_Nest *bA = (Mat_Nest*)A->data;
911: Vec *L,*R;
912: MPI_Comm comm;
913: PetscInt i,j;
917: PetscObjectGetComm((PetscObject)A,&comm);
918: if (right) {
919: /* allocate R */
920: PetscMalloc1(bA->nc, &R);
921: /* Create the right vectors */
922: for (j=0; j<bA->nc; j++) {
923: for (i=0; i<bA->nr; i++) {
924: if (bA->m[i][j]) {
925: MatCreateVecs(bA->m[i][j],&R[j],NULL);
926: break;
927: }
928: }
929: if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
930: }
931: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
932: /* hand back control to the nest vector */
933: for (j=0; j<bA->nc; j++) {
934: VecDestroy(&R[j]);
935: }
936: PetscFree(R);
937: }
939: if (left) {
940: /* allocate L */
941: PetscMalloc1(bA->nr, &L);
942: /* Create the left vectors */
943: for (i=0; i<bA->nr; i++) {
944: for (j=0; j<bA->nc; j++) {
945: if (bA->m[i][j]) {
946: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
947: break;
948: }
949: }
950: if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
951: }
953: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
954: for (i=0; i<bA->nr; i++) {
955: VecDestroy(&L[i]);
956: }
958: PetscFree(L);
959: }
960: return(0);
961: }
963: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
964: {
965: Mat_Nest *bA = (Mat_Nest*)A->data;
966: PetscBool isascii,viewSub = PETSC_FALSE;
967: PetscInt i,j;
971: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
972: if (isascii) {
974: PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
975: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
976: PetscViewerASCIIPushTab(viewer);
977: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);
979: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
980: for (i=0; i<bA->nr; i++) {
981: for (j=0; j<bA->nc; j++) {
982: MatType type;
983: char name[256] = "",prefix[256] = "";
984: PetscInt NR,NC;
985: PetscBool isNest = PETSC_FALSE;
987: if (!bA->m[i][j]) {
988: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
989: continue;
990: }
991: MatGetSize(bA->m[i][j],&NR,&NC);
992: MatGetType(bA->m[i][j], &type);
993: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
994: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
995: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
997: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
999: if (isNest || viewSub) {
1000: PetscViewerASCIIPushTab(viewer); /* push1 */
1001: MatView(bA->m[i][j],viewer);
1002: PetscViewerASCIIPopTab(viewer); /* pop1 */
1003: }
1004: }
1005: }
1006: PetscViewerASCIIPopTab(viewer); /* pop0 */
1007: }
1008: return(0);
1009: }
1011: static PetscErrorCode MatZeroEntries_Nest(Mat A)
1012: {
1013: Mat_Nest *bA = (Mat_Nest*)A->data;
1014: PetscInt i,j;
1018: for (i=0; i<bA->nr; i++) {
1019: for (j=0; j<bA->nc; j++) {
1020: if (!bA->m[i][j]) continue;
1021: MatZeroEntries(bA->m[i][j]);
1022: }
1023: }
1024: return(0);
1025: }
1027: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
1028: {
1029: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
1030: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1032: PetscBool nnzstate = PETSC_FALSE;
1035: 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);
1036: for (i=0; i<nr; i++) {
1037: for (j=0; j<nc; j++) {
1038: PetscObjectState subnnzstate = 0;
1039: if (bA->m[i][j] && bB->m[i][j]) {
1040: MatCopy(bA->m[i][j],bB->m[i][j],str);
1041: } 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);
1042: MatGetNonzeroState(bB->m[i][j],&subnnzstate);
1043: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
1044: bB->nnzstate[i*nc+j] = subnnzstate;
1045: }
1046: }
1047: if (nnzstate) B->nonzerostate++;
1048: return(0);
1049: }
1051: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
1052: {
1053: Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
1054: PetscInt i,j,nr = bY->nr,nc = bY->nc;
1056: PetscBool nnzstate = PETSC_FALSE;
1059: 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);
1060: for (i=0; i<nr; i++) {
1061: for (j=0; j<nc; j++) {
1062: PetscObjectState subnnzstate = 0;
1063: if (bY->m[i][j] && bX->m[i][j]) {
1064: MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
1065: } else if (bX->m[i][j]) {
1066: Mat M;
1068: 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);
1069: MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
1070: MatNestSetSubMat(Y,i,j,M);
1071: MatDestroy(&M);
1072: }
1073: if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
1074: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1075: bY->nnzstate[i*nc+j] = subnnzstate;
1076: }
1077: }
1078: if (nnzstate) Y->nonzerostate++;
1079: return(0);
1080: }
1082: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1083: {
1084: Mat_Nest *bA = (Mat_Nest*)A->data;
1085: Mat *b;
1086: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1090: PetscMalloc1(nr*nc,&b);
1091: for (i=0; i<nr; i++) {
1092: for (j=0; j<nc; j++) {
1093: if (bA->m[i][j]) {
1094: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1095: } else {
1096: b[i*nc+j] = NULL;
1097: }
1098: }
1099: }
1100: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1101: /* Give the new MatNest exclusive ownership */
1102: for (i=0; i<nr*nc; i++) {
1103: MatDestroy(&b[i]);
1104: }
1105: PetscFree(b);
1107: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1108: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1109: return(0);
1110: }
1112: /* nest api */
1113: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1114: {
1115: Mat_Nest *bA = (Mat_Nest*)A->data;
1118: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1119: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1120: *mat = bA->m[idxm][jdxm];
1121: return(0);
1122: }
1124: /*@
1125: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
1127: Not collective
1129: Input Parameters:
1130: + A - nest matrix
1131: . idxm - index of the matrix within the nest matrix
1132: - jdxm - index of the matrix within the nest matrix
1134: Output Parameter:
1135: . sub - matrix at index idxm,jdxm within the nest matrix
1137: Level: developer
1139: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1140: MatNestGetLocalISs(), MatNestGetISs()
1141: @*/
1142: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1143: {
1147: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1148: return(0);
1149: }
1151: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1152: {
1153: Mat_Nest *bA = (Mat_Nest*)A->data;
1154: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1158: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1159: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1160: MatGetLocalSize(mat,&m,&n);
1161: MatGetSize(mat,&M,&N);
1162: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1163: ISGetSize(bA->isglobal.row[idxm],&Mi);
1164: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1165: ISGetSize(bA->isglobal.col[jdxm],&Ni);
1166: 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);
1167: 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);
1169: /* do not increase object state */
1170: if (mat == bA->m[idxm][jdxm]) return(0);
1172: PetscObjectReference((PetscObject)mat);
1173: MatDestroy(&bA->m[idxm][jdxm]);
1174: bA->m[idxm][jdxm] = mat;
1175: PetscObjectStateIncrease((PetscObject)A);
1176: MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1177: A->nonzerostate++;
1178: return(0);
1179: }
1181: /*@
1182: MatNestSetSubMat - Set a single submatrix in the nest matrix.
1184: Logically collective on the submatrix communicator
1186: Input Parameters:
1187: + A - nest matrix
1188: . idxm - index of the matrix within the nest matrix
1189: . jdxm - index of the matrix within the nest matrix
1190: - sub - matrix at index idxm,jdxm within the nest matrix
1192: Notes:
1193: The new submatrix must have the same size and communicator as that block of the nest.
1195: This increments the reference count of the submatrix.
1197: Level: developer
1199: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1200: MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1201: @*/
1202: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1203: {
1207: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1208: return(0);
1209: }
1211: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1212: {
1213: Mat_Nest *bA = (Mat_Nest*)A->data;
1216: if (M) *M = bA->nr;
1217: if (N) *N = bA->nc;
1218: if (mat) *mat = bA->m;
1219: return(0);
1220: }
1222: /*@C
1223: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
1225: Not collective
1227: Input Parameter:
1228: . A - nest matrix
1230: Output Parameters:
1231: + M - number of rows in the nest matrix
1232: . N - number of cols in the nest matrix
1233: - mat - 2d array of matrices
1235: Notes:
1237: The user should not free the array mat.
1239: In Fortran, this routine has a calling sequence
1240: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1241: where the space allocated for the optional argument mat is assumed large enough (if provided).
1243: Level: developer
1245: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1246: MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1247: @*/
1248: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1249: {
1253: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1254: return(0);
1255: }
1257: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1258: {
1259: Mat_Nest *bA = (Mat_Nest*)A->data;
1262: if (M) *M = bA->nr;
1263: if (N) *N = bA->nc;
1264: return(0);
1265: }
1267: /*@
1268: MatNestGetSize - Returns the size of the nest matrix.
1270: Not collective
1272: Input Parameter:
1273: . A - nest matrix
1275: Output Parameters:
1276: + M - number of rows in the nested mat
1277: - N - number of cols in the nested mat
1279: Notes:
1281: Level: developer
1283: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1284: MatNestGetISs()
1285: @*/
1286: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1287: {
1291: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1292: return(0);
1293: }
1295: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1296: {
1297: Mat_Nest *vs = (Mat_Nest*)A->data;
1298: PetscInt i;
1301: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1302: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1303: return(0);
1304: }
1306: /*@C
1307: MatNestGetISs - Returns the index sets partitioning the row and column spaces
1309: Not collective
1311: Input Parameter:
1312: . A - nest matrix
1314: Output Parameters:
1315: + rows - array of row index sets
1316: - cols - array of column index sets
1318: Level: advanced
1320: Notes:
1321: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1323: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1324: MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1325: @*/
1326: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1327: {
1332: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1333: return(0);
1334: }
1336: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1337: {
1338: Mat_Nest *vs = (Mat_Nest*)A->data;
1339: PetscInt i;
1342: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1343: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1344: return(0);
1345: }
1347: /*@C
1348: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1350: Not collective
1352: Input Parameter:
1353: . A - nest matrix
1355: Output Parameters:
1356: + rows - array of row index sets (or NULL to ignore)
1357: - cols - array of column index sets (or NULL to ignore)
1359: Level: advanced
1361: Notes:
1362: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1364: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1365: MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1366: @*/
1367: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1368: {
1373: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1374: return(0);
1375: }
1377: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1378: {
1380: PetscBool flg;
1383: PetscStrcmp(vtype,VECNEST,&flg);
1384: /* In reality, this only distinguishes VECNEST and "other" */
1385: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1386: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1387: return(0);
1388: }
1390: /*@C
1391: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1393: Not collective
1395: Input Parameters:
1396: + A - nest matrix
1397: - vtype - type to use for creating vectors
1399: Notes:
1401: Level: developer
1403: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1404: @*/
1405: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1406: {
1410: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1411: return(0);
1412: }
1414: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1415: {
1416: Mat_Nest *s = (Mat_Nest*)A->data;
1417: PetscInt i,j,m,n,M,N;
1419: PetscBool cong;
1422: MatReset_Nest(A);
1424: s->nr = nr;
1425: s->nc = nc;
1427: /* Create space for submatrices */
1428: PetscMalloc1(nr,&s->m);
1429: for (i=0; i<nr; i++) {
1430: PetscMalloc1(nc,&s->m[i]);
1431: }
1432: for (i=0; i<nr; i++) {
1433: for (j=0; j<nc; j++) {
1434: s->m[i][j] = a[i*nc+j];
1435: if (a[i*nc+j]) {
1436: PetscObjectReference((PetscObject)a[i*nc+j]);
1437: }
1438: }
1439: }
1441: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1443: PetscMalloc1(nr,&s->row_len);
1444: PetscMalloc1(nc,&s->col_len);
1445: for (i=0; i<nr; i++) s->row_len[i]=-1;
1446: for (j=0; j<nc; j++) s->col_len[j]=-1;
1448: PetscCalloc1(nr*nc,&s->nnzstate);
1449: for (i=0; i<nr; i++) {
1450: for (j=0; j<nc; j++) {
1451: if (s->m[i][j]) {
1452: MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1453: }
1454: }
1455: }
1457: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1459: PetscLayoutSetSize(A->rmap,M);
1460: PetscLayoutSetLocalSize(A->rmap,m);
1461: PetscLayoutSetSize(A->cmap,N);
1462: PetscLayoutSetLocalSize(A->cmap,n);
1464: PetscLayoutSetUp(A->rmap);
1465: PetscLayoutSetUp(A->cmap);
1467: /* disable operations that are not supported for non-square matrices,
1468: or matrices for which is_row != is_col */
1469: MatHasCongruentLayouts(A,&cong);
1470: if (cong && nr != nc) cong = PETSC_FALSE;
1471: if (cong) {
1472: for (i = 0; cong && i < nr; i++) {
1473: ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1474: }
1475: }
1476: if (!cong) {
1477: A->ops->missingdiagonal = NULL;
1478: A->ops->getdiagonal = NULL;
1479: A->ops->shift = NULL;
1480: A->ops->diagonalset = NULL;
1481: }
1483: PetscCalloc2(nr,&s->left,nc,&s->right);
1484: PetscObjectStateIncrease((PetscObject)A);
1485: A->nonzerostate++;
1486: return(0);
1487: }
1489: /*@
1490: MatNestSetSubMats - Sets the nested submatrices
1492: Collective on Mat
1494: Input Parameters:
1495: + A - nested matrix
1496: . nr - number of nested row blocks
1497: . is_row - index sets for each nested row block, or NULL to make contiguous
1498: . nc - number of nested column blocks
1499: . is_col - index sets for each nested column block, or NULL to make contiguous
1500: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1502: Notes: this always resets any submatrix information previously set
1504: Level: advanced
1506: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1507: @*/
1508: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1509: {
1511: PetscInt i;
1515: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1516: if (nr && is_row) {
1519: }
1520: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1521: if (nc && is_col) {
1524: }
1526: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1527: return(0);
1528: }
1530: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1531: {
1533: PetscBool flg;
1534: PetscInt i,j,m,mi,*ix;
1537: *ltog = NULL;
1538: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1539: if (islocal[i]) {
1540: ISGetLocalSize(islocal[i],&mi);
1541: flg = PETSC_TRUE; /* We found a non-trivial entry */
1542: } else {
1543: ISGetLocalSize(isglobal[i],&mi);
1544: }
1545: m += mi;
1546: }
1547: if (!flg) return(0);
1549: PetscMalloc1(m,&ix);
1550: for (i=0,m=0; i<n; i++) {
1551: ISLocalToGlobalMapping smap = NULL;
1552: Mat sub = NULL;
1553: PetscSF sf;
1554: PetscLayout map;
1555: const PetscInt *ix2;
1557: if (!colflg) {
1558: MatNestFindNonzeroSubMatRow(A,i,&sub);
1559: } else {
1560: MatNestFindNonzeroSubMatCol(A,i,&sub);
1561: }
1562: if (sub) {
1563: if (!colflg) {
1564: MatGetLocalToGlobalMapping(sub,&smap,NULL);
1565: } else {
1566: MatGetLocalToGlobalMapping(sub,NULL,&smap);
1567: }
1568: }
1569: /*
1570: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1571: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1572: */
1573: ISGetIndices(isglobal[i],&ix2);
1574: if (islocal[i]) {
1575: PetscInt *ilocal,*iremote;
1576: PetscInt mil,nleaves;
1578: ISGetLocalSize(islocal[i],&mi);
1579: if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1580: for (j=0; j<mi; j++) ix[m+j] = j;
1581: ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);
1583: /* PetscSFSetGraphLayout does not like negative indices */
1584: PetscMalloc2(mi,&ilocal,mi,&iremote);
1585: for (j=0, nleaves = 0; j<mi; j++) {
1586: if (ix[m+j] < 0) continue;
1587: ilocal[nleaves] = j;
1588: iremote[nleaves] = ix[m+j];
1589: nleaves++;
1590: }
1591: ISGetLocalSize(isglobal[i],&mil);
1592: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1593: PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1594: PetscLayoutSetLocalSize(map,mil);
1595: PetscLayoutSetUp(map);
1596: PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1597: PetscLayoutDestroy(&map);
1598: PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1599: PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1600: PetscSFDestroy(&sf);
1601: PetscFree2(ilocal,iremote);
1602: } else {
1603: ISGetLocalSize(isglobal[i],&mi);
1604: for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1605: }
1606: ISRestoreIndices(isglobal[i],&ix2);
1607: m += mi;
1608: }
1609: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1610: return(0);
1611: }
1613: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1614: /*
1615: nprocessors = NP
1616: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1617: proc 0: => (g_0,h_0,)
1618: proc 1: => (g_1,h_1,)
1619: ...
1620: proc nprocs-1: => (g_NP-1,h_NP-1,)
1622: proc 0: proc 1: proc nprocs-1:
1623: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1625: proc 0:
1626: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1627: proc 1:
1628: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1630: proc NP-1:
1631: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1632: */
1633: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1634: {
1635: Mat_Nest *vs = (Mat_Nest*)A->data;
1636: PetscInt i,j,offset,n,nsum,bs;
1638: Mat sub = NULL;
1641: PetscMalloc1(nr,&vs->isglobal.row);
1642: PetscMalloc1(nc,&vs->isglobal.col);
1643: if (is_row) { /* valid IS is passed in */
1644: /* refs on is[] are incremented */
1645: for (i=0; i<vs->nr; i++) {
1646: PetscObjectReference((PetscObject)is_row[i]);
1648: vs->isglobal.row[i] = is_row[i];
1649: }
1650: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1651: nsum = 0;
1652: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1653: MatNestFindNonzeroSubMatRow(A,i,&sub);
1654: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1655: MatGetLocalSize(sub,&n,NULL);
1656: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1657: nsum += n;
1658: }
1659: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1660: offset -= nsum;
1661: for (i=0; i<vs->nr; i++) {
1662: MatNestFindNonzeroSubMatRow(A,i,&sub);
1663: MatGetLocalSize(sub,&n,NULL);
1664: MatGetBlockSizes(sub,&bs,NULL);
1665: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1666: ISSetBlockSize(vs->isglobal.row[i],bs);
1667: offset += n;
1668: }
1669: }
1671: if (is_col) { /* valid IS is passed in */
1672: /* refs on is[] are incremented */
1673: for (j=0; j<vs->nc; j++) {
1674: PetscObjectReference((PetscObject)is_col[j]);
1676: vs->isglobal.col[j] = is_col[j];
1677: }
1678: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1679: offset = A->cmap->rstart;
1680: nsum = 0;
1681: for (j=0; j<vs->nc; j++) {
1682: MatNestFindNonzeroSubMatCol(A,j,&sub);
1683: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1684: MatGetLocalSize(sub,NULL,&n);
1685: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1686: nsum += n;
1687: }
1688: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1689: offset -= nsum;
1690: for (j=0; j<vs->nc; j++) {
1691: MatNestFindNonzeroSubMatCol(A,j,&sub);
1692: MatGetLocalSize(sub,NULL,&n);
1693: MatGetBlockSizes(sub,NULL,&bs);
1694: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1695: ISSetBlockSize(vs->isglobal.col[j],bs);
1696: offset += n;
1697: }
1698: }
1700: /* Set up the local ISs */
1701: PetscMalloc1(vs->nr,&vs->islocal.row);
1702: PetscMalloc1(vs->nc,&vs->islocal.col);
1703: for (i=0,offset=0; i<vs->nr; i++) {
1704: IS isloc;
1705: ISLocalToGlobalMapping rmap = NULL;
1706: PetscInt nlocal,bs;
1707: MatNestFindNonzeroSubMatRow(A,i,&sub);
1708: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1709: if (rmap) {
1710: MatGetBlockSizes(sub,&bs,NULL);
1711: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1712: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1713: ISSetBlockSize(isloc,bs);
1714: } else {
1715: nlocal = 0;
1716: isloc = NULL;
1717: }
1718: vs->islocal.row[i] = isloc;
1719: offset += nlocal;
1720: }
1721: for (i=0,offset=0; i<vs->nc; i++) {
1722: IS isloc;
1723: ISLocalToGlobalMapping cmap = NULL;
1724: PetscInt nlocal,bs;
1725: MatNestFindNonzeroSubMatCol(A,i,&sub);
1726: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1727: if (cmap) {
1728: MatGetBlockSizes(sub,NULL,&bs);
1729: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1730: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1731: ISSetBlockSize(isloc,bs);
1732: } else {
1733: nlocal = 0;
1734: isloc = NULL;
1735: }
1736: vs->islocal.col[i] = isloc;
1737: offset += nlocal;
1738: }
1740: /* Set up the aggregate ISLocalToGlobalMapping */
1741: {
1742: ISLocalToGlobalMapping rmap,cmap;
1743: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1744: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1745: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1746: ISLocalToGlobalMappingDestroy(&rmap);
1747: ISLocalToGlobalMappingDestroy(&cmap);
1748: }
1750: if (PetscDefined(USE_DEBUG)) {
1751: for (i=0; i<vs->nr; i++) {
1752: for (j=0; j<vs->nc; j++) {
1753: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1754: Mat B = vs->m[i][j];
1755: if (!B) continue;
1756: MatGetSize(B,&M,&N);
1757: MatGetLocalSize(B,&m,&n);
1758: ISGetSize(vs->isglobal.row[i],&Mi);
1759: ISGetSize(vs->isglobal.col[j],&Ni);
1760: ISGetLocalSize(vs->isglobal.row[i],&mi);
1761: ISGetLocalSize(vs->isglobal.col[j],&ni);
1762: 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);
1763: 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);
1764: }
1765: }
1766: }
1768: /* Set A->assembled if all non-null blocks are currently assembled */
1769: for (i=0; i<vs->nr; i++) {
1770: for (j=0; j<vs->nc; j++) {
1771: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1772: }
1773: }
1774: A->assembled = PETSC_TRUE;
1775: return(0);
1776: }
1778: /*@C
1779: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1781: Collective on Mat
1783: Input Parameters:
1784: + comm - Communicator for the new Mat
1785: . nr - number of nested row blocks
1786: . is_row - index sets for each nested row block, or NULL to make contiguous
1787: . nc - number of nested column blocks
1788: . is_col - index sets for each nested column block, or NULL to make contiguous
1789: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1791: Output Parameter:
1792: . B - new matrix
1794: Level: advanced
1796: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1797: MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1798: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1799: @*/
1800: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1801: {
1802: Mat A;
1806: *B = NULL;
1807: MatCreate(comm,&A);
1808: MatSetType(A,MATNEST);
1809: A->preallocated = PETSC_TRUE;
1810: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1811: *B = A;
1812: return(0);
1813: }
1815: PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1816: {
1817: Mat_Nest *nest = (Mat_Nest*)A->data;
1818: Mat *trans;
1819: PetscScalar **avv;
1820: PetscScalar *vv;
1821: PetscInt **aii,**ajj;
1822: PetscInt *ii,*jj,*ci;
1823: PetscInt nr,nc,nnz,i,j;
1824: PetscBool done;
1828: MatGetSize(A,&nr,&nc);
1829: if (reuse == MAT_REUSE_MATRIX) {
1830: PetscInt rnr;
1832: MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1833: if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1834: if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1835: MatSeqAIJGetArray(*newmat,&vv);
1836: }
1837: /* extract CSR for nested SeqAIJ matrices */
1838: nnz = 0;
1839: PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1840: for (i=0; i<nest->nr; ++i) {
1841: for (j=0; j<nest->nc; ++j) {
1842: Mat B = nest->m[i][j];
1843: if (B) {
1844: PetscScalar *naa;
1845: PetscInt *nii,*njj,nnr;
1846: PetscBool istrans;
1848: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1849: if (istrans) {
1850: Mat Bt;
1852: MatTransposeGetMat(B,&Bt);
1853: MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1854: B = trans[i*nest->nc+j];
1855: }
1856: MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1857: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1858: MatSeqAIJGetArray(B,&naa);
1859: nnz += nii[nnr];
1861: aii[i*nest->nc+j] = nii;
1862: ajj[i*nest->nc+j] = njj;
1863: avv[i*nest->nc+j] = naa;
1864: }
1865: }
1866: }
1867: if (reuse != MAT_REUSE_MATRIX) {
1868: PetscMalloc1(nr+1,&ii);
1869: PetscMalloc1(nnz,&jj);
1870: PetscMalloc1(nnz,&vv);
1871: } else {
1872: if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1873: }
1875: /* new row pointer */
1876: PetscArrayzero(ii,nr+1);
1877: for (i=0; i<nest->nr; ++i) {
1878: PetscInt ncr,rst;
1880: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1881: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1882: for (j=0; j<nest->nc; ++j) {
1883: if (aii[i*nest->nc+j]) {
1884: PetscInt *nii = aii[i*nest->nc+j];
1885: PetscInt ir;
1887: for (ir=rst; ir<ncr+rst; ++ir) {
1888: ii[ir+1] += nii[1]-nii[0];
1889: nii++;
1890: }
1891: }
1892: }
1893: }
1894: for (i=0; i<nr; i++) ii[i+1] += ii[i];
1896: /* construct CSR for the new matrix */
1897: PetscCalloc1(nr,&ci);
1898: for (i=0; i<nest->nr; ++i) {
1899: PetscInt ncr,rst;
1901: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1902: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1903: for (j=0; j<nest->nc; ++j) {
1904: if (aii[i*nest->nc+j]) {
1905: PetscScalar *nvv = avv[i*nest->nc+j];
1906: PetscInt *nii = aii[i*nest->nc+j];
1907: PetscInt *njj = ajj[i*nest->nc+j];
1908: PetscInt ir,cst;
1910: ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1911: for (ir=rst; ir<ncr+rst; ++ir) {
1912: PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1914: for (ij=0;ij<rsize;ij++) {
1915: jj[ist+ij] = *njj+cst;
1916: vv[ist+ij] = *nvv;
1917: njj++;
1918: nvv++;
1919: }
1920: ci[ir] += rsize;
1921: nii++;
1922: }
1923: }
1924: }
1925: }
1926: PetscFree(ci);
1928: /* restore info */
1929: for (i=0; i<nest->nr; ++i) {
1930: for (j=0; j<nest->nc; ++j) {
1931: Mat B = nest->m[i][j];
1932: if (B) {
1933: PetscInt nnr = 0, k = i*nest->nc+j;
1935: B = (trans[k] ? trans[k] : B);
1936: MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1937: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1938: MatSeqAIJRestoreArray(B,&avv[k]);
1939: MatDestroy(&trans[k]);
1940: }
1941: }
1942: }
1943: PetscFree4(aii,ajj,avv,trans);
1945: /* finalize newmat */
1946: if (reuse == MAT_INITIAL_MATRIX) {
1947: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1948: } else if (reuse == MAT_INPLACE_MATRIX) {
1949: Mat B;
1951: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1952: MatHeaderReplace(A,&B);
1953: }
1954: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1955: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1956: {
1957: Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1958: a->free_a = PETSC_TRUE;
1959: a->free_ij = PETSC_TRUE;
1960: }
1961: return(0);
1962: }
1964: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y,PetscScalar a,Mat X)
1965: {
1966: Mat_Nest *nest = (Mat_Nest*)X->data;
1967: PetscInt i,j,k,rstart;
1968: PetscBool flg;
1972: /* Fill by row */
1973: for (j=0; j<nest->nc; ++j) {
1974: /* Using global column indices and ISAllGather() is not scalable. */
1975: IS bNis;
1976: PetscInt bN;
1977: const PetscInt *bNindices;
1978: ISAllGather(nest->isglobal.col[j], &bNis);
1979: ISGetSize(bNis,&bN);
1980: ISGetIndices(bNis,&bNindices);
1981: for (i=0; i<nest->nr; ++i) {
1982: Mat B,D=NULL;
1983: PetscInt bm, br;
1984: const PetscInt *bmindices;
1985: B = nest->m[i][j];
1986: if (!B) continue;
1987: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&flg);
1988: if (flg) {
1989: PetscTryMethod(B,"MatTransposeGetMat_C",(Mat,Mat*),(B,&D));
1990: PetscTryMethod(B,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(B,&D));
1991: MatConvert(B,((PetscObject)D)->type_name,MAT_INITIAL_MATRIX,&D);
1992: B = D;
1993: }
1994: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQSBAIJ,MATMPISBAIJ,"");
1995: if (flg) {
1996: if (D) {
1997: MatConvert(D,MATBAIJ,MAT_INPLACE_MATRIX,&D);
1998: } else {
1999: MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&D);
2000: }
2001: B = D;
2002: }
2003: ISGetLocalSize(nest->isglobal.row[i],&bm);
2004: ISGetIndices(nest->isglobal.row[i],&bmindices);
2005: MatGetOwnershipRange(B,&rstart,NULL);
2006: for (br = 0; br < bm; ++br) {
2007: PetscInt row = bmindices[br], brncols, *cols;
2008: const PetscInt *brcols;
2009: const PetscScalar *brcoldata;
2010: PetscScalar *vals = NULL;
2011: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2012: PetscMalloc1(brncols,&cols);
2013: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2014: /*
2015: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2016: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2017: */
2018: if (a != 1.0) {
2019: PetscMalloc1(brncols,&vals);
2020: for (k=0; k<brncols; k++) vals[k] = a * brcoldata[k];
2021: MatSetValues(Y,1,&row,brncols,cols,vals,ADD_VALUES);
2022: PetscFree(vals);
2023: } else {
2024: MatSetValues(Y,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2025: }
2026: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2027: PetscFree(cols);
2028: }
2029: if (D) {
2030: MatDestroy(&D);
2031: }
2032: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2033: }
2034: ISRestoreIndices(bNis,&bNindices);
2035: ISDestroy(&bNis);
2036: }
2037: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
2038: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
2039: return(0);
2040: }
2042: PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2043: {
2045: Mat_Nest *nest = (Mat_Nest*)A->data;
2046: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart,cstart,cend;
2047: PetscMPIInt size;
2048: Mat C;
2051: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
2052: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
2053: PetscInt nf;
2054: PetscBool fast;
2056: PetscStrcmp(newtype,MATAIJ,&fast);
2057: if (!fast) {
2058: PetscStrcmp(newtype,MATSEQAIJ,&fast);
2059: }
2060: for (i=0; i<nest->nr && fast; ++i) {
2061: for (j=0; j<nest->nc && fast; ++j) {
2062: Mat B = nest->m[i][j];
2063: if (B) {
2064: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
2065: if (!fast) {
2066: PetscBool istrans;
2068: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
2069: if (istrans) {
2070: Mat Bt;
2072: MatTransposeGetMat(B,&Bt);
2073: PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
2074: }
2075: }
2076: }
2077: }
2078: }
2079: for (i=0, nf=0; i<nest->nr && fast; ++i) {
2080: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
2081: if (fast) {
2082: PetscInt f,s;
2084: ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
2085: if (f != nf || s != 1) { fast = PETSC_FALSE; }
2086: else {
2087: ISGetSize(nest->isglobal.row[i],&f);
2088: nf += f;
2089: }
2090: }
2091: }
2092: for (i=0, nf=0; i<nest->nc && fast; ++i) {
2093: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
2094: if (fast) {
2095: PetscInt f,s;
2097: ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
2098: if (f != nf || s != 1) { fast = PETSC_FALSE; }
2099: else {
2100: ISGetSize(nest->isglobal.col[i],&f);
2101: nf += f;
2102: }
2103: }
2104: }
2105: if (fast) {
2106: MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
2107: return(0);
2108: }
2109: }
2110: MatGetSize(A,&M,&N);
2111: MatGetLocalSize(A,&m,&n);
2112: MatGetOwnershipRangeColumn(A,&cstart,&cend);
2113: if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2114: else {
2115: MatCreate(PetscObjectComm((PetscObject)A),&C);
2116: MatSetType(C,newtype);
2117: MatSetSizes(C,m,n,M,N);
2118: }
2119: PetscMalloc1(2*m,&dnnz);
2120: onnz = dnnz + m;
2121: for (k=0; k<m; k++) {
2122: dnnz[k] = 0;
2123: onnz[k] = 0;
2124: }
2125: for (j=0; j<nest->nc; ++j) {
2126: IS bNis;
2127: PetscInt bN;
2128: const PetscInt *bNindices;
2129: /* Using global column indices and ISAllGather() is not scalable. */
2130: ISAllGather(nest->isglobal.col[j], &bNis);
2131: ISGetSize(bNis, &bN);
2132: ISGetIndices(bNis,&bNindices);
2133: for (i=0; i<nest->nr; ++i) {
2134: PetscSF bmsf;
2135: PetscSFNode *iremote;
2136: Mat B;
2137: PetscInt bm, *sub_dnnz,*sub_onnz, br;
2138: const PetscInt *bmindices;
2139: B = nest->m[i][j];
2140: if (!B) continue;
2141: ISGetLocalSize(nest->isglobal.row[i],&bm);
2142: ISGetIndices(nest->isglobal.row[i],&bmindices);
2143: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2144: PetscMalloc1(bm,&iremote);
2145: PetscMalloc1(bm,&sub_dnnz);
2146: PetscMalloc1(bm,&sub_onnz);
2147: for (k = 0; k < bm; ++k) {
2148: sub_dnnz[k] = 0;
2149: sub_onnz[k] = 0;
2150: }
2151: /*
2152: Locate the owners for all of the locally-owned global row indices for this row block.
2153: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2154: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2155: */
2156: MatGetOwnershipRange(B,&rstart,NULL);
2157: for (br = 0; br < bm; ++br) {
2158: PetscInt row = bmindices[br], brncols, col;
2159: const PetscInt *brcols;
2160: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2161: PetscMPIInt rowowner = 0;
2162: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2163: /* how many roots */
2164: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2165: /* get nonzero pattern */
2166: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2167: for (k=0; k<brncols; k++) {
2168: col = bNindices[brcols[k]];
2169: if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2170: sub_dnnz[br]++;
2171: } else {
2172: sub_onnz[br]++;
2173: }
2174: }
2175: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2176: }
2177: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2178: /* bsf will have to take care of disposing of bedges. */
2179: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2180: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2181: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2182: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2183: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2184: PetscFree(sub_dnnz);
2185: PetscFree(sub_onnz);
2186: PetscSFDestroy(&bmsf);
2187: }
2188: ISRestoreIndices(bNis,&bNindices);
2189: ISDestroy(&bNis);
2190: }
2191: /* Resize preallocation if overestimated */
2192: for (i=0;i<m;i++) {
2193: dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2194: onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2195: }
2196: MatSeqAIJSetPreallocation(C,0,dnnz);
2197: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2198: PetscFree(dnnz);
2199: MatAXPY_Dense_Nest(C,1.0,A);
2200: if (reuse == MAT_INPLACE_MATRIX) {
2201: MatHeaderReplace(A,&C);
2202: } else *newmat = C;
2203: return(0);
2204: }
2206: PetscErrorCode MatConvert_Nest_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2207: {
2208: Mat B;
2209: PetscInt m,n,M,N;
2213: MatGetSize(A,&M,&N);
2214: MatGetLocalSize(A,&m,&n);
2215: if (reuse == MAT_REUSE_MATRIX) {
2216: B = *newmat;
2217: MatZeroEntries(B);
2218: } else {
2219: MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,&B);
2220: }
2221: MatAXPY_Dense_Nest(B,1.0,A);
2222: if (reuse == MAT_INPLACE_MATRIX) {
2223: MatHeaderReplace(A,&B);
2224: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2225: return(0);
2226: }
2228: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2229: {
2230: Mat_Nest *bA = (Mat_Nest*)mat->data;
2231: MatOperation opAdd;
2232: PetscInt i,j,nr = bA->nr,nc = bA->nc;
2233: PetscBool flg;
2237: *has = PETSC_FALSE;
2238: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2239: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2240: for (j=0; j<nc; j++) {
2241: for (i=0; i<nr; i++) {
2242: if (!bA->m[i][j]) continue;
2243: MatHasOperation(bA->m[i][j],opAdd,&flg);
2244: if (!flg) return(0);
2245: }
2246: }
2247: }
2248: if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2249: return(0);
2250: }
2252: /*MC
2253: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2255: Level: intermediate
2257: Notes:
2258: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2259: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2260: It is usually used with DMComposite and DMCreateMatrix()
2262: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2263: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2264: than the nest matrix.
2266: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2267: VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2268: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2269: M*/
2270: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2271: {
2272: Mat_Nest *s;
2276: PetscNewLog(A,&s);
2277: A->data = (void*)s;
2279: s->nr = -1;
2280: s->nc = -1;
2281: s->m = NULL;
2282: s->splitassembly = PETSC_FALSE;
2284: PetscMemzero(A->ops,sizeof(*A->ops));
2286: A->ops->mult = MatMult_Nest;
2287: A->ops->multadd = MatMultAdd_Nest;
2288: A->ops->multtranspose = MatMultTranspose_Nest;
2289: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2290: A->ops->transpose = MatTranspose_Nest;
2291: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2292: A->ops->assemblyend = MatAssemblyEnd_Nest;
2293: A->ops->zeroentries = MatZeroEntries_Nest;
2294: A->ops->copy = MatCopy_Nest;
2295: A->ops->axpy = MatAXPY_Nest;
2296: A->ops->duplicate = MatDuplicate_Nest;
2297: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2298: A->ops->destroy = MatDestroy_Nest;
2299: A->ops->view = MatView_Nest;
2300: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2301: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2302: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2303: A->ops->getdiagonal = MatGetDiagonal_Nest;
2304: A->ops->diagonalscale = MatDiagonalScale_Nest;
2305: A->ops->scale = MatScale_Nest;
2306: A->ops->shift = MatShift_Nest;
2307: A->ops->diagonalset = MatDiagonalSet_Nest;
2308: A->ops->setrandom = MatSetRandom_Nest;
2309: A->ops->hasoperation = MatHasOperation_Nest;
2310: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2312: A->spptr = NULL;
2313: A->assembled = PETSC_FALSE;
2315: /* expose Nest api's */
2316: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2317: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2318: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2319: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
2320: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
2321: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2322: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
2323: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2324: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2325: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2326: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2327: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);
2328: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",MatConvert_Nest_Dense);
2329: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",MatConvert_Nest_Dense);
2330: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2331: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2332: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);
2334: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2335: return(0);
2336: }