Actual source code: matnest.c

petsc-master 2019-11-13
Report Typos and Errors

  2:  #include <../src/mat/impls/nest/matnestimpl.h>
  3:  #include <../src/mat/impls/aij/seq/aij.h>
  4:  #include <petscsf.h>

  6: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
  7: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left);
  8: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*);

 10: /* private functions */
 11: static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
 12: {
 13:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 14:   PetscInt       i,j;

 18:   *m = *n = *M = *N = 0;
 19:   for (i=0; i<bA->nr; i++) {  /* rows */
 20:     PetscInt sm,sM;
 21:     ISGetLocalSize(bA->isglobal.row[i],&sm);
 22:     ISGetSize(bA->isglobal.row[i],&sM);
 23:     *m  += sm;
 24:     *M  += sM;
 25:   }
 26:   for (j=0; j<bA->nc; j++) {  /* cols */
 27:     PetscInt sn,sN;
 28:     ISGetLocalSize(bA->isglobal.col[j],&sn);
 29:     ISGetSize(bA->isglobal.col[j],&sN);
 30:     *n  += sn;
 31:     *N  += sN;
 32:   }
 33:   return(0);
 34: }

 36: /* operations */
 37: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
 38: {
 39:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 40:   Vec            *bx = bA->right,*by = bA->left;
 41:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

 45:   for (i=0; i<nr; i++) {VecGetSubVector(y,bA->isglobal.row[i],&by[i]);}
 46:   for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
 47:   for (i=0; i<nr; i++) {
 48:     VecZeroEntries(by[i]);
 49:     for (j=0; j<nc; j++) {
 50:       if (!bA->m[i][j]) continue;
 51:       /* y[i] <- y[i] + A[i][j] * x[j] */
 52:       MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);
 53:     }
 54:   }
 55:   for (i=0; i<nr; i++) {VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);}
 56:   for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
 57:   return(0);
 58: }

 60: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
 61: {
 62:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 63:   Vec            *bx = bA->right,*bz = bA->left;
 64:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

 68:   for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
 69:   for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
 70:   for (i=0; i<nr; i++) {
 71:     if (y != z) {
 72:       Vec by;
 73:       VecGetSubVector(y,bA->isglobal.row[i],&by);
 74:       VecCopy(by,bz[i]);
 75:       VecRestoreSubVector(y,bA->isglobal.row[i],&by);
 76:     }
 77:     for (j=0; j<nc; j++) {
 78:       if (!bA->m[i][j]) continue;
 79:       /* y[i] <- y[i] + A[i][j] * x[j] */
 80:       MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
 81:     }
 82:   }
 83:   for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
 84:   for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
 85:   return(0);
 86: }

 88: typedef struct {
 89:   Mat          *workC;    /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
 90:   PetscScalar  *tarray;   /* buffer for storing all temporary products A[i][j] B[j] */
 91:   PetscInt     *dm,*dn,k; /* displacements and number of submatrices */
 92: } Nest_Dense;

 94: PETSC_INTERN PetscErrorCode MatMatMultNumeric_Nest_Dense(Mat A,Mat B,Mat C)
 95: {
 96:   Mat_Nest          *bA = (Mat_Nest*)A->data;
 97:   PetscContainer    container;
 98:   Nest_Dense        *contents;
 99:   Mat               viewB,viewC,seq;
100:   const PetscScalar *barray;
101:   PetscScalar       *carray;
102:   PetscInt          i,j,M,N,nr = bA->nr,nc = bA->nc,ldb,ldc;
103:   PetscErrorCode    ierr;

106:   PetscObjectQuery((PetscObject)C,"workC",(PetscObject*)&container);
107:   if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exist");
108:   PetscContainerGetPointer(container,(void**)&contents);
109:   MatDenseGetLDA(B,&ldb);
110:   MatDenseGetLDA(C,&ldc);
111:   MatGetSize(B,NULL,&N);
112:   MatZeroEntries(C);
113:   MatDenseGetArrayRead(B,&barray);
114:   MatDenseGetArray(C,&carray);
115:   for (i=0; i<nr; i++) {
116:     ISGetSize(bA->isglobal.row[i],&M);
117:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
118:     MatDenseGetLocalMatrix(viewC,&seq);
119:     MatSeqDenseSetLDA(seq,ldc);
120:     for (j=0; j<nc; j++) {
121:       if (!bA->m[i][j]) continue;
122:       ISGetSize(bA->isglobal.col[j],&M);
123:       MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
124:       MatDenseGetLocalMatrix(viewB,&seq);
125:       MatSeqDenseSetLDA(seq,ldb);
126:       /* workC <- A[i][j] * B[j] */
127:       MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]);
128:       MatDestroy(&viewB);
129:        /* C[i] <- workC + C[i] */
130:       MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
131:     }
132:     MatDestroy(&viewC);
133:   }
134:   MatDenseRestoreArray(C,&carray);
135:   MatDenseRestoreArrayRead(B,&barray);
136:   return(0);
137: }

139: PetscErrorCode MatNest_DenseDestroy(void *ctx)
140: {
141:   Nest_Dense     *contents = (Nest_Dense*)ctx;
142:   PetscInt       i;

146:   PetscFree(contents->tarray);
147:   for (i=0; i<contents->k; i++) {
148:     MatDestroy(contents->workC + i);
149:   }
150:   PetscFree3(contents->dm,contents->dn,contents->workC);
151:   PetscFree(contents);
152:   return(0);
153: }

