Actual source code: submatfree.c

  1: #include "submatfree.h"                /*I  "mat.h"  I*/

  3: int ISCreateComplement(IS, Vec, IS *);
  4: int VecISSetToConstant(IS, PetscScalar, Vec);


  9: /*@C
 10:   MatCreateSubMatrixFree - Creates a reduced matrix by masking a
 11:   full matrix.

 13:    Collective on matrix

 15:    Input Parameters:
 16: +  mat - matrix of arbitrary type
 17: .  RowMask - the rows that will be zero
 18: -  ColMask - the columns that will be zero

 20:    Output Parameters:
 21: .  J - New matrix

 23:    Notes: 
 24:    The user provides the input data and is responsible for destroying
 25:    this data after matrix J has been destroyed.  
 26:  
 27:    Level: developer

 29: .seealso: MatCreate()
 30: @*/
 31: int MatCreateSubMatrixFree(Mat mat,IS RowMask, IS ColMask, Mat *J)
 32: {
 33:   MPI_Comm     comm=((PetscObject)mat)->comm;
 34:   MatSubMatFreeCtx ctx;
 35:   int          info;
 36:   PetscInt mloc,nloc,m,n;

 39:   info = PetscNew(_p_MatSubMatFreeCtx,&ctx);CHKERRQ(info);

 41:   ctx->A=mat;
 42:   //  ctx->Row=Row;
 43:   //  ctx->Col=Col;

 45:   info = MatGetSize(mat,&m,&n);CHKERRQ(info);
 46:   info = MatGetLocalSize(mat,&mloc,&nloc);CHKERRQ(info);

 48:   info = VecCreateMPI(comm,nloc,n,&ctx->VC);CHKERRQ(info);
 49:   //  info = ISCreateComplement(Col, ctx->VC, &ctx->ColComplement);CHKERRQ(info);
 50:   //  ctx->RowComplement=ctx->ColComplement;
 51:   ctx->VR=ctx->VC;
 52:   info =  PetscObjectReference((PetscObject)mat);CHKERRQ(info);

 54:   info=ISCreateComplement(RowMask,ctx->VC,&ctx->RowComplement);CHKERRQ(info);
 55:   info=ISCreateComplement(ColMask,ctx->VC,&ctx->ColComplement);CHKERRQ(info);
 56:   /*
 57:   info =  PetscObjectReference((PetscObject)ctx->RowComplement);CHKERRQ(info);
 58:   info =  PetscObjectReference((PetscObject)ctx->ColComplement);CHKERRQ(info);
 59:   */
 60:   info = MatCreateShell(comm,mloc,nloc,m,n,ctx,J);CHKERRQ(info);

 62:   //  info = MatShellSetOperation(*J,MATOP_GET_ROW,(void(*)())MatGetRow_SMF);CHKERRQ(info);
 63:   //  info = MatShellSetOperation(*J,MATOP_RESTORE_ROW,(void(*)())MatRestoreRow_SMF);CHKERRQ(info);
 64:   info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_SMF);CHKERRQ(info);
 65:   info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_SMF);CHKERRQ(info);
 66:   info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_SMF);CHKERRQ(info);
 67:   info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_SMF);CHKERRQ(info);
 68:   info = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)())MatDiagonalSet_SMF);CHKERRQ(info);
 69:   info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_SMF);CHKERRQ(info);
 70:   info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_SMF);CHKERRQ(info);
 71:   info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_SMF);CHKERRQ(info);
 72:   info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_SMF);CHKERRQ(info);
 73:   info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_SMF);CHKERRQ(info);
 74:   info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_SMF);CHKERRQ(info);
 75:   info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_SMF);CHKERRQ(info);
 76:   info = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)())MatDuplicate_SMF);CHKERRQ(info);
 77:   info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_SMF);CHKERRQ(info);
 78:   info = MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)())MatDuplicate_SMF);CHKERRQ(info);

 80:   info = PetscLogObjectParent(mat,*J); CHKERRQ(info);

 82:   return(0);  
 83: }

 87: int MatSMFResetRowColumn(Mat mat,IS RowMask,IS ColMask){
 88:   MatSubMatFreeCtx ctx;
 89:   int           info;

 92:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
 93:   info = PetscObjectReference((PetscObject)RowMask);CHKERRQ(info);
 94:   info = PetscObjectReference((PetscObject)ColMask);CHKERRQ(info);
 95:   info = ISDestroy(ctx->RowComplement);CHKERRQ(info);
 96:   info = ISDestroy(ctx->ColComplement);CHKERRQ(info);
 97:   ctx->ColComplement=ColMask;
 98:   ctx->RowComplement=RowMask;
 99:   return(0);  
100: }

