Actual source code: submatfree.c

      1: #include <petsctao.h>
      2: #include <../src/tao/matrix/submatfree.h>

      4: /*@C
      5:   MatCreateSubMatrixFree - Creates a reduced matrix by masking a
      6:   full matrix.

      8:    Collective on matrix

     10:    Input Parameters:
     11: +  mat - matrix of arbitrary type
     12: .  Rows - the rows that will be in the submatrix
     13: -  Cols - the columns that will be in the submatrix

     15:    Output Parameters:
     16: .  J - New matrix

     18:    Notes:
     19:    The caller is responsible for destroying the input objects after matrix J has been destroyed.

     21:    Level: developer

     23: .seealso: MatCreate()
     24: @*/
     25: PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J)
     26: {
     27:   MPI_Comm         comm=PetscObjectComm((PetscObject)mat);
     28:   MatSubMatFreeCtx ctx;
     29:   PetscErrorCode   ierr;
     30:   PetscInt         mloc,nloc,m,n;

     33:   PetscNew(&ctx);
     34:   ctx->A  = mat;
     35:   MatGetSize(mat,&m,&n);
     36:   MatGetLocalSize(mat,&mloc,&nloc);
     37:   MatCreateVecs(mat,NULL,&ctx->VC);
     38:   ctx->VR = ctx->VC;
     39:    PetscObjectReference((PetscObject)mat);

     42:   ctx->Rows = Rows;
     43:   ctx->Cols = Cols;
     44:   PetscObjectReference((PetscObject)Rows);
     45:   PetscObjectReference((PetscObject)Cols);
     46:   MatCreateShell(comm,mloc,nloc,m,n,ctx,J);
     47:   MatShellSetManageScalingShifts(*J);
     48:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);
     49:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);
     50:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);
     51:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);
     52:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);
     53:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);
     54:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);
     55:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);
     56:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);
     57:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);
     58:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_SMF);
     59:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);
     60:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);
     61:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_SMF);
     62:   MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);

     64:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
     65:   return(0);
     66: }

     68: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){
     69:   MatSubMatFreeCtx ctx;
     70:   PetscErrorCode   ierr;

     73:   MatShellGetContext(mat,(void **)&ctx);
     74:   ISDestroy(&ctx->Rows);
     75:   ISDestroy(&ctx->Cols);
     76:   PetscObjectReference((PetscObject)Rows);
     77:   PetscObjectReference((PetscObject)Cols);
     78:   ctx->Cols=Cols;
     79:   ctx->Rows=Rows;
     80:   return(0);
     81: }

     83: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
     84: {
     85:   MatSubMatFreeCtx ctx;
     86:   PetscErrorCode   ierr;

     89:   MatShellGetContext(mat,(void **)&ctx);
     90:   VecCopy(a,ctx->VR);
     91:   VecISSet(ctx->VR,ctx->Cols,0.0);
     92:   MatMult(ctx->A,ctx->VR,y);
     93:   VecISSet(y,ctx->Rows,0.0);
     94:   return(0);
     95: }

     97: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
     98: {
     99:   MatSubMatFreeCtx ctx;
    100:   PetscErrorCode   ierr;

    103:   MatShellGetContext(mat,(void **)&ctx);
    104:   VecCopy(a,ctx->VC);
    105:   VecISSet(ctx->VC,ctx->Rows,0.0);
    106:   MatMultTranspose(ctx->A,ctx->VC,y);
    107:   VecISSet(y,ctx->Cols,0.0);
    108:   return(0);
    109: }

    111: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
    112: {
    113:   MatSubMatFreeCtx ctx;
    114:   PetscErrorCode   ierr;

    117:   MatShellGetContext(M,(void **)&ctx);
    118:   MatDiagonalSet(ctx->A,D,is);
    119:   return(0);
    120: }

    122: PetscErrorCode MatDestroy_SMF(Mat mat)
    123: {
    124:   PetscErrorCode   ierr;
    125:   MatSubMatFreeCtx ctx;

    128:   MatShellGetContext(mat,(void **)&ctx);
    129:   MatDestroy(&ctx->A);
    130:   ISDestroy(&ctx->Rows);
    131:   ISDestroy(&ctx->Cols);
    132:   VecDestroy(&ctx->VC);
    133:   PetscFree(ctx);
    134:   return(0);
    135: }

    139: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
    140: {
    141:   PetscErrorCode   ierr;
    142:   MatSubMatFreeCtx ctx;

    145:   MatShellGetContext(mat,(void **)&ctx);
    146:   MatView(ctx->A,viewer);
    147:   return(0);
    148: }

    150: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
    151: {
    152:   PetscErrorCode   ierr;
    153:   MatSubMatFreeCtx ctx;

    156:   MatShellGetContext(Y,(void **)&ctx);
    157:   MatShift(ctx->A,a);
    158:   return(0);
    159: }

    161: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
    162: {
    163:   PetscErrorCode   ierr;
    164:   MatSubMatFreeCtx ctx;

    167:   MatShellGetContext(mat,(void **)&ctx);
    168:   MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
    169:   return(0);
    170: }

    172: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
    173: {
    174:   PetscErrorCode    ierr;
    175:   MatSubMatFreeCtx  ctx1,ctx2;
    176:   PetscBool         flg1,flg2,flg3;

    179:   MatShellGetContext(A,(void **)&ctx1);
    180:   MatShellGetContext(B,(void **)&ctx2);
    181:   ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
    182:   ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
    183:   if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
    184:     *flg=PETSC_FALSE;
    185:   } else {
    186:     MatEqual(ctx1->A,ctx2->A,&flg1);
    187:     if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
    188:     else { *flg=PETSC_TRUE;}
    189:   }
    190:   return(0);
    191: }

    193: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
    194: {
    195:   PetscErrorCode   ierr;
    196:   MatSubMatFreeCtx ctx;

    199:   MatShellGetContext(mat,(void **)&ctx);
    200:   MatScale(ctx->A,a);
    201:   return(0);
    202: }

    204: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
    205: {
    207:   PetscFunctionReturn(1);
    208: }

    210: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
    211: {
    212:   PetscErrorCode   ierr;
    213:   MatSubMatFreeCtx ctx;

    216:   MatShellGetContext(mat,(void **)&ctx);
    217:   MatGetDiagonal(ctx->A,v);
    218:   return(0);
    219: }

    221: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
    222: {
    223:   MatSubMatFreeCtx ctx;
    224:   PetscErrorCode   ierr;

    227:   MatShellGetContext(M,(void **)&ctx);
    228:   MatGetRowMax(ctx->A,D,NULL);
    229:   return(0);
    230: }

    232: PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
    233: {
    235:   PetscInt       i;

    238:   if (scall == MAT_INITIAL_MATRIX) {
    239:     PetscCalloc1(n+1,B);
    240:   }

    242:   for (i=0; i<n; i++) {
    243:     MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
    244:   }
    245:   return(0);
    246: }

    248: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
    249:                         Mat *newmat)
    250: {
    251:   PetscErrorCode   ierr;
    252:   MatSubMatFreeCtx ctx;

    255:   MatShellGetContext(mat,(void **)&ctx);
    256:   if (newmat){
    257:     MatDestroy(&*newmat);
    258:   }
    259:   MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
    260:   return(0);
    261: }

    263: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
    264: {
    265:   PetscErrorCode   ierr;
    266:   MatSubMatFreeCtx ctx;

    269:   MatShellGetContext(mat,(void **)&ctx);
    270:   MatGetRow(ctx->A,row,ncols,cols,vals);
    271:   return(0);
    272: }

    274: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
    275: {
    276:   PetscErrorCode   ierr;
    277:   MatSubMatFreeCtx ctx;

    280:   MatShellGetContext(mat,(void **)&ctx);
    281:   MatRestoreRow(ctx->A,row,ncols,cols,vals);
    282:   return(0);
    283: }

    285: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
    286: {
    287:   PetscErrorCode   ierr;
    288:   MatSubMatFreeCtx ctx;

    291:   MatShellGetContext(mat,(void **)&ctx);
    292:   MatGetColumnVector(ctx->A,Y,col);
    293:   return(0);
    294: }

    296: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
    297: {
    298:   PetscErrorCode    ierr;
    299:   MatSubMatFreeCtx  ctx;

    302:   MatShellGetContext(mat,(void **)&ctx);
    303:   if (type == NORM_FROBENIUS) {
    304:     *norm = 1.0;
    305:   } else if (type == NORM_1 || type == NORM_INFINITY) {
    306:     *norm = 1.0;
    307:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
    308:   return(0);
    309: }