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: MatDenseGetArrayWrite(C,&carray);
126: for (i=0; i<nr; i++) {
127: ISGetSize(bA->isglobal.row[i],&M);
128: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
129: MatDenseSetLDA(viewC,ldc);
130: for (j=0; j<nc; j++) {
131: if (!bA->m[i][j]) continue;
132: ISGetSize(bA->isglobal.col[j],&M);
133: MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
134: MatDenseSetLDA(viewB,ldb);
136: /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
137: workC = contents->workC[i*nc + j];
138: productB = workC->product->B;
139: workC->product->B = viewB; /* use newly created dense matrix viewB */
140: MatProductNumeric(workC);
141: MatDestroy(&viewB);
142: workC->product->B = productB; /* resume original B */
144: /* C[i] <- workC + C[i] */
145: MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
146: }
147: MatDestroy(&viewC);
148: }
149: MatDenseRestoreArrayWrite(C,&carray);
150: MatDenseRestoreArrayRead(B,&barray);
152: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
154: return(0);
155: }
157: PetscErrorCode MatNest_DenseDestroy(void *ctx)
158: {
159: Nest_Dense *contents = (Nest_Dense*)ctx;
160: PetscInt i;
164: PetscFree(contents->tarray);
165: for (i=0; i<contents->k; i++) {
166: MatDestroy(contents->workC + i);
167: }
168: PetscFree3(contents->dm,contents->dn,contents->workC);
169: PetscFree(contents);
170: return(0);
171: }
173: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
174: {
175: Mat_Nest *bA;
176: Mat viewB,workC;
177: const PetscScalar *barray;
178: PetscInt i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
179: Nest_Dense *contents=NULL;
180: PetscBool cisdense;
181: PetscErrorCode ierr;
182: Mat A,B;
183: PetscReal fill;
186: MatCheckProduct(C,4);
187: if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
188: A = C->product->A;
189: B = C->product->B;
190: fill = C->product->fill;
191: bA = (Mat_Nest*)A->data;
192: nr = bA->nr;
193: nc = bA->nc;
194: MatGetLocalSize(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: {
438: MatReset_Nest(A);
439: 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: }
622: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A,PetscInt i,PetscInt j,Mat *B)
623: {
624: Mat_Nest *vs = (Mat_Nest*)A->data;
625: PetscInt lr,lc;
629: MatCreate(PetscObjectComm((PetscObject)A),B);
630: ISGetLocalSize(vs->isglobal.row[i],&lr);
631: ISGetLocalSize(vs->isglobal.col[j],&lc);
632: MatSetSizes(*B,lr,lc,PETSC_DECIDE,PETSC_DECIDE);
633: MatSetType(*B,MATAIJ);
634: MatSeqAIJSetPreallocation(*B,0,NULL);
635: MatMPIAIJSetPreallocation(*B,0,NULL,0,NULL);
636: MatSetUp(*B);
637: MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
638: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
639: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
640: return(0);
641: }
643: static PetscErrorCode MatNestGetBlock_Private(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *B)
644: {
645: Mat_Nest *vs = (Mat_Nest*)A->data;
646: Mat *a;
647: PetscInt i,j,k,l,nr=rend-rbegin,nc=cend-cbegin;
648: char keyname[256];
649: PetscBool *b;
650: PetscBool flg;
654: *B = NULL;
655: PetscSNPrintf(keyname,sizeof(keyname),"NestBlock_%D-%Dx%D-%D",rbegin,rend,cbegin,cend);
656: PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
657: if (*B) return(0);
659: PetscMalloc2(nr*nc,&a,nr*nc,&b);
660: for (i=0; i<nr; i++) {
661: for (j=0; j<nc; j++) {
662: a[i*nc + j] = vs->m[rbegin+i][cbegin+j];
663: b[i*nc + j] = PETSC_FALSE;
664: }
665: }
666: if (nc!=vs->nc&&nr!=vs->nr) {
667: for (i=0; i<nr; i++) {
668: for (j=0; j<nc; j++) {
669: flg = PETSC_FALSE;
670: for (k=0; (k<nr&&!flg); k++) {
671: if (a[j + k*nc]) flg = PETSC_TRUE;
672: }
673: if (flg) {
674: flg = PETSC_FALSE;
675: for (l=0; (l<nc&&!flg); l++) {
676: if (a[i*nc + l]) flg = PETSC_TRUE;
677: }
678: }
679: if (!flg) {
680: b[i*nc + j] = PETSC_TRUE;
681: MatNestFillEmptyMat_Private(A,rbegin+i,cbegin+j,a + i*nc + j);
682: }
683: }
684: }
685: }
686: MatCreateNest(PetscObjectComm((PetscObject)A),nr,nr!=vs->nr?NULL:vs->isglobal.row,nc,nc!=vs->nc?NULL:vs->isglobal.col,a,B);
687: for (i=0; i<nr; i++) {
688: for (j=0; j<nc; j++) {
689: if (b[i*nc + j]) {
690: MatDestroy(a + i*nc + j);
691: }
692: }
693: }
694: PetscFree2(a,b);
695: (*B)->assembled = A->assembled;
696: PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
697: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
698: return(0);
699: }
701: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
702: {
703: Mat_Nest *vs = (Mat_Nest*)A->data;
704: PetscInt rbegin,rend,cbegin,cend;
708: MatNestFindISRange(A,vs->nr,is->row,isrow,&rbegin,&rend);
709: MatNestFindISRange(A,vs->nc,is->col,iscol,&cbegin,&cend);
710: if (rend == rbegin + 1 && cend == cbegin + 1) {
711: if (!vs->m[rbegin][cbegin]) {
712: MatNestFillEmptyMat_Private(A,rbegin,cbegin,vs->m[rbegin] + cbegin);
713: }
714: *B = vs->m[rbegin][cbegin];
715: } else if (rbegin != -1 && cbegin != -1) {
716: MatNestGetBlock_Private(A,rbegin,rend,cbegin,cend,B);
717: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
718: return(0);
719: }
721: /*
722: TODO: This does not actually returns a submatrix we can modify
723: */
724: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
725: {
727: Mat_Nest *vs = (Mat_Nest*)A->data;
728: Mat sub;
731: MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
732: switch (reuse) {
733: case MAT_INITIAL_MATRIX:
734: if (sub) { PetscObjectReference((PetscObject)sub); }
735: *B = sub;
736: break;
737: case MAT_REUSE_MATRIX:
738: if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
739: break;
740: case MAT_IGNORE_MATRIX: /* Nothing to do */
741: break;
742: case MAT_INPLACE_MATRIX: /* Nothing to do */
743: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
744: }
745: return(0);
746: }
748: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
749: {
751: Mat_Nest *vs = (Mat_Nest*)A->data;
752: Mat sub;
755: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
756: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
757: if (sub) {PetscObjectReference((PetscObject)sub);}
758: *B = sub;
759: return(0);
760: }
762: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
763: {
765: Mat_Nest *vs = (Mat_Nest*)A->data;
766: Mat sub;
769: MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
770: if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
771: if (sub) {
772: if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
773: MatDestroy(B);
774: }
775: return(0);
776: }
778: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
779: {
780: Mat_Nest *bA = (Mat_Nest*)A->data;
781: PetscInt i;
785: for (i=0; i<bA->nr; i++) {
786: Vec bv;
787: VecGetSubVector(v,bA->isglobal.row[i],&bv);
788: if (bA->m[i][i]) {
789: MatGetDiagonal(bA->m[i][i],bv);
790: } else {
791: VecSet(bv,0.0);
792: }
793: VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
794: }
795: return(0);
796: }
798: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
799: {
800: Mat_Nest *bA = (Mat_Nest*)A->data;
801: Vec bl,*br;
802: PetscInt i,j;
806: PetscCalloc1(bA->nc,&br);
807: if (r) {
808: for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
809: }
810: bl = NULL;
811: for (i=0; i<bA->nr; i++) {
812: if (l) {
813: VecGetSubVector(l,bA->isglobal.row[i],&bl);
814: }
815: for (j=0; j<bA->nc; j++) {
816: if (bA->m[i][j]) {
817: MatDiagonalScale(bA->m[i][j],bl,br[j]);
818: }
819: }
820: if (l) {
821: VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
822: }
823: }
824: if (r) {
825: for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
826: }
827: PetscFree(br);
828: return(0);
829: }
831: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
832: {
833: Mat_Nest *bA = (Mat_Nest*)A->data;
834: PetscInt i,j;
838: for (i=0; i<bA->nr; i++) {
839: for (j=0; j<bA->nc; j++) {
840: if (bA->m[i][j]) {
841: MatScale(bA->m[i][j],a);
842: }
843: }
844: }
845: return(0);
846: }
848: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
849: {
850: Mat_Nest *bA = (Mat_Nest*)A->data;
851: PetscInt i;
853: PetscBool nnzstate = PETSC_FALSE;
856: for (i=0; i<bA->nr; i++) {
857: PetscObjectState subnnzstate = 0;
858: 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);
859: MatShift(bA->m[i][i],a);
860: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
861: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
862: bA->nnzstate[i*bA->nc+i] = subnnzstate;
863: }
864: if (nnzstate) A->nonzerostate++;
865: return(0);
866: }
868: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
869: {
870: Mat_Nest *bA = (Mat_Nest*)A->data;
871: PetscInt i;
873: PetscBool nnzstate = PETSC_FALSE;
876: for (i=0; i<bA->nr; i++) {
877: PetscObjectState subnnzstate = 0;
878: Vec bv;
879: VecGetSubVector(D,bA->isglobal.row[i],&bv);
880: if (bA->m[i][i]) {
881: MatDiagonalSet(bA->m[i][i],bv,is);
882: MatGetNonzeroState(bA->m[i][i],&subnnzstate);
883: }
884: VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
885: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
886: bA->nnzstate[i*bA->nc+i] = subnnzstate;
887: }
888: if (nnzstate) A->nonzerostate++;
889: return(0);
890: }
892: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
893: {
894: Mat_Nest *bA = (Mat_Nest*)A->data;
895: PetscInt i,j;
899: for (i=0; i<bA->nr; i++) {
900: for (j=0; j<bA->nc; j++) {
901: if (bA->m[i][j]) {
902: MatSetRandom(bA->m[i][j],rctx);
903: }
904: }
905: }
906: return(0);
907: }
909: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
910: {
911: Mat_Nest *bA = (Mat_Nest*)A->data;
912: Vec *L,*R;
913: MPI_Comm comm;
914: PetscInt i,j;
918: PetscObjectGetComm((PetscObject)A,&comm);
919: if (right) {
920: /* allocate R */
921: PetscMalloc1(bA->nc, &R);
922: /* Create the right vectors */
923: for (j=0; j<bA->nc; j++) {
924: for (i=0; i<bA->nr; i++) {
925: if (bA->m[i][j]) {
926: MatCreateVecs(bA->m[i][j],&R[j],NULL);
927: break;
928: }
929: }
930: if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
931: }
932: VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
933: /* hand back control to the nest vector */
934: for (j=0; j<bA->nc; j++) {
935: VecDestroy(&R[j]);
936: }
937: PetscFree(R);
938: }
940: if (left) {
941: /* allocate L */
942: PetscMalloc1(bA->nr, &L);
943: /* Create the left vectors */
944: for (i=0; i<bA->nr; i++) {
945: for (j=0; j<bA->nc; j++) {
946: if (bA->m[i][j]) {
947: MatCreateVecs(bA->m[i][j],NULL,&L[i]);
948: break;
949: }
950: }
951: if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
952: }
954: VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
955: for (i=0; i<bA->nr; i++) {
956: VecDestroy(&L[i]);
957: }
959: PetscFree(L);
960: }
961: return(0);
962: }
964: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
965: {
966: Mat_Nest *bA = (Mat_Nest*)A->data;
967: PetscBool isascii,viewSub = PETSC_FALSE;
968: PetscInt i,j;
972: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
973: if (isascii) {
975: PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
976: PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
977: PetscViewerASCIIPushTab(viewer);
978: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);
980: PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
981: for (i=0; i<bA->nr; i++) {
982: for (j=0; j<bA->nc; j++) {
983: MatType type;
984: char name[256] = "",prefix[256] = "";
985: PetscInt NR,NC;
986: PetscBool isNest = PETSC_FALSE;
988: if (!bA->m[i][j]) {
989: PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
990: continue;
991: }
992: MatGetSize(bA->m[i][j],&NR,&NC);
993: MatGetType(bA->m[i][j], &type);
994: if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
995: if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
996: PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);
998: PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);
1000: if (isNest || viewSub) {
1001: PetscViewerASCIIPushTab(viewer); /* push1 */
1002: MatView(bA->m[i][j],viewer);
1003: PetscViewerASCIIPopTab(viewer); /* pop1 */
1004: }
1005: }
1006: }
1007: PetscViewerASCIIPopTab(viewer); /* pop0 */
1008: }
1009: return(0);
1010: }
1012: static PetscErrorCode MatZeroEntries_Nest(Mat A)
1013: {
1014: Mat_Nest *bA = (Mat_Nest*)A->data;
1015: PetscInt i,j;
1019: for (i=0; i<bA->nr; i++) {
1020: for (j=0; j<bA->nc; j++) {
1021: if (!bA->m[i][j]) continue;
1022: MatZeroEntries(bA->m[i][j]);
1023: }
1024: }
1025: return(0);
1026: }
1028: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
1029: {
1030: Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
1031: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1033: PetscBool nnzstate = PETSC_FALSE;
1036: 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);
1037: for (i=0; i<nr; i++) {
1038: for (j=0; j<nc; j++) {
1039: PetscObjectState subnnzstate = 0;
1040: if (bA->m[i][j] && bB->m[i][j]) {
1041: MatCopy(bA->m[i][j],bB->m[i][j],str);
1042: } 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);
1043: MatGetNonzeroState(bB->m[i][j],&subnnzstate);
1044: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
1045: bB->nnzstate[i*nc+j] = subnnzstate;
1046: }
1047: }
1048: if (nnzstate) B->nonzerostate++;
1049: return(0);
1050: }
1052: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
1053: {
1054: Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
1055: PetscInt i,j,nr = bY->nr,nc = bY->nc;
1057: PetscBool nnzstate = PETSC_FALSE;
1060: 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);
1061: for (i=0; i<nr; i++) {
1062: for (j=0; j<nc; j++) {
1063: PetscObjectState subnnzstate = 0;
1064: if (bY->m[i][j] && bX->m[i][j]) {
1065: MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
1066: } else if (bX->m[i][j]) {
1067: Mat M;
1069: 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);
1070: MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
1071: MatNestSetSubMat(Y,i,j,M);
1072: MatDestroy(&M);
1073: }
1074: if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
1075: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1076: bY->nnzstate[i*nc+j] = subnnzstate;
1077: }
1078: }
1079: if (nnzstate) Y->nonzerostate++;
1080: return(0);
1081: }
1083: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1084: {
1085: Mat_Nest *bA = (Mat_Nest*)A->data;
1086: Mat *b;
1087: PetscInt i,j,nr = bA->nr,nc = bA->nc;
1091: PetscMalloc1(nr*nc,&b);
1092: for (i=0; i<nr; i++) {
1093: for (j=0; j<nc; j++) {
1094: if (bA->m[i][j]) {
1095: MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1096: } else {
1097: b[i*nc+j] = NULL;
1098: }
1099: }
1100: }
1101: MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1102: /* Give the new MatNest exclusive ownership */
1103: for (i=0; i<nr*nc; i++) {
1104: MatDestroy(&b[i]);
1105: }
1106: PetscFree(b);
1108: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1109: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1110: return(0);
1111: }
1113: /* nest api */
1114: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1115: {
1116: Mat_Nest *bA = (Mat_Nest*)A->data;
1119: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1120: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1121: *mat = bA->m[idxm][jdxm];
1122: return(0);
1123: }
1125: /*@
1126: MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
1128: Not collective
1130: Input Parameters:
1131: + A - nest matrix
1132: . idxm - index of the matrix within the nest matrix
1133: - jdxm - index of the matrix within the nest matrix
1135: Output Parameter:
1136: . sub - matrix at index idxm,jdxm within the nest matrix
1138: Level: developer
1140: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1141: MatNestGetLocalISs(), MatNestGetISs()
1142: @*/
1143: PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1144: {
1148: PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1149: return(0);
1150: }
1152: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1153: {
1154: Mat_Nest *bA = (Mat_Nest*)A->data;
1155: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1159: if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1160: if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1161: MatGetLocalSize(mat,&m,&n);
1162: MatGetSize(mat,&M,&N);
1163: ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1164: ISGetSize(bA->isglobal.row[idxm],&Mi);
1165: ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1166: ISGetSize(bA->isglobal.col[jdxm],&Ni);
1167: 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);
1168: 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);
1170: /* do not increase object state */
1171: if (mat == bA->m[idxm][jdxm]) return(0);
1173: PetscObjectReference((PetscObject)mat);
1174: MatDestroy(&bA->m[idxm][jdxm]);
1175: bA->m[idxm][jdxm] = mat;
1176: PetscObjectStateIncrease((PetscObject)A);
1177: MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1178: A->nonzerostate++;
1179: return(0);
1180: }
1182: /*@
1183: MatNestSetSubMat - Set a single submatrix in the nest matrix.
1185: Logically collective on the submatrix communicator
1187: Input Parameters:
1188: + A - nest matrix
1189: . idxm - index of the matrix within the nest matrix
1190: . jdxm - index of the matrix within the nest matrix
1191: - sub - matrix at index idxm,jdxm within the nest matrix
1193: Notes:
1194: The new submatrix must have the same size and communicator as that block of the nest.
1196: This increments the reference count of the submatrix.
1198: Level: developer
1200: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1201: MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1202: @*/
1203: PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1204: {
1208: PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1209: return(0);
1210: }
1212: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1213: {
1214: Mat_Nest *bA = (Mat_Nest*)A->data;
1217: if (M) *M = bA->nr;
1218: if (N) *N = bA->nc;
1219: if (mat) *mat = bA->m;
1220: return(0);
1221: }
1223: /*@C
1224: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
1226: Not collective
1228: Input Parameters:
1229: . A - nest matrix
1231: Output Parameter:
1232: + M - number of rows in the nest matrix
1233: . N - number of cols in the nest matrix
1234: - mat - 2d array of matrices
1236: Notes:
1238: The user should not free the array mat.
1240: In Fortran, this routine has a calling sequence
1241: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1242: where the space allocated for the optional argument mat is assumed large enough (if provided).
1244: Level: developer
1246: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1247: MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1248: @*/
1249: PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1250: {
1254: PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1255: return(0);
1256: }
1258: PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1259: {
1260: Mat_Nest *bA = (Mat_Nest*)A->data;
1263: if (M) *M = bA->nr;
1264: if (N) *N = bA->nc;
1265: return(0);
1266: }
1268: /*@
1269: MatNestGetSize - Returns the size of the nest matrix.
1271: Not collective
1273: Input Parameters:
1274: . A - nest matrix
1276: Output Parameter:
1277: + M - number of rows in the nested mat
1278: - N - number of cols in the nested mat
1280: Notes:
1282: Level: developer
1284: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1285: MatNestGetISs()
1286: @*/
1287: PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1288: {
1292: PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1293: return(0);
1294: }
1296: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1297: {
1298: Mat_Nest *vs = (Mat_Nest*)A->data;
1299: PetscInt i;
1302: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1303: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1304: return(0);
1305: }
1307: /*@C
1308: MatNestGetISs - Returns the index sets partitioning the row and column spaces
1310: Not collective
1312: Input Parameters:
1313: . A - nest matrix
1315: Output Parameter:
1316: + rows - array of row index sets
1317: - cols - array of column index sets
1319: Level: advanced
1321: Notes:
1322: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1324: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1325: MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1326: @*/
1327: PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1328: {
1333: PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1334: return(0);
1335: }
1337: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1338: {
1339: Mat_Nest *vs = (Mat_Nest*)A->data;
1340: PetscInt i;
1343: if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1344: if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1345: return(0);
1346: }
1348: /*@C
1349: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1351: Not collective
1353: Input Parameters:
1354: . A - nest matrix
1356: Output Parameter:
1357: + rows - array of row index sets (or NULL to ignore)
1358: - cols - array of column index sets (or NULL to ignore)
1360: Level: advanced
1362: Notes:
1363: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1365: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1366: MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1367: @*/
1368: PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1369: {
1374: PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1375: return(0);
1376: }
1378: PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1379: {
1381: PetscBool flg;
1384: PetscStrcmp(vtype,VECNEST,&flg);
1385: /* In reality, this only distinguishes VECNEST and "other" */
1386: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1387: else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1388: return(0);
1389: }
1391: /*@C
1392: MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1394: Not collective
1396: Input Parameters:
1397: + A - nest matrix
1398: - vtype - type to use for creating vectors
1400: Notes:
1402: Level: developer
1404: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1405: @*/
1406: PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1407: {
1411: PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1412: return(0);
1413: }
1415: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1416: {
1417: Mat_Nest *s = (Mat_Nest*)A->data;
1418: PetscInt i,j,m,n,M,N;
1420: PetscBool cong;
1423: MatReset_Nest(A);
1425: s->nr = nr;
1426: s->nc = nc;
1428: /* Create space for submatrices */
1429: PetscMalloc1(nr,&s->m);
1430: for (i=0; i<nr; i++) {
1431: PetscMalloc1(nc,&s->m[i]);
1432: }
1433: for (i=0; i<nr; i++) {
1434: for (j=0; j<nc; j++) {
1435: s->m[i][j] = a[i*nc+j];
1436: if (a[i*nc+j]) {
1437: PetscObjectReference((PetscObject)a[i*nc+j]);
1438: }
1439: }
1440: }
1442: MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);
1444: PetscMalloc1(nr,&s->row_len);
1445: PetscMalloc1(nc,&s->col_len);
1446: for (i=0; i<nr; i++) s->row_len[i]=-1;
1447: for (j=0; j<nc; j++) s->col_len[j]=-1;
1449: PetscCalloc1(nr*nc,&s->nnzstate);
1450: for (i=0; i<nr; i++) {
1451: for (j=0; j<nc; j++) {
1452: if (s->m[i][j]) {
1453: MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1454: }
1455: }
1456: }
1458: MatNestGetSizes_Private(A,&m,&n,&M,&N);
1460: PetscLayoutSetSize(A->rmap,M);
1461: PetscLayoutSetLocalSize(A->rmap,m);
1462: PetscLayoutSetSize(A->cmap,N);
1463: PetscLayoutSetLocalSize(A->cmap,n);
1465: PetscLayoutSetUp(A->rmap);
1466: PetscLayoutSetUp(A->cmap);
1468: /* disable operations that are not supported for non-square matrices,
1469: or matrices for which is_row != is_col */
1470: MatHasCongruentLayouts(A,&cong);
1471: if (cong && nr != nc) cong = PETSC_FALSE;
1472: if (cong) {
1473: for (i = 0; cong && i < nr; i++) {
1474: ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1475: }
1476: }
1477: if (!cong) {
1478: A->ops->missingdiagonal = NULL;
1479: A->ops->getdiagonal = NULL;
1480: A->ops->shift = NULL;
1481: A->ops->diagonalset = NULL;
1482: }
1484: PetscCalloc2(nr,&s->left,nc,&s->right);
1485: PetscObjectStateIncrease((PetscObject)A);
1486: A->nonzerostate++;
1487: return(0);
1488: }
1490: /*@
1491: MatNestSetSubMats - Sets the nested submatrices
1493: Collective on Mat
1495: Input Parameter:
1496: + A - nested matrix
1497: . nr - number of nested row blocks
1498: . is_row - index sets for each nested row block, or NULL to make contiguous
1499: . nc - number of nested column blocks
1500: . is_col - index sets for each nested column block, or NULL to make contiguous
1501: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1503: Notes: this always resets any submatrix information previously set
1505: Level: advanced
1507: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1508: @*/
1509: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1510: {
1512: PetscInt i;
1516: if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1517: if (nr && is_row) {
1520: }
1521: if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1522: if (nc && is_col) {
1525: }
1527: PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1528: return(0);
1529: }
1531: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1532: {
1534: PetscBool flg;
1535: PetscInt i,j,m,mi,*ix;
1538: *ltog = NULL;
1539: for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1540: if (islocal[i]) {
1541: ISGetLocalSize(islocal[i],&mi);
1542: flg = PETSC_TRUE; /* We found a non-trivial entry */
1543: } else {
1544: ISGetLocalSize(isglobal[i],&mi);
1545: }
1546: m += mi;
1547: }
1548: if (!flg) return(0);
1550: PetscMalloc1(m,&ix);
1551: for (i=0,m=0; i<n; i++) {
1552: ISLocalToGlobalMapping smap = NULL;
1553: Mat sub = NULL;
1554: PetscSF sf;
1555: PetscLayout map;
1556: const PetscInt *ix2;
1558: if (!colflg) {
1559: MatNestFindNonzeroSubMatRow(A,i,&sub);
1560: } else {
1561: MatNestFindNonzeroSubMatCol(A,i,&sub);
1562: }
1563: if (sub) {
1564: if (!colflg) {
1565: MatGetLocalToGlobalMapping(sub,&smap,NULL);
1566: } else {
1567: MatGetLocalToGlobalMapping(sub,NULL,&smap);
1568: }
1569: }
1570: /*
1571: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1572: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1573: */
1574: ISGetIndices(isglobal[i],&ix2);
1575: if (islocal[i]) {
1576: PetscInt *ilocal,*iremote;
1577: PetscInt mil,nleaves;
1579: ISGetLocalSize(islocal[i],&mi);
1580: if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1581: for (j=0; j<mi; j++) ix[m+j] = j;
1582: ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);
1584: /* PetscSFSetGraphLayout does not like negative indices */
1585: PetscMalloc2(mi,&ilocal,mi,&iremote);
1586: for (j=0, nleaves = 0; j<mi; j++) {
1587: if (ix[m+j] < 0) continue;
1588: ilocal[nleaves] = j;
1589: iremote[nleaves] = ix[m+j];
1590: nleaves++;
1591: }
1592: ISGetLocalSize(isglobal[i],&mil);
1593: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1594: PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1595: PetscLayoutSetLocalSize(map,mil);
1596: PetscLayoutSetUp(map);
1597: PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1598: PetscLayoutDestroy(&map);
1599: PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1600: PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1601: PetscSFDestroy(&sf);
1602: PetscFree2(ilocal,iremote);
1603: } else {
1604: ISGetLocalSize(isglobal[i],&mi);
1605: for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1606: }
1607: ISRestoreIndices(isglobal[i],&ix2);
1608: m += mi;
1609: }
1610: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1611: return(0);
1612: }
1615: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1616: /*
1617: nprocessors = NP
1618: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1619: proc 0: => (g_0,h_0,)
1620: proc 1: => (g_1,h_1,)
1621: ...
1622: proc nprocs-1: => (g_NP-1,h_NP-1,)
1624: proc 0: proc 1: proc nprocs-1:
1625: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1627: proc 0:
1628: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1629: proc 1:
1630: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1632: proc NP-1:
1633: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1634: */
1635: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1636: {
1637: Mat_Nest *vs = (Mat_Nest*)A->data;
1638: PetscInt i,j,offset,n,nsum,bs;
1640: Mat sub = NULL;
1643: PetscMalloc1(nr,&vs->isglobal.row);
1644: PetscMalloc1(nc,&vs->isglobal.col);
1645: if (is_row) { /* valid IS is passed in */
1646: /* refs on is[] are incremeneted */
1647: for (i=0; i<vs->nr; i++) {
1648: PetscObjectReference((PetscObject)is_row[i]);
1650: vs->isglobal.row[i] = is_row[i];
1651: }
1652: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1653: nsum = 0;
1654: for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1655: MatNestFindNonzeroSubMatRow(A,i,&sub);
1656: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1657: MatGetLocalSize(sub,&n,NULL);
1658: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1659: nsum += n;
1660: }
1661: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1662: offset -= nsum;
1663: for (i=0; i<vs->nr; i++) {
1664: MatNestFindNonzeroSubMatRow(A,i,&sub);
1665: MatGetLocalSize(sub,&n,NULL);
1666: MatGetBlockSizes(sub,&bs,NULL);
1667: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1668: ISSetBlockSize(vs->isglobal.row[i],bs);
1669: offset += n;
1670: }
1671: }
1673: if (is_col) { /* valid IS is passed in */
1674: /* refs on is[] are incremeneted */
1675: for (j=0; j<vs->nc; j++) {
1676: PetscObjectReference((PetscObject)is_col[j]);
1678: vs->isglobal.col[j] = is_col[j];
1679: }
1680: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1681: offset = A->cmap->rstart;
1682: nsum = 0;
1683: for (j=0; j<vs->nc; j++) {
1684: MatNestFindNonzeroSubMatCol(A,j,&sub);
1685: if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1686: MatGetLocalSize(sub,NULL,&n);
1687: if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1688: nsum += n;
1689: }
1690: MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1691: offset -= nsum;
1692: for (j=0; j<vs->nc; j++) {
1693: MatNestFindNonzeroSubMatCol(A,j,&sub);
1694: MatGetLocalSize(sub,NULL,&n);
1695: MatGetBlockSizes(sub,NULL,&bs);
1696: ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1697: ISSetBlockSize(vs->isglobal.col[j],bs);
1698: offset += n;
1699: }
1700: }
1702: /* Set up the local ISs */
1703: PetscMalloc1(vs->nr,&vs->islocal.row);
1704: PetscMalloc1(vs->nc,&vs->islocal.col);
1705: for (i=0,offset=0; i<vs->nr; i++) {
1706: IS isloc;
1707: ISLocalToGlobalMapping rmap = NULL;
1708: PetscInt nlocal,bs;
1709: MatNestFindNonzeroSubMatRow(A,i,&sub);
1710: if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1711: if (rmap) {
1712: MatGetBlockSizes(sub,&bs,NULL);
1713: ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1714: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1715: ISSetBlockSize(isloc,bs);
1716: } else {
1717: nlocal = 0;
1718: isloc = NULL;
1719: }
1720: vs->islocal.row[i] = isloc;
1721: offset += nlocal;
1722: }
1723: for (i=0,offset=0; i<vs->nc; i++) {
1724: IS isloc;
1725: ISLocalToGlobalMapping cmap = NULL;
1726: PetscInt nlocal,bs;
1727: MatNestFindNonzeroSubMatCol(A,i,&sub);
1728: if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1729: if (cmap) {
1730: MatGetBlockSizes(sub,NULL,&bs);
1731: ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1732: ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1733: ISSetBlockSize(isloc,bs);
1734: } else {
1735: nlocal = 0;
1736: isloc = NULL;
1737: }
1738: vs->islocal.col[i] = isloc;
1739: offset += nlocal;
1740: }
1742: /* Set up the aggregate ISLocalToGlobalMapping */
1743: {
1744: ISLocalToGlobalMapping rmap,cmap;
1745: MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1746: MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1747: if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1748: ISLocalToGlobalMappingDestroy(&rmap);
1749: ISLocalToGlobalMappingDestroy(&cmap);
1750: }
1752: if (PetscDefined(USE_DEBUG)) {
1753: for (i=0; i<vs->nr; i++) {
1754: for (j=0; j<vs->nc; j++) {
1755: PetscInt m,n,M,N,mi,ni,Mi,Ni;
1756: Mat B = vs->m[i][j];
1757: if (!B) continue;
1758: MatGetSize(B,&M,&N);
1759: MatGetLocalSize(B,&m,&n);
1760: ISGetSize(vs->isglobal.row[i],&Mi);
1761: ISGetSize(vs->isglobal.col[j],&Ni);
1762: ISGetLocalSize(vs->isglobal.row[i],&mi);
1763: ISGetLocalSize(vs->isglobal.col[j],&ni);
1764: 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);
1765: 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);
1766: }
1767: }
1768: }
1770: /* Set A->assembled if all non-null blocks are currently assembled */
1771: for (i=0; i<vs->nr; i++) {
1772: for (j=0; j<vs->nc; j++) {
1773: if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1774: }
1775: }
1776: A->assembled = PETSC_TRUE;
1777: return(0);
1778: }
1780: /*@C
1781: MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1783: Collective on Mat
1785: Input Parameter:
1786: + comm - Communicator for the new Mat
1787: . nr - number of nested row blocks
1788: . is_row - index sets for each nested row block, or NULL to make contiguous
1789: . nc - number of nested column blocks
1790: . is_col - index sets for each nested column block, or NULL to make contiguous
1791: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1793: Output Parameter:
1794: . B - new matrix
1796: Level: advanced
1798: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1799: MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1800: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1801: @*/
1802: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1803: {
1804: Mat A;
1808: *B = NULL;
1809: MatCreate(comm,&A);
1810: MatSetType(A,MATNEST);
1811: A->preallocated = PETSC_TRUE;
1812: MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1813: *B = A;
1814: return(0);
1815: }
1817: PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1818: {
1819: Mat_Nest *nest = (Mat_Nest*)A->data;
1820: Mat *trans;
1821: PetscScalar **avv;
1822: PetscScalar *vv;
1823: PetscInt **aii,**ajj;
1824: PetscInt *ii,*jj,*ci;
1825: PetscInt nr,nc,nnz,i,j;
1826: PetscBool done;
1830: MatGetSize(A,&nr,&nc);
1831: if (reuse == MAT_REUSE_MATRIX) {
1832: PetscInt rnr;
1834: MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1835: if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1836: if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1837: MatSeqAIJGetArray(*newmat,&vv);
1838: }
1839: /* extract CSR for nested SeqAIJ matrices */
1840: nnz = 0;
1841: PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1842: for (i=0; i<nest->nr; ++i) {
1843: for (j=0; j<nest->nc; ++j) {
1844: Mat B = nest->m[i][j];
1845: if (B) {
1846: PetscScalar *naa;
1847: PetscInt *nii,*njj,nnr;
1848: PetscBool istrans;
1850: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1851: if (istrans) {
1852: Mat Bt;
1854: MatTransposeGetMat(B,&Bt);
1855: MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1856: B = trans[i*nest->nc+j];
1857: }
1858: MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1859: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1860: MatSeqAIJGetArray(B,&naa);
1861: nnz += nii[nnr];
1863: aii[i*nest->nc+j] = nii;
1864: ajj[i*nest->nc+j] = njj;
1865: avv[i*nest->nc+j] = naa;
1866: }
1867: }
1868: }
1869: if (reuse != MAT_REUSE_MATRIX) {
1870: PetscMalloc1(nr+1,&ii);
1871: PetscMalloc1(nnz,&jj);
1872: PetscMalloc1(nnz,&vv);
1873: } else {
1874: if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1875: }
1877: /* new row pointer */
1878: PetscArrayzero(ii,nr+1);
1879: for (i=0; i<nest->nr; ++i) {
1880: PetscInt ncr,rst;
1882: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1883: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1884: for (j=0; j<nest->nc; ++j) {
1885: if (aii[i*nest->nc+j]) {
1886: PetscInt *nii = aii[i*nest->nc+j];
1887: PetscInt ir;
1889: for (ir=rst; ir<ncr+rst; ++ir) {
1890: ii[ir+1] += nii[1]-nii[0];
1891: nii++;
1892: }
1893: }
1894: }
1895: }
1896: for (i=0; i<nr; i++) ii[i+1] += ii[i];
1898: /* construct CSR for the new matrix */
1899: PetscCalloc1(nr,&ci);
1900: for (i=0; i<nest->nr; ++i) {
1901: PetscInt ncr,rst;
1903: ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1904: ISGetLocalSize(nest->isglobal.row[i],&ncr);
1905: for (j=0; j<nest->nc; ++j) {
1906: if (aii[i*nest->nc+j]) {
1907: PetscScalar *nvv = avv[i*nest->nc+j];
1908: PetscInt *nii = aii[i*nest->nc+j];
1909: PetscInt *njj = ajj[i*nest->nc+j];
1910: PetscInt ir,cst;
1912: ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1913: for (ir=rst; ir<ncr+rst; ++ir) {
1914: PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1916: for (ij=0;ij<rsize;ij++) {
1917: jj[ist+ij] = *njj+cst;
1918: vv[ist+ij] = *nvv;
1919: njj++;
1920: nvv++;
1921: }
1922: ci[ir] += rsize;
1923: nii++;
1924: }
1925: }
1926: }
1927: }
1928: PetscFree(ci);
1930: /* restore info */
1931: for (i=0; i<nest->nr; ++i) {
1932: for (j=0; j<nest->nc; ++j) {
1933: Mat B = nest->m[i][j];
1934: if (B) {
1935: PetscInt nnr = 0, k = i*nest->nc+j;
1937: B = (trans[k] ? trans[k] : B);
1938: MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1939: if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1940: MatSeqAIJRestoreArray(B,&avv[k]);
1941: MatDestroy(&trans[k]);
1942: }
1943: }
1944: }
1945: PetscFree4(aii,ajj,avv,trans);
1947: /* finalize newmat */
1948: if (reuse == MAT_INITIAL_MATRIX) {
1949: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1950: } else if (reuse == MAT_INPLACE_MATRIX) {
1951: Mat B;
1953: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1954: MatHeaderReplace(A,&B);
1955: }
1956: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1957: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1958: {
1959: Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1960: a->free_a = PETSC_TRUE;
1961: a->free_ij = PETSC_TRUE;
1962: }
1963: return(0);
1964: }
1966: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y,PetscScalar a,Mat X)
1967: {
1968: Mat_Nest *nest = (Mat_Nest*)X->data;
1969: PetscInt i,j,k,rstart;
1970: PetscBool flg;
1974: /* Fill by row */
1975: for (j=0; j<nest->nc; ++j) {
1976: /* Using global column indices and ISAllGather() is not scalable. */
1977: IS bNis;
1978: PetscInt bN;
1979: const PetscInt *bNindices;
1980: ISAllGather(nest->isglobal.col[j], &bNis);
1981: ISGetSize(bNis,&bN);
1982: ISGetIndices(bNis,&bNindices);
1983: for (i=0; i<nest->nr; ++i) {
1984: Mat B,D=NULL;
1985: PetscInt bm, br;
1986: const PetscInt *bmindices;
1987: B = nest->m[i][j];
1988: if (!B) continue;
1989: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&flg);
1990: if (flg) {
1991: PetscTryMethod(B,"MatTransposeGetMat_C",(Mat,Mat*),(B,&D));
1992: PetscTryMethod(B,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(B,&D));
1993: MatConvert(B,((PetscObject)D)->type_name,MAT_INITIAL_MATRIX,&D);
1994: B = D;
1995: }
1996: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQSBAIJ,MATMPISBAIJ,"");
1997: if (flg) {
1998: if (D) {
1999: MatConvert(D,MATBAIJ,MAT_INPLACE_MATRIX,&D);
2000: } else {
2001: MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&D);
2002: }
2003: B = D;
2004: }
2005: ISGetLocalSize(nest->isglobal.row[i],&bm);
2006: ISGetIndices(nest->isglobal.row[i],&bmindices);
2007: MatGetOwnershipRange(B,&rstart,NULL);
2008: for (br = 0; br < bm; ++br) {
2009: PetscInt row = bmindices[br], brncols, *cols;
2010: const PetscInt *brcols;
2011: const PetscScalar *brcoldata;
2012: PetscScalar *vals = NULL;
2013: MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2014: PetscMalloc1(brncols,&cols);
2015: for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2016: /*
2017: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2018: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2019: */
2020: if (a != 1.0) {
2021: PetscMalloc1(brncols,&vals);
2022: for (k=0; k<brncols; k++) vals[k] = a * brcoldata[k];
2023: MatSetValues(Y,1,&row,brncols,cols,vals,ADD_VALUES);
2024: PetscFree(vals);
2025: } else {
2026: MatSetValues(Y,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2027: }
2028: MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2029: PetscFree(cols);
2030: }
2031: if (D) {
2032: MatDestroy(&D);
2033: }
2034: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2035: }
2036: ISRestoreIndices(bNis,&bNindices);
2037: ISDestroy(&bNis);
2038: }
2039: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
2040: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
2041: return(0);
2042: }
2044: PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2045: {
2047: Mat_Nest *nest = (Mat_Nest*)A->data;
2048: PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart,cstart,cend;
2049: PetscMPIInt size;
2050: Mat C;
2053: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
2054: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
2055: PetscInt nf;
2056: PetscBool fast;
2058: PetscStrcmp(newtype,MATAIJ,&fast);
2059: if (!fast) {
2060: PetscStrcmp(newtype,MATSEQAIJ,&fast);
2061: }
2062: for (i=0; i<nest->nr && fast; ++i) {
2063: for (j=0; j<nest->nc && fast; ++j) {
2064: Mat B = nest->m[i][j];
2065: if (B) {
2066: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
2067: if (!fast) {
2068: PetscBool istrans;
2070: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
2071: if (istrans) {
2072: Mat Bt;
2074: MatTransposeGetMat(B,&Bt);
2075: PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
2076: }
2077: }
2078: }
2079: }
2080: }
2081: for (i=0, nf=0; i<nest->nr && fast; ++i) {
2082: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
2083: if (fast) {
2084: PetscInt f,s;
2086: ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
2087: if (f != nf || s != 1) { fast = PETSC_FALSE; }
2088: else {
2089: ISGetSize(nest->isglobal.row[i],&f);
2090: nf += f;
2091: }
2092: }
2093: }
2094: for (i=0, nf=0; i<nest->nc && fast; ++i) {
2095: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
2096: if (fast) {
2097: PetscInt f,s;
2099: ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
2100: if (f != nf || s != 1) { fast = PETSC_FALSE; }
2101: else {
2102: ISGetSize(nest->isglobal.col[i],&f);
2103: nf += f;
2104: }
2105: }
2106: }
2107: if (fast) {
2108: MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
2109: return(0);
2110: }
2111: }
2112: MatGetSize(A,&M,&N);
2113: MatGetLocalSize(A,&m,&n);
2114: MatGetOwnershipRangeColumn(A,&cstart,&cend);
2115: if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2116: else {
2117: MatCreate(PetscObjectComm((PetscObject)A),&C);
2118: MatSetType(C,newtype);
2119: MatSetSizes(C,m,n,M,N);
2120: }
2121: PetscMalloc1(2*m,&dnnz);
2122: onnz = dnnz + m;
2123: for (k=0; k<m; k++) {
2124: dnnz[k] = 0;
2125: onnz[k] = 0;
2126: }
2127: for (j=0; j<nest->nc; ++j) {
2128: IS bNis;
2129: PetscInt bN;
2130: const PetscInt *bNindices;
2131: /* Using global column indices and ISAllGather() is not scalable. */
2132: ISAllGather(nest->isglobal.col[j], &bNis);
2133: ISGetSize(bNis, &bN);
2134: ISGetIndices(bNis,&bNindices);
2135: for (i=0; i<nest->nr; ++i) {
2136: PetscSF bmsf;
2137: PetscSFNode *iremote;
2138: Mat B;
2139: PetscInt bm, *sub_dnnz,*sub_onnz, br;
2140: const PetscInt *bmindices;
2141: B = nest->m[i][j];
2142: if (!B) continue;
2143: ISGetLocalSize(nest->isglobal.row[i],&bm);
2144: ISGetIndices(nest->isglobal.row[i],&bmindices);
2145: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2146: PetscMalloc1(bm,&iremote);
2147: PetscMalloc1(bm,&sub_dnnz);
2148: PetscMalloc1(bm,&sub_onnz);
2149: for (k = 0; k < bm; ++k){
2150: sub_dnnz[k] = 0;
2151: sub_onnz[k] = 0;
2152: }
2153: /*
2154: Locate the owners for all of the locally-owned global row indices for this row block.
2155: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2156: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2157: */
2158: MatGetOwnershipRange(B,&rstart,NULL);
2159: for (br = 0; br < bm; ++br) {
2160: PetscInt row = bmindices[br], brncols, col;
2161: const PetscInt *brcols;
2162: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2163: PetscMPIInt rowowner = 0;
2164: PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2165: /* how many roots */
2166: iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2167: /* get nonzero pattern */
2168: MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2169: for (k=0; k<brncols; k++) {
2170: col = bNindices[brcols[k]];
2171: if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2172: sub_dnnz[br]++;
2173: } else {
2174: sub_onnz[br]++;
2175: }
2176: }
2177: MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2178: }
2179: ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2180: /* bsf will have to take care of disposing of bedges. */
2181: PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2182: PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2183: PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2184: PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2185: PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2186: PetscFree(sub_dnnz);
2187: PetscFree(sub_onnz);
2188: PetscSFDestroy(&bmsf);
2189: }
2190: ISRestoreIndices(bNis,&bNindices);
2191: ISDestroy(&bNis);
2192: }
2193: /* Resize preallocation if overestimated */
2194: for (i=0;i<m;i++) {
2195: dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2196: onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2197: }
2198: MatSeqAIJSetPreallocation(C,0,dnnz);
2199: MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2200: PetscFree(dnnz);
2201: MatAXPY_Dense_Nest(C,1.0,A);
2202: if (reuse == MAT_INPLACE_MATRIX) {
2203: MatHeaderReplace(A,&C);
2204: } else *newmat = C;
2205: return(0);
2206: }
2208: PetscErrorCode MatConvert_Nest_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2209: {
2210: Mat B;
2211: PetscInt m,n,M,N;
2215: MatGetSize(A,&M,&N);
2216: MatGetLocalSize(A,&m,&n);
2217: if (reuse == MAT_REUSE_MATRIX) {
2218: B = *newmat;
2219: MatZeroEntries(B);
2220: } else {
2221: MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,&B);
2222: }
2223: MatAXPY_Dense_Nest(B,1.0,A);
2224: if (reuse == MAT_INPLACE_MATRIX) {
2225: MatHeaderReplace(A,&B);
2226: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2227: return(0);
2228: }
2230: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2231: {
2232: Mat_Nest *bA = (Mat_Nest*)mat->data;
2233: MatOperation opAdd;
2234: PetscInt i,j,nr = bA->nr,nc = bA->nc;
2235: PetscBool flg;
2239: *has = PETSC_FALSE;
2240: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2241: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2242: for (j=0; j<nc; j++) {
2243: for (i=0; i<nr; i++) {
2244: if (!bA->m[i][j]) continue;
2245: MatHasOperation(bA->m[i][j],opAdd,&flg);
2246: if (!flg) return(0);
2247: }
2248: }
2249: }
2250: if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2251: return(0);
2252: }
2254: /*MC
2255: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2257: Level: intermediate
2259: Notes:
2260: This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2261: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2262: It is usually used with DMComposite and DMCreateMatrix()
2264: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2265: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2266: than the nest matrix.
2268: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2269: VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2270: MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2271: M*/
2272: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2273: {
2274: Mat_Nest *s;
2278: PetscNewLog(A,&s);
2279: A->data = (void*)s;
2281: s->nr = -1;
2282: s->nc = -1;
2283: s->m = NULL;
2284: s->splitassembly = PETSC_FALSE;
2286: PetscMemzero(A->ops,sizeof(*A->ops));
2288: A->ops->mult = MatMult_Nest;
2289: A->ops->multadd = MatMultAdd_Nest;
2290: A->ops->multtranspose = MatMultTranspose_Nest;
2291: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2292: A->ops->transpose = MatTranspose_Nest;
2293: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2294: A->ops->assemblyend = MatAssemblyEnd_Nest;
2295: A->ops->zeroentries = MatZeroEntries_Nest;
2296: A->ops->copy = MatCopy_Nest;
2297: A->ops->axpy = MatAXPY_Nest;
2298: A->ops->duplicate = MatDuplicate_Nest;
2299: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2300: A->ops->destroy = MatDestroy_Nest;
2301: A->ops->view = MatView_Nest;
2302: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2303: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2304: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2305: A->ops->getdiagonal = MatGetDiagonal_Nest;
2306: A->ops->diagonalscale = MatDiagonalScale_Nest;
2307: A->ops->scale = MatScale_Nest;
2308: A->ops->shift = MatShift_Nest;
2309: A->ops->diagonalset = MatDiagonalSet_Nest;
2310: A->ops->setrandom = MatSetRandom_Nest;
2311: A->ops->hasoperation = MatHasOperation_Nest;
2312: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2314: A->spptr = NULL;
2315: A->assembled = PETSC_FALSE;
2317: /* expose Nest api's */
2318: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2319: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2320: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2321: PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);
2322: PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);
2323: PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2324: PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);
2325: PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2326: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2327: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2328: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2329: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);
2330: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",MatConvert_Nest_Dense);
2331: PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",MatConvert_Nest_Dense);
2332: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2333: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2334: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);
2336: PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2337: return(0);
2338: }