155: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_Nest_Dense(Mat A,Mat B,PetscReal fill,Mat *C)
156: {
157:   Mat_Nest          *bA = (Mat_Nest*)A->data;
158:   Mat               viewB,viewSeq;
159:   const PetscScalar *barray;
160:   PetscInt          i,j,M,N,m,nr = bA->nr,nc = bA->nc,maxm = 0,ldb;
161:   PetscContainer    container;
162:   Nest_Dense        *contents;
163:   PetscErrorCode    ierr;

166:   MatGetSize(B,NULL,&N);
167:   if (!(*C)) {
168:     MatGetLocalSize(A,&m,NULL);
169:     MatGetSize(A,&M,NULL);
170:     MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,C);
171:   } else {
172:     if ((*C)->rmap->n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrix local row dimensions are incompatible, %D != %D",(*C)->rmap->n,A->rmap->n);
173:     if ((*C)->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrix local column dimensions are incompatible, %D != %D",(*C)->cmap->n,B->cmap->n);
174:   }

176:   PetscNew(&contents);
177:   PetscContainerCreate(PetscObjectComm((PetscObject)A),&container);
178:   PetscContainerSetPointer(container,contents);
179:   PetscContainerSetUserDestroy(container,MatNest_DenseDestroy);
180:   PetscObjectCompose((PetscObject)*C,"workC",(PetscObject)container);
181:   PetscContainerDestroy(&container);
182:   PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
183:   contents->k = nr*nc;
184:   for (i=0; i<nr; i++) {
185:     ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
186:     maxm = PetscMax(maxm,contents->dm[i+1]);
187:     contents->dm[i+1] += contents->dm[i];
188:   }
189:   for (i=0; i<nc; i++) {
190:     ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
191:     contents->dn[i+1] += contents->dn[i];
192:   }
193:   PetscMalloc1(maxm*N,&contents->tarray);
194:   MatDenseGetLDA(B,&ldb);
195:   MatGetSize(B,NULL,&N);
196:   MatDenseGetArrayRead(B,&barray);
197:   /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
198:   for (j=0; j<nc; j++) {
199:     ISGetSize(bA->isglobal.col[j],&M);
200:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
201:     MatDenseGetLocalMatrix(viewB,&viewSeq);
202:     MatSeqDenseSetLDA(viewSeq,ldb);
203:     for (i=0; i<nr; i++) {
204:       if (!bA->m[i][j]) continue;
205:       /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
206:       MatMatMultSymbolic(bA->m[i][j],viewB,fill,contents->workC + i*nc + j);
207:       MatDenseGetLocalMatrix(contents->workC[i*nc + j],&viewSeq);
208:       /* free the memory allocated in MatMatMultSymbolic, since tarray will be shared by all Mat */
209:       MatSeqDenseSetPreallocation(viewSeq,contents->tarray);
210:     }
211:     MatDestroy(&viewB);
212:   }
213:   MatDenseRestoreArrayRead(B,&barray);

215:   (*C)->ops->matmultnumeric = MatMatMultNumeric_Nest_Dense;
216:   return(0);
217: }

219: PETSC_INTERN PetscErrorCode MatMatMult_Nest_Dense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
220: {

224:   if (scall == MAT_INITIAL_MATRIX) {
225:     *C = NULL;
226:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
227:     MatMatMultSymbolic_Nest_Dense(A,B,fill,C);
228:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
229:   }
230:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
231:   MatMatMultNumeric_Nest_Dense(A,B,*C);
232:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
233:   return(0);
234: }

236: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
237: {
238:   Mat_Nest       *bA = (Mat_Nest*)A->data;
239:   Vec            *bx = bA->left,*by = bA->right;
240:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

244:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
245:   for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
246:   for (j=0; j<nc; j++) {
247:     VecZeroEntries(by[j]);
248:     for (i=0; i<nr; i++) {
249:       if (!bA->m[i][j]) continue;
250:       /* y[j] <- y[j] + (A[i][j])^T * x[i] */
251:       MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
252:     }
253:   }
254:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
255:   for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
256:   return(0);
257: }

259: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
260: {
261:   Mat_Nest       *bA = (Mat_Nest*)A->data;
262:   Vec            *bx = bA->left,*bz = bA->right;
263:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

267:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
268:   for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
269:   for (j=0; j<nc; j++) {
270:     if (y != z) {
271:       Vec by;
272:       VecGetSubVector(y,bA->isglobal.col[j],&by);
273:       VecCopy(by,bz[j]);
274:       VecRestoreSubVector(y,bA->isglobal.col[j],&by);
275:     }
276:     for (i=0; i<nr; i++) {
277:       if (!bA->m[i][j]) continue;
278:       /* z[j] <- y[j] + (A[i][j])^T * x[i] */
279:       MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
280:     }
281:   }
282:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
283:   for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
284:   return(0);
285: }

287: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
288: {
289:   Mat_Nest       *bA = (Mat_Nest*)A->data, *bC;
290:   Mat            C;
291:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

295:   if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");

297:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
298:     Mat *subs;
299:     IS  *is_row,*is_col;

301:     PetscCalloc1(nr * nc,&subs);
302:     PetscMalloc2(nr,&is_row,nc,&is_col);
303:     MatNestGetISs(A,is_row,is_col);
304:     if (reuse == MAT_INPLACE_MATRIX) {
305:       for (i=0; i<nr; i++) {
306:         for (j=0; j<nc; j++) {
307:           subs[i + nr * j] = bA->m[i][j];
308:         }
309:       }
310:     }

312:     MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
313:     PetscFree(subs);
314:     PetscFree2(is_row,is_col);
315:   } else {
316:     C = *B;
317:   }

319:   bC = (Mat_Nest*)C->data;
320:   for (i=0; i<nr; i++) {
321:     for (j=0; j<nc; j++) {
322:       if (bA->m[i][j]) {
323:         MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
324:       } else {
325:         bC->m[j][i] = NULL;
326:       }
327:     }
328:   }

330:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
331:     *B = C;
332:   } else {
333:     MatHeaderMerge(A, &C);
334:   }
335:   return(0);
336: }

338: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
339: {
341:   IS             *lst = *list;
342:   PetscInt       i;

345:   if (!lst) return(0);
346:   for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
347:   PetscFree(lst);
348:   *list = NULL;
349:   return(0);
350: }

352: static PetscErrorCode MatDestroy_Nest(Mat A)
353: {
354:   Mat_Nest       *vs = (Mat_Nest*)A->data;
355:   PetscInt       i,j;

359:   /* release the matrices and the place holders */
360:   MatNestDestroyISList(vs->nr,&vs->isglobal.row);
361:   MatNestDestroyISList(vs->nc,&vs->isglobal.col);
362:   MatNestDestroyISList(vs->nr,&vs->islocal.row);
363:   MatNestDestroyISList(vs->nc,&vs->islocal.col);

365:   PetscFree(vs->row_len);
366:   PetscFree(vs->col_len);

368:   PetscFree2(vs->left,vs->right);

370:   /* release the matrices and the place holders */
371:   if (vs->m) {
372:     for (i=0; i<vs->nr; i++) {
373:       for (j=0; j<vs->nc; j++) {
374:         MatDestroy(&vs->m[i][j]);
375:       }
376:       PetscFree(vs->m[i]);
377:     }
378:     PetscFree(vs->m);
379:   }
380:   PetscFree(A->data);

382:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
383:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
384:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
385:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
386:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
387:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
388:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
389:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
390:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
391:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
392:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
393:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
394:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_mpidense_C",0);
395:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_seqdense_C",0);
396:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_dense_C",0);
397:   return(0);
398: }

400: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
401: {
402:   Mat_Nest       *vs = (Mat_Nest*)A->data;
403:   PetscInt       i,j;

407:   for (i=0; i<vs->nr; i++) {
408:     for (j=0; j<vs->nc; j++) {
409:       if (vs->m[i][j]) {
410:         MatAssemblyBegin(vs->m[i][j],type);
411:         if (!vs->splitassembly) {
412:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
413:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
414:            * already performing an assembly, but the result would by more complicated and appears to offer less
415:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
416:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
417:            */
418:           MatAssemblyEnd(vs->m[i][j],type);
419:         }
420:       }
421:     }
422:   }
423:   return(0);
424: }

426: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
427: {
428:   Mat_Nest       *vs = (Mat_Nest*)A->data;
429:   PetscInt       i,j;

433:   for (i=0; i<vs->nr; i++) {
434:     for (j=0; j<vs->nc; j++) {
435:       if (vs->m[i][j]) {
436:         if (vs->splitassembly) {
437:           MatAssemblyEnd(vs->m[i][j],type);
438:         }
439:       }
440:     }
441:   }
442:   return(0);
443: }

