Actual source code: submat.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/matimpl.h>          /*I "petscmat.h" I*/

  4: typedef struct {
  5:   IS          isrow,iscol;      /* rows and columns in submatrix, only used to check consistency */
  6:   Vec         left,right;       /* optional scaling */
  7:   Vec         olwork,orwork;    /* work vectors outside the scatters, only touched by PreScale and only created if needed*/
  8:   Vec         lwork,rwork;      /* work vectors inside the scatters */
  9:   VecScatter  lrestrict,rprolong;
 10:   Mat         A;
 11:   PetscScalar scale;
 12: } Mat_SubMatrix;

 16: static PetscErrorCode PreScaleLeft(Mat N,Vec x,Vec *xx)
 17: {
 18:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;

 22:   if (!Na->left) {
 23:     *xx = x;
 24:   } else {
 25:     if (!Na->olwork) {
 26:       VecDuplicate(Na->left,&Na->olwork);
 27:     }
 28:     VecPointwiseMult(Na->olwork,x,Na->left);
 29:     *xx  = Na->olwork;
 30:   }
 31:   return(0);
 32: }

 36: static PetscErrorCode PreScaleRight(Mat N,Vec x,Vec *xx)
 37: {
 38:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;

 42:   if (!Na->right) {
 43:     *xx = x;
 44:   } else {
 45:     if (!Na->orwork) {
 46:       VecDuplicate(Na->right,&Na->orwork);
 47:     }
 48:     VecPointwiseMult(Na->orwork,x,Na->right);
 49:     *xx  = Na->orwork;
 50:   }
 51:   return(0);
 52: }

 56: static PetscErrorCode PostScaleLeft(Mat N,Vec x)
 57: {
 58:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;

 62:   if (Na->left) {
 63:     VecPointwiseMult(x,x,Na->left);
 64:   }
 65:   return(0);
 66: }

 70: static PetscErrorCode PostScaleRight(Mat N,Vec x)
 71: {
 72:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;

 76:   if (Na->right) {
 77:     VecPointwiseMult(x,x,Na->right);
 78:   }
 79:   return(0);
 80: }

 84: static PetscErrorCode MatScale_SubMatrix(Mat N,PetscScalar scale)
 85: {
 86:   Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;

 89:   Na->scale *= scale;
 90:   return(0);
 91: }

 95: static PetscErrorCode MatDiagonalScale_SubMatrix(Mat N,Vec left,Vec right)
 96: {
 97:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;

101:   if (left) {
102:     if (!Na->left) {
103:       VecDuplicate(left,&Na->left);
104:       VecCopy(left,Na->left);
105:     } else {
106:       VecPointwiseMult(Na->left,left,Na->left);
107:     }
108:   }
109:   if (right) {
110:     if (!Na->right) {
111:       VecDuplicate(right,&Na->right);
112:       VecCopy(right,Na->right);
113:     } else {
114:       VecPointwiseMult(Na->right,right,Na->right);
115:     }
116:   }
117:   return(0);
118: }