104: int MatMult_SMF(Mat mat,Vec a,Vec y)
105: {
106:   MatSubMatFreeCtx ctx;
107:   int           info;

110:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
111:   info = VecCopy(a,ctx->VR);CHKERRQ(info);
112:   info = VecISSetToConstant(ctx->ColComplement,0.0,ctx->VR);CHKERRQ(info);
113:   info = MatMult(ctx->A,ctx->VR,y);CHKERRQ(info);
114:   info = VecISSetToConstant(ctx->RowComplement,0.0,y);CHKERRQ(info);
115:   return(0);
116: }

120: int MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
121: {
122:   MatSubMatFreeCtx ctx;
123:   int info;

126:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
127:   info = VecCopy(a,ctx->VC);CHKERRQ(info);
128:   info = VecISSetToConstant(ctx->RowComplement,0.0,ctx->VC);CHKERRQ(info);
129:   info = MatMultTranspose(ctx->A,ctx->VC,y);CHKERRQ(info);
130:   info = VecISSetToConstant(ctx->ColComplement,0.0,y);CHKERRQ(info);
131:   return(0);
132: } 

136: int MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
137: {
138:   MatSubMatFreeCtx ctx;
139:   int           info;

142:   info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
143:   info = MatDiagonalSet(ctx->A,D,is);
144:   return(0);
145: } 

149: int MatDestroy_SMF(Mat mat)
150: {
151:   int          info;
152:   MatSubMatFreeCtx ctx;

155:   info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
156:   //  info=ISDestroy(ctx->Row);CHKERRQ(info);
157:   //  info=ISDestroy(ctx->Col);CHKERRQ(info);
158:   info =MatDestroy(ctx->A);CHKERRQ(info);
159:   info=ISDestroy(ctx->RowComplement);CHKERRQ(info);
160:   info=ISDestroy(ctx->ColComplement);CHKERRQ(info);
161:   info=VecDestroy(ctx->VC);CHKERRQ(info);
162:   info = PetscFree(ctx); CHKERRQ(info);
163:   return(0);
164: }



170: int MatView_SMF(Mat mat,PetscViewer viewer)
171: {
172:   int          info;
173:   MatSubMatFreeCtx ctx;

176:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
177:   info = MatView(ctx->A,viewer);CHKERRQ(info);
178:   //  info = ISView(ctx->Row,viewer);CHKERRQ(info);
179:   //  info = ISView(ctx->Col,viewer);CHKERRQ(info);
180:   return(0);
181: }

185: int MatShift_SMF(Mat Y, PetscScalar a)
186: {
187:   int          info;
188:   MatSubMatFreeCtx ctx;

191:   info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
192:   info = MatShift(ctx->A,a);CHKERRQ(info);
193:   return(0);
194: }

198: int MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
199: {
200:   int          info;
201:   MatSubMatFreeCtx ctx;

204:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
205:   info = MatCreateSubMatrixFree(ctx->A,ctx->RowComplement,ctx->ColComplement,M);CHKERRQ(info);
206:   return(0);
207: }

