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=mat->comm;
     34:   MatSubMatFreeCtx ctx;
     35:   int          info,mloc,nloc,m,n;

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

     41:   ctx->A=mat;
     42:   //  ctx->Row=Row;
     43:   //  ctx->Col=Col;
     44:   ctx->ColComplement=ColMask;
     45:   ctx->RowComplement=RowMask;

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

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

     54:   ctx->VR=ctx->VC;
     55:   info =  PetscObjectReference((PetscObject)mat);CHKERRQ(info);
     56:   info =  PetscObjectReference((PetscObject)RowMask);CHKERRQ(info);
     57:   info =  PetscObjectReference((PetscObject)ColMask);CHKERRQ(info);

     59:   info = MatCreateShell(comm,mloc,nloc,m,n,ctx,J);CHKERRQ(info);

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

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

     81:   return(0);
     82: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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