122: static PetscErrorCode MatMult_SubMatrix(Mat N,Vec x,Vec y)
123: {
124:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;
125:   Vec            xx  = 0;

129:   PreScaleRight(N,x,&xx);
130:   VecZeroEntries(Na->rwork);
131:   VecScatterBegin(Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
132:   VecScatterEnd  (Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
133:   MatMult(Na->A,Na->rwork,Na->lwork);
134:   VecScatterBegin(Na->lrestrict,Na->lwork,y,INSERT_VALUES,SCATTER_FORWARD);
135:   VecScatterEnd  (Na->lrestrict,Na->lwork,y,INSERT_VALUES,SCATTER_FORWARD);
136:   PostScaleLeft(N,y);
137:   VecScale(y,Na->scale);
138:   return(0);
139: }

143: static PetscErrorCode MatMultAdd_SubMatrix(Mat N,Vec v1,Vec v2,Vec v3)
144: {
145:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;
146:   Vec            xx  = 0;

150:   PreScaleRight(N,v1,&xx);
151:   VecZeroEntries(Na->rwork);
152:   VecScatterBegin(Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
153:   VecScatterEnd  (Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
154:   MatMult(Na->A,Na->rwork,Na->lwork);
155:   if (v2 == v3) {
156:     if (Na->scale == (PetscScalar)1.0 && !Na->left) {
157:       VecScatterBegin(Na->lrestrict,Na->lwork,v3,ADD_VALUES,SCATTER_FORWARD);
158:       VecScatterEnd  (Na->lrestrict,Na->lwork,v3,ADD_VALUES,SCATTER_FORWARD);
159:     } else {
160:       if (!Na->olwork) {VecDuplicate(v3,&Na->olwork);}
161:       VecScatterBegin(Na->lrestrict,Na->lwork,Na->olwork,INSERT_VALUES,SCATTER_FORWARD);
162:       VecScatterEnd  (Na->lrestrict,Na->lwork,Na->olwork,INSERT_VALUES,SCATTER_FORWARD);
163:       PostScaleLeft(N,Na->olwork);
164:       VecAXPY(v3,Na->scale,Na->olwork);
165:     }
166:   } else {
167:     VecScatterBegin(Na->lrestrict,Na->lwork,v3,INSERT_VALUES,SCATTER_FORWARD);
168:     VecScatterEnd  (Na->lrestrict,Na->lwork,v3,INSERT_VALUES,SCATTER_FORWARD);
169:     PostScaleLeft(N,v3);
170:     VecAYPX(v3,Na->scale,v2);
171:   }
172:   return(0);
173: }

177: static PetscErrorCode MatMultTranspose_SubMatrix(Mat N,Vec x,Vec y)
178: {
179:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;
180:   Vec            xx  = 0;

184:   PreScaleLeft(N,x,&xx);
185:   VecZeroEntries(Na->lwork);
186:   VecScatterBegin(Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
187:   VecScatterEnd  (Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
188:   MatMultTranspose(Na->A,Na->lwork,Na->rwork);
189:   VecScatterBegin(Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);
190:   VecScatterEnd  (Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);
191:   PostScaleRight(N,y);
192:   VecScale(y,Na->scale);
193:   return(0);
194: }

198: static PetscErrorCode MatMultTransposeAdd_SubMatrix(Mat N,Vec v1,Vec v2,Vec v3)
199: {
200:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;
201:   Vec            xx  = 0;

205:   PreScaleLeft(N,v1,&xx);
206:   VecZeroEntries(Na->lwork);
207:   VecScatterBegin(Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
208:   VecScatterEnd  (Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
209:   MatMultTranspose(Na->A,Na->lwork,Na->rwork);
210:   if (v2 == v3) {
211:     if (Na->scale == (PetscScalar)1.0 && !Na->right) {
212:       VecScatterBegin(Na->rprolong,Na->rwork,v3,ADD_VALUES,SCATTER_REVERSE);
213:       VecScatterEnd  (Na->rprolong,Na->rwork,v3,ADD_VALUES,SCATTER_REVERSE);
214:     } else {
215:       if (!Na->orwork) {VecDuplicate(v3,&Na->orwork);}
216:       VecScatterBegin(Na->rprolong,Na->rwork,Na->orwork,INSERT_VALUES,SCATTER_REVERSE);
217:       VecScatterEnd  (Na->rprolong,Na->rwork,Na->orwork,INSERT_VALUES,SCATTER_REVERSE);
218:       PostScaleRight(N,Na->orwork);
219:       VecAXPY(v3,Na->scale,Na->orwork);
220:     }
221:   } else {
222:     VecScatterBegin(Na->rprolong,Na->rwork,v3,INSERT_VALUES,SCATTER_REVERSE);
223:     VecScatterEnd  (Na->rprolong,Na->rwork,v3,INSERT_VALUES,SCATTER_REVERSE);
224:     PostScaleRight(N,v3);
225:     VecAYPX(v3,Na->scale,v2);
226:   }
227:   return(0);
228: }

232: static PetscErrorCode MatDestroy_SubMatrix(Mat N)
233: {
234:   Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;

238:   ISDestroy(&Na->isrow);
239:   ISDestroy(&Na->iscol);
240:   VecDestroy(&Na->left);
241:   VecDestroy(&Na->right);
242:   VecDestroy(&Na->olwork);
243:   VecDestroy(&Na->orwork);
244:   VecDestroy(&Na->lwork);
245:   VecDestroy(&Na->rwork);
246:   VecScatterDestroy(&Na->lrestrict);
247:   VecScatterDestroy(&Na->rprolong);
248:   MatDestroy(&Na->A);
249:   PetscFree(N->data);
250:   return(0);
251: }

255: /*@
256:    MatCreateSubMatrix - Creates a composite matrix that acts as a submatrix

258:    Collective on Mat

260:    Input Parameters:
261: +  A - matrix that we will extract a submatrix of
262: .  isrow - rows to be present in the submatrix
263: -  iscol - columns to be present in the submatrix

265:    Output Parameters:
266: .  newmat - new matrix

268:    Level: developer

270:    Notes:
271:    Most will use MatGetSubMatrix which provides a more efficient representation if it is available.

273: .seealso: MatGetSubMatrix(), MatSubMatrixUpdate()
274: @*/
275: PetscErrorCode  MatCreateSubMatrix(Mat A,IS isrow,IS iscol,Mat *newmat)
276: {
277:   Vec            left,right;
278:   PetscInt       m,n;
279:   Mat            N;
280:   Mat_SubMatrix  *Na;

288:   *newmat = 0;

290:   MatCreate(PetscObjectComm((PetscObject)A),&N);
291:   ISGetLocalSize(isrow,&m);
292:   ISGetLocalSize(iscol,&n);
293:   MatSetSizes(N,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
294:   PetscObjectChangeTypeName((PetscObject)N,MATSUBMATRIX);

296:   PetscNewLog(N,Mat_SubMatrix,&Na);
297:   N->data   = (void*)Na;
298:   PetscObjectReference((PetscObject)A);
299:   PetscObjectReference((PetscObject)isrow);
300:   PetscObjectReference((PetscObject)iscol);
301:   Na->A     = A;
302:   Na->isrow = isrow;
303:   Na->iscol = iscol;
304:   Na->scale = 1.0;

306:   N->ops->destroy          = MatDestroy_SubMatrix;
307:   N->ops->mult             = MatMult_SubMatrix;
308:   N->ops->multadd          = MatMultAdd_SubMatrix;
309:   N->ops->multtranspose    = MatMultTranspose_SubMatrix;
310:   N->ops->multtransposeadd = MatMultTransposeAdd_SubMatrix;
311:   N->ops->scale            = MatScale_SubMatrix;
312:   N->ops->diagonalscale    = MatDiagonalScale_SubMatrix;

314:   PetscLayoutSetBlockSize(N->rmap,A->rmap->bs);
315:   PetscLayoutSetBlockSize(N->cmap,A->cmap->bs);
316:   PetscLayoutSetUp(N->rmap);
317:   PetscLayoutSetUp(N->cmap);

319:   MatGetVecs(A,&Na->rwork,&Na->lwork);
320:   VecCreate(PetscObjectComm((PetscObject)isrow),&left);
321:   VecCreate(PetscObjectComm((PetscObject)iscol),&right);
322:   VecSetSizes(left,m,PETSC_DETERMINE);
323:   VecSetSizes(right,n,PETSC_DETERMINE);
324:   VecSetUp(left);
325:   VecSetUp(right);
326:   VecScatterCreate(Na->lwork,isrow,left,NULL,&Na->lrestrict);
327:   VecScatterCreate(right,NULL,Na->rwork,iscol,&Na->rprolong);
328:   VecDestroy(&left);
329:   VecDestroy(&right);

331:   N->assembled = PETSC_TRUE;

333:   MatSetUp(N);

335:   *newmat      = N;
336:   return(0);
337: }


342: /*@
343:    MatSubMatrixUpdate - Updates a submatrix

345:    Collective on Mat

347:    Input Parameters:
348: +  N - submatrix to update
349: .  A - full matrix in the submatrix
350: .  isrow - rows in the update (same as the first time the submatrix was created)
351: -  iscol - columns in the update (same as the first time the submatrix was created)

353:    Level: developer

355:    Notes:
356:    Most will use MatGetSubMatrix which provides a more efficient representation if it is available.

358: .seealso: MatGetSubMatrix(), MatCreateSubMatrix()
359: @*/
360: PetscErrorCode  MatSubMatrixUpdate(Mat N,Mat A,IS isrow,IS iscol)
361: {
363:   PetscBool      flg;
364:   Mat_SubMatrix  *Na;

371:   PetscObjectTypeCompare((PetscObject)N,MATSUBMATRIX,&flg);
372:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix has wrong type");

374:   Na   = (Mat_SubMatrix*)N->data;
375:   ISEqual(isrow,Na->isrow,&flg);
376:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different row indices");
377:   ISEqual(iscol,Na->iscol,&flg);
378:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different column indices");

380:   PetscObjectReference((PetscObject)A);
381:   MatDestroy(&Na->A);
382:   Na->A = A;

384:   Na->scale = 1.0;
385:   VecDestroy(&Na->left);
386:   VecDestroy(&Na->right);
387:   return(0);
388: }