445: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
446: {
448:   Mat_Nest       *vs = (Mat_Nest*)A->data;
449:   PetscInt       j;
450:   Mat            sub;

453:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
454:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
455:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
456:   *B = sub;
457:   return(0);
458: }

460: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
461: {
463:   Mat_Nest       *vs = (Mat_Nest*)A->data;
464:   PetscInt       i;
465:   Mat            sub;

468:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
469:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
470:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
471:   *B = sub;
472:   return(0);
473: }

475: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
476: {
478:   PetscInt       i;
479:   PetscBool      flg;

485:   *found = -1;
486:   for (i=0; i<n; i++) {
487:     if (!list[i]) continue;
488:     ISEqual(list[i],is,&flg);
489:     if (flg) {
490:       *found = i;
491:       return(0);
492:     }
493:   }
494:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
495:   return(0);
496: }

498: /* Get a block row as a new MatNest */
499: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
500: {
501:   Mat_Nest       *vs = (Mat_Nest*)A->data;
502:   char           keyname[256];

506:   *B   = NULL;
507:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
508:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
509:   if (*B) return(0);

511:   MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);

513:   (*B)->assembled = A->assembled;

515:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
516:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
517:   return(0);
518: }

520: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
521: {
522:   Mat_Nest       *vs = (Mat_Nest*)A->data;
524:   PetscInt       row,col;
525:   PetscBool      same,isFullCol,isFullColGlobal;

528:   /* Check if full column space. This is a hack */
529:   isFullCol = PETSC_FALSE;
530:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
531:   if (same) {
532:     PetscInt n,first,step,i,an,am,afirst,astep;
533:     ISStrideGetInfo(iscol,&first,&step);
534:     ISGetLocalSize(iscol,&n);
535:     isFullCol = PETSC_TRUE;
536:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
537:       ISStrideGetInfo(is->col[i],&afirst,&astep);
538:       ISGetLocalSize(is->col[i],&am);
539:       if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
540:       an += am;
541:     }
542:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
543:   }
544:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

546:   if (isFullColGlobal && vs->nc > 1) {
547:     PetscInt row;
548:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
549:     MatNestGetRow(A,row,B);
550:   } else {
551:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
552:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
553:     if (!vs->m[row][col]) {
554:       PetscInt lr,lc;

556:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
557:       ISGetLocalSize(vs->isglobal.row[row],&lr);
558:       ISGetLocalSize(vs->isglobal.col[col],&lc);
559:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
560:       MatSetUp(vs->m[row][col]);
561:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
562:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
563:     }
564:     *B = vs->m[row][col];
565:   }
566:   return(0);
567: }

569: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
570: {
572:   Mat_Nest       *vs = (Mat_Nest*)A->data;
573:   Mat            sub;

576:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
577:   switch (reuse) {
578:   case MAT_INITIAL_MATRIX:
579:     if (sub) { PetscObjectReference((PetscObject)sub); }
580:     *B = sub;
581:     break;
582:   case MAT_REUSE_MATRIX:
583:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
584:     break;
585:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
586:     break;
587:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
588:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
589:     break;
590:   }
591:   return(0);
592: }

594: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
595: {
597:   Mat_Nest       *vs = (Mat_Nest*)A->data;
598:   Mat            sub;

601:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
602:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
603:   if (sub) {PetscObjectReference((PetscObject)sub);}
604:   *B = sub;
605:   return(0);
606: }

608: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
609: {
611:   Mat_Nest       *vs = (Mat_Nest*)A->data;
612:   Mat            sub;

615:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
616:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
617:   if (sub) {
618:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
619:     MatDestroy(B);
620:   }
621:   return(0);
622: }

624: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
625: {
626:   Mat_Nest       *bA = (Mat_Nest*)A->data;
627:   PetscInt       i;

631:   for (i=0; i<bA->nr; i++) {
632:     Vec bv;
633:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
634:     if (bA->m[i][i]) {
635:       MatGetDiagonal(bA->m[i][i],bv);
636:     } else {
637:       VecSet(bv,0.0);
638:     }
639:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
640:   }
641:   return(0);
642: }

644: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
645: {
646:   Mat_Nest       *bA = (Mat_Nest*)A->data;
647:   Vec            bl,*br;
648:   PetscInt       i,j;

652:   PetscCalloc1(bA->nc,&br);
653:   if (r) {
654:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
655:   }
656:   bl = NULL;
657:   for (i=0; i<bA->nr; i++) {
658:     if (l) {
659:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
660:     }
661:     for (j=0; j<bA->nc; j++) {
662:       if (bA->m[i][j]) {
663:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
664:       }
665:     }
666:     if (l) {
667:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
668:     }
669:   }
670:   if (r) {
671:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
672:   }
673:   PetscFree(br);
674:   return(0);
675: }

677: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
678: {
679:   Mat_Nest       *bA = (Mat_Nest*)A->data;
680:   PetscInt       i,j;

684:   for (i=0; i<bA->nr; i++) {
685:     for (j=0; j<bA->nc; j++) {
686:       if (bA->m[i][j]) {
687:         MatScale(bA->m[i][j],a);
688:       }
689:     }
690:   }
691:   return(0);
692: }

694: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
695: {
696:   Mat_Nest       *bA = (Mat_Nest*)A->data;
697:   PetscInt       i;

701:   for (i=0; i<bA->nr; i++) {
702:     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);
703:     MatShift(bA->m[i][i],a);
704:   }
705:   return(0);
706: }

708: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
709: {
710:   Mat_Nest       *bA = (Mat_Nest*)A->data;
711:   PetscInt       i;

715:   for (i=0; i<bA->nr; i++) {
716:     Vec bv;
717:     VecGetSubVector(D,bA->isglobal.row[i],&bv);
718:     if (bA->m[i][i]) {
719:       MatDiagonalSet(bA->m[i][i],bv,is);
720:     }
721:     VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
722:   }
723:   return(0);
724: }

726: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
727: {
728:   Mat_Nest       *bA = (Mat_Nest*)A->data;
729:   PetscInt       i,j;

733:   for (i=0; i<bA->nr; i++) {
734:     for (j=0; j<bA->nc; j++) {
735:       if (bA->m[i][j]) {
736:         MatSetRandom(bA->m[i][j],rctx);
737:       }
738:     }
739:   }
740:   return(0);
741: }

