Actual source code: matnest.c

petsc-main 2021-04-20
Report Typos and Errors
  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:   switch (reuse) {
2116:   case MAT_INITIAL_MATRIX:
2117:     MatCreate(PetscObjectComm((PetscObject)A),&C);
2118:     MatSetType(C,newtype);
2119:     MatSetSizes(C,m,n,M,N);
2120:     *newmat = C;
2121:     break;
2122:   case MAT_REUSE_MATRIX:
2123:     C = *newmat;
2124:     break;
2125:   default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
2126:   }
2127:   PetscMalloc1(2*m,&dnnz);
2128:   onnz = dnnz + m;
2129:   for (k=0; k<m; k++) {
2130:     dnnz[k] = 0;
2131:     onnz[k] = 0;
2132:   }
2133:   for (j=0; j<nest->nc; ++j) {
2134:     IS             bNis;
2135:     PetscInt       bN;
2136:     const PetscInt *bNindices;
2137:     /* Using global column indices and ISAllGather() is not scalable. */
2138:     ISAllGather(nest->isglobal.col[j], &bNis);
2139:     ISGetSize(bNis, &bN);
2140:     ISGetIndices(bNis,&bNindices);
2141:     for (i=0; i<nest->nr; ++i) {
2142:       PetscSF        bmsf;
2143:       PetscSFNode    *iremote;
2144:       Mat            B;
2145:       PetscInt       bm, *sub_dnnz,*sub_onnz, br;
2146:       const PetscInt *bmindices;
2147:       B = nest->m[i][j];
2148:       if (!B) continue;
2149:       ISGetLocalSize(nest->isglobal.row[i],&bm);
2150:       ISGetIndices(nest->isglobal.row[i],&bmindices);
2151:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2152:       PetscMalloc1(bm,&iremote);
2153:       PetscMalloc1(bm,&sub_dnnz);
2154:       PetscMalloc1(bm,&sub_onnz);
2155:       for (k = 0; k < bm; ++k){
2156:         sub_dnnz[k] = 0;
2157:         sub_onnz[k] = 0;
2158:       }
2159:       /*
2160:        Locate the owners for all of the locally-owned global row indices for this row block.
2161:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
2162:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2163:        */
2164:       MatGetOwnershipRange(B,&rstart,NULL);
2165:       for (br = 0; br < bm; ++br) {
2166:         PetscInt       row = bmindices[br], brncols, col;
2167:         const PetscInt *brcols;
2168:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
2169:         PetscMPIInt    rowowner = 0;
2170:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2171:         /* how many roots  */
2172:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
2173:         /* get nonzero pattern */
2174:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2175:         for (k=0; k<brncols; k++) {
2176:           col  = bNindices[brcols[k]];
2177:           if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2178:             sub_dnnz[br]++;
2179:           } else {
2180:             sub_onnz[br]++;
2181:           }
2182:         }
2183:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2184:       }
2185:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2186:       /* bsf will have to take care of disposing of bedges. */
2187:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2188:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2189:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2190:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2191:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2192:       PetscFree(sub_dnnz);
2193:       PetscFree(sub_onnz);
2194:       PetscSFDestroy(&bmsf);
2195:     }
2196:     ISRestoreIndices(bNis,&bNindices);
2197:     ISDestroy(&bNis);
2198:   }
2199:   /* Resize preallocation if overestimated */
2200:   for (i=0;i<m;i++) {
2201:     dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2202:     onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2203:   }
2204:   MatSeqAIJSetPreallocation(C,0,dnnz);
2205:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2206:   PetscFree(dnnz);
2207:   MatAXPY_Dense_Nest(C,1.0,A);
2208:   return(0);
2209: }