211: int MatEqual_SMF(Mat A,Mat B,PetscTruth *flg)
212: {
213:   int          info;
214:   MatSubMatFreeCtx  ctx1,ctx2;
215:   PetscTruth flg1,flg2,flg3;

218:   info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
219:   info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
220:   info = ISEqual(ctx1->RowComplement,ctx2->RowComplement,&flg2);CHKERRQ(info);
221:   info = ISEqual(ctx1->ColComplement,ctx2->ColComplement,&flg3);CHKERRQ(info);
222:   if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
223:     *flg=PETSC_FALSE;
224:   } else {
225:     info = MatEqual(ctx1->A,ctx2->A,&flg1);CHKERRQ(info);
226:     if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;} 
227:     else { *flg=PETSC_TRUE;} 
228:   }
229:   return(0);
230: }

234: int MatScale_SMF(Mat mat, PetscScalar a)
235: {
236:   int          info;
237:   MatSubMatFreeCtx ctx;

240:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
241:   info = MatScale(ctx->A,a);CHKERRQ(info);
242:   return(0);
243: }

247: int MatTranspose_SMF(Mat mat,Mat *B)
248: {
250:   PetscFunctionReturn(1);
251: }

255: int MatGetDiagonal_SMF(Mat mat,Vec v)
256: {
257:   int          info;
258:   MatSubMatFreeCtx ctx;

261:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
262:   info = MatGetDiagonal(ctx->A,v);CHKERRQ(info);
263:   //  info = VecISSetToConstant(ctx->RowComplement,0.0,v);CHKERRQ(info);
264:   return(0);
265: }

269: int MatGetRowMax_SMF(Mat M, Vec D)
270: {
271:   MatSubMatFreeCtx ctx;
272:   int           info;

275:   info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
276:   info = MatGetRowMax(ctx->A,D,PETSC_NULL);
277:   //  info = VecISSetToConstant(ctx->RowComplement,0.0,D);CHKERRQ(info);
278:   return(0);
279: } 


284: int MatGetSubMatrices_SMF(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
285: {
286:   int info,i;

289:   if (scall == MAT_INITIAL_MATRIX) {
290:     info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
291:   }

293:   for ( i=0; i<n; i++ ) {
294:     info = MatGetSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(info);
295:   }
296:   return(0);
297: }

301: int MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
302:                         Mat *newmat)
303: {
304:   int          info;
305:   MatSubMatFreeCtx ctx;

308:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
309:   if (newmat){
310:     info=MatDestroy(*newmat);CHKERRQ(info);
311:   }
312:   info = MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);CHKERRQ(info);
313:   return(0);
314: }

318: int MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
319: {
320:   int          info;
321:   MatSubMatFreeCtx ctx;

324:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
325:   info = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(info);
326:   return(0);
327: }

331: int MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
332: {
333:   int          info;
334:   MatSubMatFreeCtx ctx;

337:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
338:   info = MatRestoreRow(ctx->A,row,ncols,cols,vals);CHKERRQ(info);
339:   return(0);
340: }

344: int MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
345: {
346:   int    info;
347:   MatSubMatFreeCtx  ctx;

350:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
351:   info = MatGetColumnVector(ctx->A,Y,col);CHKERRQ(info);
352:   //  info = VecISSetToConstant(ctx->RowComplement,0.0,Y);CHKERRQ(info);
353:   return(0);
354: }

358: int MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat)
359: {
360:   int info;
361:   int size;
362:   MatSubMatFreeCtx  ctx;

365:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
366:   MPI_Comm_size(((PetscObject)mat)->comm,&size);
367:   PetscFunctionReturn(1);
368: }

372: int MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
373: {
374:   int info;
375:   MatSubMatFreeCtx  ctx;

378:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

380:   if (type == NORM_FROBENIUS) {
381:     *norm = 1.0;
382:   } else if (type == NORM_1 || type == NORM_INFINITY) {
383:     *norm = 1.0;
384:   } else {
385:     SETERRQ(PETSC_ERR_SUP,"No two norm");
386:   }
387:   return(0);
388: }