743: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
744: {
745:   Mat_Nest       *bA = (Mat_Nest*)A->data;
746:   Vec            *L,*R;
747:   MPI_Comm       comm;
748:   PetscInt       i,j;

752:   PetscObjectGetComm((PetscObject)A,&comm);
753:   if (right) {
754:     /* allocate R */
755:     PetscMalloc1(bA->nc, &R);
756:     /* Create the right vectors */
757:     for (j=0; j<bA->nc; j++) {
758:       for (i=0; i<bA->nr; i++) {
759:         if (bA->m[i][j]) {
760:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
761:           break;
762:         }
763:       }
764:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
765:     }
766:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
767:     /* hand back control to the nest vector */
768:     for (j=0; j<bA->nc; j++) {
769:       VecDestroy(&R[j]);
770:     }
771:     PetscFree(R);
772:   }

774:   if (left) {
775:     /* allocate L */
776:     PetscMalloc1(bA->nr, &L);
777:     /* Create the left vectors */
778:     for (i=0; i<bA->nr; i++) {
779:       for (j=0; j<bA->nc; j++) {
780:         if (bA->m[i][j]) {
781:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
782:           break;
783:         }
784:       }
785:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
786:     }

788:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
789:     for (i=0; i<bA->nr; i++) {
790:       VecDestroy(&L[i]);
791:     }

793:     PetscFree(L);
794:   }
795:   return(0);
796: }

798: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
799: {
800:   Mat_Nest       *bA = (Mat_Nest*)A->data;
801:   PetscBool      isascii,viewSub = PETSC_FALSE;
802:   PetscInt       i,j;

806:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
807:   if (isascii) {

809:     PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
810:     PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
811:     PetscViewerASCIIPushTab(viewer);
812:     PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);

814:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
815:     for (i=0; i<bA->nr; i++) {
816:       for (j=0; j<bA->nc; j++) {
817:         MatType   type;
818:         char      name[256] = "",prefix[256] = "";
819:         PetscInt  NR,NC;
820:         PetscBool isNest = PETSC_FALSE;

822:         if (!bA->m[i][j]) {
823:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
824:           continue;
825:         }
826:         MatGetSize(bA->m[i][j],&NR,&NC);
827:         MatGetType(bA->m[i][j], &type);
828:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
829:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
830:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

832:         PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);

834:         if (isNest || viewSub) {
835:           PetscViewerASCIIPushTab(viewer);  /* push1 */
836:           MatView(bA->m[i][j],viewer);
837:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
838:         }
839:       }
840:     }
841:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
842:   }
843:   return(0);
844: }

846: static PetscErrorCode MatZeroEntries_Nest(Mat A)
847: {
848:   Mat_Nest       *bA = (Mat_Nest*)A->data;
849:   PetscInt       i,j;

853:   for (i=0; i<bA->nr; i++) {
854:     for (j=0; j<bA->nc; j++) {
855:       if (!bA->m[i][j]) continue;
856:       MatZeroEntries(bA->m[i][j]);
857:     }
858:   }
859:   return(0);
860: }

862: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
863: {
864:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
865:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

869:   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);
870:   for (i=0; i<nr; i++) {
871:     for (j=0; j<nc; j++) {
872:       if (bA->m[i][j] && bB->m[i][j]) {
873:         MatCopy(bA->m[i][j],bB->m[i][j],str);
874:       } 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);
875:     }
876:   }
877:   PetscObjectStateIncrease((PetscObject)B);
878:   return(0);
879: }

881: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
882: {
883:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
884:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;

888:   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);
889:   for (i=0; i<nr; i++) {
890:     for (j=0; j<nc; j++) {
891:       if (bY->m[i][j] && bX->m[i][j]) {
892:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
893:       } else if (bX->m[i][j]) {
894:         Mat M;

896:         if (str != DIFFERENT_NONZERO_PATTERN) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j);
897:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
898:         MatNestSetSubMat(Y,i,j,M);
899:         MatDestroy(&M);
900:       }
901:     }
902:   }
903:   return(0);
904: }

906: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
907: {
908:   Mat_Nest       *bA = (Mat_Nest*)A->data;
909:   Mat            *b;
910:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

914:   PetscMalloc1(nr*nc,&b);
915:   for (i=0; i<nr; i++) {
916:     for (j=0; j<nc; j++) {
917:       if (bA->m[i][j]) {
918:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
919:       } else {
920:         b[i*nc+j] = NULL;
921:       }
922:     }
923:   }
924:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
925:   /* Give the new MatNest exclusive ownership */
926:   for (i=0; i<nr*nc; i++) {
927:     MatDestroy(&b[i]);
928:   }
929:   PetscFree(b);

931:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
932:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
933:   return(0);
934: }

936: /* nest api */
937: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
938: {
939:   Mat_Nest *bA = (Mat_Nest*)A->data;

942:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
943:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
944:   *mat = bA->m[idxm][jdxm];
945:   return(0);
946: }

948: /*@
949:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

951:  Not collective

953:  Input Parameters:
954: +   A  - nest matrix
955: .   idxm - index of the matrix within the nest matrix
956: -   jdxm - index of the matrix within the nest matrix

958:  Output Parameter:
959: .   sub - matrix at index idxm,jdxm within the nest matrix

961:  Level: developer

963: .seealso: MatNestGetSize(), MatNestGetSubMats()
964: @*/
965: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
966: {

970:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
971:   return(0);
972: }

974: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
975: {
976:   Mat_Nest       *bA = (Mat_Nest*)A->data;
977:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

981:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
982:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
983:   MatGetLocalSize(mat,&m,&n);
984:   MatGetSize(mat,&M,&N);
985:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
986:   ISGetSize(bA->isglobal.row[idxm],&Mi);
987:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
988:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
989:   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);
990:   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);

992:   PetscObjectReference((PetscObject)mat);
993:   MatDestroy(&bA->m[idxm][jdxm]);
994:   bA->m[idxm][jdxm] = mat;
995:   return(0);
996: }

998: /*@
999:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

1001:  Logically collective on the submatrix communicator

1003:  Input Parameters:
1004: +   A  - nest matrix
1005: .   idxm - index of the matrix within the nest matrix
1006: .   jdxm - index of the matrix within the nest matrix
1007: -   sub - matrix at index idxm,jdxm within the nest matrix

1009:  Notes:
1010:  The new submatrix must have the same size and communicator as that block of the nest.

1012:  This increments the reference count of the submatrix.

1014:  Level: developer

1016: .seealso: MatNestSetSubMats(), MatNestGetSubMats()
1017: @*/
1018: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1019: {

1023:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1024:   return(0);
1025: }

1027: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1028: {
1029:   Mat_Nest *bA = (Mat_Nest*)A->data;

1032:   if (M)   *M   = bA->nr;
1033:   if (N)   *N   = bA->nc;
1034:   if (mat) *mat = bA->m;
1035:   return(0);
1036: }

1038: /*@C
1039:  MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.

1041:  Not collective

1043:  Input Parameters:
1044: .   A  - nest matrix

1046:  Output Parameter:
1047: +   M - number of rows in the nest matrix
1048: .   N - number of cols in the nest matrix
1049: -   mat - 2d array of matrices

1051:  Notes:

1053:  The user should not free the array mat.

1055:  In Fortran, this routine has a calling sequence
1056: $   call MatNestGetSubMats(A, M, N, mat, ierr)
1057:  where the space allocated for the optional argument mat is assumed large enough (if provided).

1059:  Level: developer

1061: .seealso: MatNestGetSize(), MatNestGetSubMat()
1062: @*/
1063: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1064: {

1068:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1069:   return(0);
1070: }