2211: PetscErrorCode MatConvert_Nest_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2212: {
2213:   Mat            B;
2214:   PetscInt       m,n,M,N;

2218:   MatGetSize(A,&M,&N);
2219:   MatGetLocalSize(A,&m,&n);
2220:   if (reuse == MAT_REUSE_MATRIX) {
2221:     B = *newmat;
2222:     MatZeroEntries(B);
2223:   } else {
2224:     MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,&B);
2225:   }
2226:   MatAXPY_Dense_Nest(B,1.0,A);
2227:   if (reuse == MAT_INPLACE_MATRIX) {
2228:     MatHeaderReplace(A,&B);
2229:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2230:   return(0);
2231: }

2233: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2234: {
2235:   Mat_Nest       *bA = (Mat_Nest*)mat->data;
2236:   MatOperation   opAdd;
2237:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
2238:   PetscBool      flg;

2242:   *has = PETSC_FALSE;
2243:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2244:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2245:     for (j=0; j<nc; j++) {
2246:       for (i=0; i<nr; i++) {
2247:         if (!bA->m[i][j]) continue;
2248:         MatHasOperation(bA->m[i][j],opAdd,&flg);
2249:         if (!flg) return(0);
2250:       }
2251:     }
2252:   }
2253:   if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2254:   return(0);
2255: }

2257: /*MC
2258:   MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.

2260:   Level: intermediate

2262:   Notes:
2263:   This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2264:   It allows the use of symmetric and block formats for parts of multi-physics simulations.
2265:   It is usually used with DMComposite and DMCreateMatrix()

2267:   Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2268:   rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2269:   than the nest matrix.

2271: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2272:           VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2273:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2274: M*/
2275: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2276: {
2277:   Mat_Nest       *s;

2281:   PetscNewLog(A,&s);
2282:   A->data = (void*)s;

2284:   s->nr            = -1;
2285:   s->nc            = -1;
2286:   s->m             = NULL;
2287:   s->splitassembly = PETSC_FALSE;

2289:   PetscMemzero(A->ops,sizeof(*A->ops));

2291:   A->ops->mult                  = MatMult_Nest;
2292:   A->ops->multadd               = MatMultAdd_Nest;
2293:   A->ops->multtranspose         = MatMultTranspose_Nest;
2294:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2295:   A->ops->transpose             = MatTranspose_Nest;
2296:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2297:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2298:   A->ops->zeroentries           = MatZeroEntries_Nest;
2299:   A->ops->copy                  = MatCopy_Nest;
2300:   A->ops->axpy                  = MatAXPY_Nest;
2301:   A->ops->duplicate             = MatDuplicate_Nest;
2302:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2303:   A->ops->destroy               = MatDestroy_Nest;
2304:   A->ops->view                  = MatView_Nest;
2305:   A->ops->getvecs               = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2306:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2307:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2308:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2309:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2310:   A->ops->scale                 = MatScale_Nest;
2311:   A->ops->shift                 = MatShift_Nest;
2312:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2313:   A->ops->setrandom             = MatSetRandom_Nest;
2314:   A->ops->hasoperation          = MatHasOperation_Nest;
2315:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2317:   A->spptr        = NULL;
2318:   A->assembled    = PETSC_FALSE;

2320:   /* expose Nest api's */
2321:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",        MatNestGetSubMat_Nest);
2322:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",        MatNestSetSubMat_Nest);
2323:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",       MatNestGetSubMats_Nest);
2324:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",          MatNestGetSize_Nest);
2325:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",           MatNestGetISs_Nest);
2326:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",      MatNestGetLocalISs_Nest);
2327:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",       MatNestSetVecType_Nest);
2328:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",       MatNestSetSubMats_Nest);
2329:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",  MatConvert_Nest_AIJ);
2330:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",  MatConvert_Nest_AIJ);
2331:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",     MatConvert_Nest_AIJ);
2332:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",      MatConvert_Nest_IS);
2333:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",MatConvert_Nest_Dense);
2334:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",MatConvert_Nest_Dense);
2335:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2336:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2337:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);

2339:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2340:   return(0);
2341: }