1072: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1073: {
1074:   Mat_Nest *bA = (Mat_Nest*)A->data;

1077:   if (M) *M = bA->nr;
1078:   if (N) *N = bA->nc;
1079:   return(0);
1080: }

1082: /*@
1083:  MatNestGetSize - Returns the size of the nest matrix.

1085:  Not collective

1087:  Input Parameters:
1088: .   A  - nest matrix

1090:  Output Parameter:
1091: +   M - number of rows in the nested mat
1092: -   N - number of cols in the nested mat

1094:  Notes:

1096:  Level: developer

1098: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
1099: @*/
1100: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1101: {

1105:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1106:   return(0);
1107: }

1109: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1110: {
1111:   Mat_Nest *vs = (Mat_Nest*)A->data;
1112:   PetscInt i;

1115:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1116:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1117:   return(0);
1118: }

1120: /*@C
1121:  MatNestGetISs - Returns the index sets partitioning the row and column spaces

1123:  Not collective

1125:  Input Parameters:
1126: .   A  - nest matrix

1128:  Output Parameter:
1129: +   rows - array of row index sets
1130: -   cols - array of column index sets

1132:  Level: advanced

1134:  Notes:
1135:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.

1137: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
1138: @*/
1139: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
1140: {

1145:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1146:   return(0);
1147: }

1149: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1150: {
1151:   Mat_Nest *vs = (Mat_Nest*)A->data;
1152:   PetscInt i;

1155:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1156:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1157:   return(0);
1158: }

1160: /*@C
1161:  MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces

1163:  Not collective

1165:  Input Parameters:
1166: .   A  - nest matrix

1168:  Output Parameter:
1169: +   rows - array of row index sets (or NULL to ignore)
1170: -   cols - array of column index sets (or NULL to ignore)

1172:  Level: advanced

1174:  Notes:
1175:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.

1177: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
1178: @*/
1179: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1180: {

1185:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1186:   return(0);
1187: }

1189: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1190: {
1192:   PetscBool      flg;

1195:   PetscStrcmp(vtype,VECNEST,&flg);
1196:   /* In reality, this only distinguishes VECNEST and "other" */
1197:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1198:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1199:   return(0);
1200: }

1202: /*@C
1203:  MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()

1205:  Not collective

1207:  Input Parameters:
1208: +  A  - nest matrix
1209: -  vtype - type to use for creating vectors

1211:  Notes:

1213:  Level: developer

1215: .seealso: MatCreateVecs()
1216: @*/
1217: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1218: {

1222:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1223:   return(0);
1224: }

1226: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1227: {
1228:   Mat_Nest       *s = (Mat_Nest*)A->data;
1229:   PetscInt       i,j,m,n,M,N;

1233:   s->nr = nr;
1234:   s->nc = nc;

1236:   /* Create space for submatrices */
1237:   PetscMalloc1(nr,&s->m);
1238:   for (i=0; i<nr; i++) {
1239:     PetscMalloc1(nc,&s->m[i]);
1240:   }
1241:   for (i=0; i<nr; i++) {
1242:     for (j=0; j<nc; j++) {
1243:       s->m[i][j] = a[i*nc+j];
1244:       if (a[i*nc+j]) {
1245:         PetscObjectReference((PetscObject)a[i*nc+j]);
1246:       }
1247:     }
1248:   }

1250:   MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);

1252:   PetscMalloc1(nr,&s->row_len);
1253:   PetscMalloc1(nc,&s->col_len);
1254:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1255:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1257:   MatNestGetSizes_Private(A,&m,&n,&M,&N);

1259:   PetscLayoutSetSize(A->rmap,M);
1260:   PetscLayoutSetLocalSize(A->rmap,m);
1261:   PetscLayoutSetSize(A->cmap,N);
1262:   PetscLayoutSetLocalSize(A->cmap,n);

1264:   PetscLayoutSetUp(A->rmap);
1265:   PetscLayoutSetUp(A->cmap);

1267:   PetscCalloc2(nr,&s->left,nc,&s->right);
1268:   return(0);
1269: }

1271: /*@
1272:    MatNestSetSubMats - Sets the nested submatrices

1274:    Collective on Mat

1276:    Input Parameter:
1277: +  A - nested matrix
1278: .  nr - number of nested row blocks
1279: .  is_row - index sets for each nested row block, or NULL to make contiguous
1280: .  nc - number of nested column blocks
1281: .  is_col - index sets for each nested column block, or NULL to make contiguous
1282: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1284:    Level: advanced

1286: .seealso: MatCreateNest(), MATNEST
1287: @*/
1288: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1289: {
1291:   PetscInt       i,nr_nc;

1295:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1296:   if (nr && is_row) {
1299:   }
1300:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1301:   if (nc && is_col) {
1304:   }
1305:   nr_nc=nr*nc;
1307:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1308:   return(0);
1309: }

1311: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1312: {
1314:   PetscBool      flg;
1315:   PetscInt       i,j,m,mi,*ix;

1318:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1319:     if (islocal[i]) {
1320:       ISGetSize(islocal[i],&mi);
1321:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1322:     } else {
1323:       ISGetSize(isglobal[i],&mi);
1324:     }
1325:     m += mi;
1326:   }
1327:   if (flg) {
1328:     PetscMalloc1(m,&ix);
1329:     for (i=0,m=0; i<n; i++) {
1330:       ISLocalToGlobalMapping smap = NULL;
1331:       Mat                    sub = NULL;
1332:       PetscSF                sf;
1333:       PetscLayout            map;
1334:       PetscInt               *ix2;

1336:       if (!colflg) {
1337:         MatNestFindNonzeroSubMatRow(A,i,&sub);
1338:       } else {
1339:         MatNestFindNonzeroSubMatCol(A,i,&sub);
1340:       }
1341:       if (sub) {
1342:         if (!colflg) {
1343:           MatGetLocalToGlobalMapping(sub,&smap,NULL);
1344:         } else {
1345:           MatGetLocalToGlobalMapping(sub,NULL,&smap);
1346:         }
1347:       }
1348:       if (islocal[i]) {
1349:         ISGetSize(islocal[i],&mi);
1350:       } else {
1351:         ISGetSize(isglobal[i],&mi);
1352:       }
1353:       for (j=0; j<mi; j++) ix[m+j] = j;
1354:       if (smap) {ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);}

1356:       /*
1357:         Now we need to extract the monolithic global indices that correspond to the given split global indices.
1358:         In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1359:        */
1360:       PetscMalloc1(mi,&ix2);
1361:       PetscSFCreate(((PetscObject)isglobal[i])->comm,&sf);
1362:       PetscLayoutCreate(((PetscObject)isglobal[i])->comm,&map);
1363:       PetscLayoutSetLocalSize(map,mi);
1364:       PetscLayoutSetUp(map);
1365:       PetscSFSetGraphLayout(sf,map,mi,NULL,PETSC_USE_POINTER,ix+m);
1366:       PetscLayoutDestroy(&map);
1367:       for (j=0; j<mi; j++) ix2[j] = ix[m+j];
1368:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1369:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1370:       PetscSFDestroy(&sf);
1371:       PetscFree(ix2);
1372:       m   += mi;
1373:     }
1374:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1375:   } else {
1376:     *ltog = NULL;
1377:   }
1378:   return(0);
1379: }


1382: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1383: /*
1384:   nprocessors = NP
1385:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1386:        proc 0: => (g_0,h_0,)
1387:        proc 1: => (g_1,h_1,)
1388:        ...
1389:        proc nprocs-1: => (g_NP-1,h_NP-1,)

1391:             proc 0:                      proc 1:                    proc nprocs-1:
1392:     is[0] = (0,1,2,...,nlocal(g_0)-1)  (0,1,...,nlocal(g_1)-1)  (0,1,...,nlocal(g_NP-1))

1394:             proc 0:
1395:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1396:             proc 1:
1397:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1399:             proc NP-1:
1400:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1401: */
1402: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1403: {
1404:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1405:   PetscInt       i,j,offset,n,nsum,bs;
1407:   Mat            sub = NULL;

1410:   PetscMalloc1(nr,&vs->isglobal.row);
1411:   PetscMalloc1(nc,&vs->isglobal.col);
1412:   if (is_row) { /* valid IS is passed in */
1413:     /* refs on is[] are incremeneted */
1414:     for (i=0; i<vs->nr; i++) {
1415:       PetscObjectReference((PetscObject)is_row[i]);

1417:       vs->isglobal.row[i] = is_row[i];
1418:     }
1419:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1420:     nsum = 0;
1421:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1422:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1423:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1424:       MatGetLocalSize(sub,&n,NULL);
1425:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1426:       nsum += n;
1427:     }
1428:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1429:     offset -= nsum;
1430:     for (i=0; i<vs->nr; i++) {
1431:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1432:       MatGetLocalSize(sub,&n,NULL);
1433:       MatGetBlockSize(sub,&bs);
1434:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1435:       ISSetBlockSize(vs->isglobal.row[i],bs);
1436:       offset += n;
1437:     }
1438:   }

1440:   if (is_col) { /* valid IS is passed in */
1441:     /* refs on is[] are incremeneted */
1442:     for (j=0; j<vs->nc; j++) {
1443:       PetscObjectReference((PetscObject)is_col[j]);

1445:       vs->isglobal.col[j] = is_col[j];
1446:     }
1447:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1448:     offset = A->cmap->rstart;
1449:     nsum   = 0;
1450:     for (j=0; j<vs->nc; j++) {
1451:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1452:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1453:       MatGetLocalSize(sub,NULL,&n);
1454:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1455:       nsum += n;
1456:     }
1457:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1458:     offset -= nsum;
1459:     for (j=0; j<vs->nc; j++) {
1460:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1461:       MatGetLocalSize(sub,NULL,&n);
1462:       MatGetBlockSize(sub,&bs);
1463:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1464:       ISSetBlockSize(vs->isglobal.col[j],bs);
1465:       offset += n;
1466:     }
1467:   }

1469:   /* Set up the local ISs */
1470:   PetscMalloc1(vs->nr,&vs->islocal.row);
1471:   PetscMalloc1(vs->nc,&vs->islocal.col);
1472:   for (i=0,offset=0; i<vs->nr; i++) {
1473:     IS                     isloc;
1474:     ISLocalToGlobalMapping rmap = NULL;
1475:     PetscInt               nlocal,bs;
1476:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1477:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1478:     if (rmap) {
1479:       MatGetBlockSize(sub,&bs);
1480:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1481:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1482:       ISSetBlockSize(isloc,bs);
1483:     } else {
1484:       nlocal = 0;
1485:       isloc  = NULL;
1486:     }
1487:     vs->islocal.row[i] = isloc;
1488:     offset            += nlocal;
1489:   }
1490:   for (i=0,offset=0; i<vs->nc; i++) {
1491:     IS                     isloc;
1492:     ISLocalToGlobalMapping cmap = NULL;
1493:     PetscInt               nlocal,bs;
1494:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1495:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1496:     if (cmap) {
1497:       MatGetBlockSize(sub,&bs);
1498:       ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1499:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1500:       ISSetBlockSize(isloc,bs);
1501:     } else {
1502:       nlocal = 0;
1503:       isloc  = NULL;
1504:     }
1505:     vs->islocal.col[i] = isloc;
1506:     offset            += nlocal;
1507:   }

1509:   /* Set up the aggregate ISLocalToGlobalMapping */
1510:   {
1511:     ISLocalToGlobalMapping rmap,cmap;
1512:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1513:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1514:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1515:     ISLocalToGlobalMappingDestroy(&rmap);
1516:     ISLocalToGlobalMappingDestroy(&cmap);
1517:   }

1519: #if defined(PETSC_USE_DEBUG)
1520:   for (i=0; i<vs->nr; i++) {
1521:     for (j=0; j<vs->nc; j++) {
1522:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1523:       Mat      B = vs->m[i][j];
1524:       if (!B) continue;
1525:       MatGetSize(B,&M,&N);
1526:       MatGetLocalSize(B,&m,&n);
1527:       ISGetSize(vs->isglobal.row[i],&Mi);
1528:       ISGetSize(vs->isglobal.col[j],&Ni);
1529:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1530:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1531:       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);
1532:       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);
1533:     }
1534:   }
1535: #endif

1537:   /* Set A->assembled if all non-null blocks are currently assembled */
1538:   for (i=0; i<vs->nr; i++) {
1539:     for (j=0; j<vs->nc; j++) {
1540:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1541:     }
1542:   }
1543:   A->assembled = PETSC_TRUE;
1544:   return(0);
1545: }

1547: /*@C
1548:    MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately

1550:    Collective on Mat

1552:    Input Parameter:
1553: +  comm - Communicator for the new Mat
1554: .  nr - number of nested row blocks
1555: .  is_row - index sets for each nested row block, or NULL to make contiguous
1556: .  nc - number of nested column blocks
1557: .  is_col - index sets for each nested column block, or NULL to make contiguous
1558: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1560:    Output Parameter:
1561: .  B - new matrix

1563:    Level: advanced

1565: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1566: @*/
1567: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1568: {
1569:   Mat            A;

1573:   *B   = 0;
1574:   MatCreate(comm,&A);
1575:   MatSetType(A,MATNEST);
1576:   A->preallocated = PETSC_TRUE;
1577:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1578:   *B   = A;
1579:   return(0);
1580: }

1582: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1583: {
1584:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1585:   Mat            *trans;
1586:   PetscScalar    **avv;
1587:   PetscScalar    *vv;
1588:   PetscInt       **aii,**ajj;
1589:   PetscInt       *ii,*jj,*ci;
1590:   PetscInt       nr,nc,nnz,i,j;
1591:   PetscBool      done;

1595:   MatGetSize(A,&nr,&nc);
1596:   if (reuse == MAT_REUSE_MATRIX) {
1597:     PetscInt rnr;

1599:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1600:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1601:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1602:     MatSeqAIJGetArray(*newmat,&vv);
1603:   }
1604:   /* extract CSR for nested SeqAIJ matrices */
1605:   nnz  = 0;
1606:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1607:   for (i=0; i<nest->nr; ++i) {
1608:     for (j=0; j<nest->nc; ++j) {
1609:       Mat B = nest->m[i][j];
1610:       if (B) {
1611:         PetscScalar *naa;
1612:         PetscInt    *nii,*njj,nnr;
1613:         PetscBool   istrans;

1615:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1616:         if (istrans) {
1617:           Mat Bt;

1619:           MatTransposeGetMat(B,&Bt);
1620:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1621:           B    = trans[i*nest->nc+j];
1622:         }
1623:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1624:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1625:         MatSeqAIJGetArray(B,&naa);
1626:         nnz += nii[nnr];

1628:         aii[i*nest->nc+j] = nii;
1629:         ajj[i*nest->nc+j] = njj;
1630:         avv[i*nest->nc+j] = naa;
1631:       }
1632:     }
1633:   }
1634:   if (reuse != MAT_REUSE_MATRIX) {
1635:     PetscMalloc1(nr+1,&ii);
1636:     PetscMalloc1(nnz,&jj);
1637:     PetscMalloc1(nnz,&vv);
1638:   } else {
1639:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1640:   }

1642:   /* new row pointer */
1643:   PetscArrayzero(ii,nr+1);
1644:   for (i=0; i<nest->nr; ++i) {
1645:     PetscInt       ncr,rst;

1647:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1648:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1649:     for (j=0; j<nest->nc; ++j) {
1650:       if (aii[i*nest->nc+j]) {
1651:         PetscInt    *nii = aii[i*nest->nc+j];
1652:         PetscInt    ir;

1654:         for (ir=rst; ir<ncr+rst; ++ir) {
1655:           ii[ir+1] += nii[1]-nii[0];
1656:           nii++;
1657:         }
1658:       }
1659:     }
1660:   }
1661:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1663:   /* construct CSR for the new matrix */
1664:   PetscCalloc1(nr,&ci);
1665:   for (i=0; i<nest->nr; ++i) {
1666:     PetscInt       ncr,rst;

1668:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1669:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1670:     for (j=0; j<nest->nc; ++j) {
1671:       if (aii[i*nest->nc+j]) {
1672:         PetscScalar *nvv = avv[i*nest->nc+j];
1673:         PetscInt    *nii = aii[i*nest->nc+j];
1674:         PetscInt    *njj = ajj[i*nest->nc+j];
1675:         PetscInt    ir,cst;

1677:         ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1678:         for (ir=rst; ir<ncr+rst; ++ir) {
1679:           PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];

1681:           for (ij=0;ij<rsize;ij++) {
1682:             jj[ist+ij] = *njj+cst;
1683:             vv[ist+ij] = *nvv;
1684:             njj++;
1685:             nvv++;
1686:           }
1687:           ci[ir] += rsize;
1688:           nii++;
1689:         }
1690:       }
1691:     }
1692:   }
1693:   PetscFree(ci);

1695:   /* restore info */
1696:   for (i=0; i<nest->nr; ++i) {
1697:     for (j=0; j<nest->nc; ++j) {
1698:       Mat B = nest->m[i][j];
1699:       if (B) {
1700:         PetscInt nnr = 0, k = i*nest->nc+j;

1702:         B    = (trans[k] ? trans[k] : B);
1703:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1704:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1705:         MatSeqAIJRestoreArray(B,&avv[k]);
1706:         MatDestroy(&trans[k]);
1707:       }
1708:     }
1709:   }
1710:   PetscFree4(aii,ajj,avv,trans);

1712:   /* finalize newmat */
1713:   if (reuse == MAT_INITIAL_MATRIX) {
1714:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1715:   } else if (reuse == MAT_INPLACE_MATRIX) {
1716:     Mat B;

1718:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1719:     MatHeaderReplace(A,&B);
1720:   }
1721:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1722:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1723:   {
1724:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1725:     a->free_a     = PETSC_TRUE;
1726:     a->free_ij    = PETSC_TRUE;
1727:   }
1728:   return(0);
1729: }

1731: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1732: {
1734:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1735:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1736:   PetscInt       cstart,cend;
1737:   PetscMPIInt    size;
1738:   Mat            C;

1741:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1742:   if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1743:     PetscInt  nf;
1744:     PetscBool fast;

1746:     PetscStrcmp(newtype,MATAIJ,&fast);
1747:     if (!fast) {
1748:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1749:     }
1750:     for (i=0; i<nest->nr && fast; ++i) {
1751:       for (j=0; j<nest->nc && fast; ++j) {
1752:         Mat B = nest->m[i][j];
1753:         if (B) {
1754:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1755:           if (!fast) {
1756:             PetscBool istrans;

1758:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1759:             if (istrans) {
1760:               Mat Bt;

1762:               MatTransposeGetMat(B,&Bt);
1763:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1764:             }
1765:           }
1766:         }
1767:       }
1768:     }
1769:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1770:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1771:       if (fast) {
1772:         PetscInt f,s;

1774:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1775:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1776:         else {
1777:           ISGetSize(nest->isglobal.row[i],&f);
1778:           nf  += f;
1779:         }
1780:       }
1781:     }
1782:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1783:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1784:       if (fast) {
1785:         PetscInt f,s;

1787:         ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1788:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1789:         else {
1790:           ISGetSize(nest->isglobal.col[i],&f);
1791:           nf  += f;
1792:         }
1793:       }
1794:     }
1795:     if (fast) {
1796:       MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1797:       return(0);
1798:     }
1799:   }
1800:   MatGetSize(A,&M,&N);
1801:   MatGetLocalSize(A,&m,&n);
1802:   MatGetOwnershipRangeColumn(A,&cstart,&cend);
1803:   switch (reuse) {
1804:   case MAT_INITIAL_MATRIX:
1805:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1806:     MatSetType(C,newtype);
1807:     MatSetSizes(C,m,n,M,N);
1808:     *newmat = C;
1809:     break;
1810:   case MAT_REUSE_MATRIX:
1811:     C = *newmat;
1812:     break;
1813:   default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1814:   }
1815:   PetscMalloc1(2*m,&dnnz);
1816:   onnz = dnnz + m;
1817:   for (k=0; k<m; k++) {
1818:     dnnz[k] = 0;
1819:     onnz[k] = 0;
1820:   }
1821:   for (j=0; j<nest->nc; ++j) {
1822:     IS             bNis;
1823:     PetscInt       bN;
1824:     const PetscInt *bNindices;
1825:     /* Using global column indices and ISAllGather() is not scalable. */
1826:     ISAllGather(nest->isglobal.col[j], &bNis);
1827:     ISGetSize(bNis, &bN);
1828:     ISGetIndices(bNis,&bNindices);
1829:     for (i=0; i<nest->nr; ++i) {
1830:       PetscSF        bmsf;
1831:       PetscSFNode    *iremote;
1832:       Mat            B;
1833:       PetscInt       bm, *sub_dnnz,*sub_onnz, br;
1834:       const PetscInt *bmindices;
1835:       B = nest->m[i][j];
1836:       if (!B) continue;
1837:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1838:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1839:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1840:       PetscMalloc1(bm,&iremote);
1841:       PetscMalloc1(bm,&sub_dnnz);
1842:       PetscMalloc1(bm,&sub_onnz);
1843:       for (k = 0; k < bm; ++k){
1844:             sub_dnnz[k] = 0;
1845:             sub_onnz[k] = 0;
1846:       }
1847:       /*
1848:        Locate the owners for all of the locally-owned global row indices for this row block.
1849:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
1850:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1851:        */
1852:       MatGetOwnershipRange(B,&rstart,NULL);
1853:       for (br = 0; br < bm; ++br) {
1854:         PetscInt       row = bmindices[br], brncols, col;
1855:         const PetscInt *brcols;
1856:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
1857:         PetscMPIInt    rowowner = 0;
1858:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1859:         /* how many roots  */
1860:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
1861:         /* get nonzero pattern */
1862:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
1863:         for (k=0; k<brncols; k++) {
1864:           col  = bNindices[brcols[k]];
1865:           if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
1866:             sub_dnnz[br]++;
1867:           } else {
1868:             sub_onnz[br]++;
1869:           }
1870:         }
1871:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
1872:       }
1873:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1874:       /* bsf will have to take care of disposing of bedges. */
1875:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1876:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1877:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1878:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1879:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1880:       PetscFree(sub_dnnz);
1881:       PetscFree(sub_onnz);
1882:       PetscSFDestroy(&bmsf);
1883:     }
1884:     ISRestoreIndices(bNis,&bNindices);
1885:     ISDestroy(&bNis);
1886:   }
1887:   /* Resize preallocation if overestimated */
1888:   for (i=0;i<m;i++) {
1889:     dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
1890:     onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
1891:   }
1892:   MatSeqAIJSetPreallocation(C,0,dnnz);
1893:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1894:   PetscFree(dnnz);

1896:   /* Fill by row */
1897:   for (j=0; j<nest->nc; ++j) {
1898:     /* Using global column indices and ISAllGather() is not scalable. */
1899:     IS             bNis;
1900:     PetscInt       bN;
1901:     const PetscInt *bNindices;
1902:     ISAllGather(nest->isglobal.col[j], &bNis);
1903:     ISGetSize(bNis,&bN);
1904:     ISGetIndices(bNis,&bNindices);
1905:     for (i=0; i<nest->nr; ++i) {
1906:       Mat            B;
1907:       PetscInt       bm, br;
1908:       const PetscInt *bmindices;
1909:       B = nest->m[i][j];
1910:       if (!B) continue;
1911:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1912:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1913:       MatGetOwnershipRange(B,&rstart,NULL);
1914:       for (br = 0; br < bm; ++br) {
1915:         PetscInt          row = bmindices[br], brncols,  *cols;
1916:         const PetscInt    *brcols;
1917:         const PetscScalar *brcoldata;
1918:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1919:         PetscMalloc1(brncols,&cols);
1920:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1921:         /*
1922:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1923:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1924:          */
1925:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1926:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1927:         PetscFree(cols);
1928:       }
1929:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1930:     }
1931:     ISRestoreIndices(bNis,&bNindices);
1932:     ISDestroy(&bNis);
1933:   }
1934:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1935:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1936:   return(0);
1937: }

1939: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
1940: {
1941:   Mat_Nest       *bA = (Mat_Nest*)mat->data;
1942:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
1943:   PetscBool      flg;

1947:   *has = PETSC_FALSE;
1948:   if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MAT_MULT) {
1949:     for (j=0; j<nc; j++) {
1950:       for (i=0; i<nr; i++) {
1951:         if (!bA->m[i][j]) continue;
1952:         MatHasOperation(bA->m[i][j],op,&flg);
1953:         if (!flg) return(0);
1954:       }
1955:     }
1956:   }
1957:   if (((void**)mat->ops)[op] || (op == MATOP_MAT_MULT && flg)) *has = PETSC_TRUE;
1958:   return(0);
1959: }

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

1964:   Level: intermediate

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

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

1975: .seealso: MatCreate(), MatType, MatCreateNest()
1976: M*/
1977: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1978: {
1979:   Mat_Nest       *s;

1983:   PetscNewLog(A,&s);
1984:   A->data = (void*)s;

1986:   s->nr            = -1;
1987:   s->nc            = -1;
1988:   s->m             = NULL;
1989:   s->splitassembly = PETSC_FALSE;

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

1993:   A->ops->mult                  = MatMult_Nest;
1994:   A->ops->multadd               = MatMultAdd_Nest;
1995:   A->ops->multtranspose         = MatMultTranspose_Nest;
1996:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
1997:   A->ops->transpose             = MatTranspose_Nest;
1998:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
1999:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2000:   A->ops->zeroentries           = MatZeroEntries_Nest;
2001:   A->ops->copy                  = MatCopy_Nest;
2002:   A->ops->axpy                  = MatAXPY_Nest;
2003:   A->ops->duplicate             = MatDuplicate_Nest;
2004:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2005:   A->ops->destroy               = MatDestroy_Nest;
2006:   A->ops->view                  = MatView_Nest;
2007:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2008:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2009:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2010:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2011:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2012:   A->ops->scale                 = MatScale_Nest;
2013:   A->ops->shift                 = MatShift_Nest;
2014:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2015:   A->ops->setrandom             = MatSetRandom_Nest;
2016:   A->ops->hasoperation          = MatHasOperation_Nest;

2018:   A->spptr        = 0;
2019:   A->assembled    = PETSC_FALSE;

2021:   /* expose Nest api's */
2022:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",        MatNestGetSubMat_Nest);
2023:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",        MatNestSetSubMat_Nest);
2024:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",       MatNestGetSubMats_Nest);
2025:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",          MatNestGetSize_Nest);
2026:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",           MatNestGetISs_Nest);
2027:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",      MatNestGetLocalISs_Nest);
2028:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",       MatNestSetVecType_Nest);
2029:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",       MatNestSetSubMats_Nest);
2030:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",  MatConvert_Nest_AIJ);
2031:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",  MatConvert_Nest_AIJ);
2032:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",     MatConvert_Nest_AIJ);
2033:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",      MatConvert_Nest_IS);
2034:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_mpidense_C",MatMatMult_Nest_Dense);
2035:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_seqdense_C",MatMatMult_Nest_Dense);
2036:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_dense_C",   MatMatMult_Nest_Dense);

2038:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2039:   return(0);
2040: }