Actual source code: submatfree.c

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

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

      8:   Collective

     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 Parameter:
     16: . J - New matrix

     18:   Level: developer

     20:   Note:
     21:   The caller is responsible for destroying the input objects after matrix J has been destroyed.

     23:   Developer Note:
     24:   This should be moved/supported in `Mat`

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

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

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

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

     71:   PetscFunctionBegin;
     72:   PetscCall(MatShellGetContext(mat, &ctx));
     73:   PetscCall(ISDestroy(&ctx->Rows));
     74:   PetscCall(ISDestroy(&ctx->Cols));
     75:   PetscCall(PetscObjectReference((PetscObject)Rows));
     76:   PetscCall(PetscObjectReference((PetscObject)Cols));
     77:   ctx->Cols = Cols;
     78:   ctx->Rows = Rows;
     79:   PetscFunctionReturn(PETSC_SUCCESS);
     80: }

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

     86:   PetscFunctionBegin;
     87:   PetscCall(MatShellGetContext(mat, &ctx));
     88:   PetscCall(VecCopy(a, ctx->VR));
     89:   PetscCall(VecISSet(ctx->VR, ctx->Cols, 0.0));
     90:   PetscCall(MatMult(ctx->A, ctx->VR, y));
     91:   PetscCall(VecISSet(y, ctx->Rows, 0.0));
     92:   PetscFunctionReturn(PETSC_SUCCESS);
     93: }

     95: PetscErrorCode MatMultTranspose_SMF(Mat mat, Vec a, Vec y)
     96: {
     97:   MatSubMatFreeCtx ctx;

     99:   PetscFunctionBegin;
    100:   PetscCall(MatShellGetContext(mat, &ctx));
    101:   PetscCall(VecCopy(a, ctx->VC));
    102:   PetscCall(VecISSet(ctx->VC, ctx->Rows, 0.0));
    103:   PetscCall(MatMultTranspose(ctx->A, ctx->VC, y));
    104:   PetscCall(VecISSet(y, ctx->Cols, 0.0));
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D, InsertMode is)
    109: {
    110:   MatSubMatFreeCtx ctx;

    112:   PetscFunctionBegin;
    113:   PetscCall(MatShellGetContext(M, &ctx));
    114:   PetscCall(MatDiagonalSet(ctx->A, D, is));
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: PetscErrorCode MatDestroy_SMF(Mat mat)
    119: {
    120:   MatSubMatFreeCtx ctx;

    122:   PetscFunctionBegin;
    123:   PetscCall(MatShellGetContext(mat, &ctx));
    124:   PetscCall(MatDestroy(&ctx->A));
    125:   PetscCall(ISDestroy(&ctx->Rows));
    126:   PetscCall(ISDestroy(&ctx->Cols));
    127:   PetscCall(VecDestroy(&ctx->VC));
    128:   PetscCall(PetscFree(ctx));
    129:   PetscFunctionReturn(PETSC_SUCCESS);
    130: }

    132: PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer)
    133: {
    134:   MatSubMatFreeCtx ctx;

    136:   PetscFunctionBegin;
    137:   PetscCall(MatShellGetContext(mat, &ctx));
    138:   PetscCall(MatView(ctx->A, viewer));
    139:   PetscFunctionReturn(PETSC_SUCCESS);
    140: }

    142: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
    143: {
    144:   MatSubMatFreeCtx ctx;

    146:   PetscFunctionBegin;
    147:   PetscCall(MatShellGetContext(Y, &ctx));
    148:   PetscCall(MatShift(ctx->A, a));
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: PetscErrorCode MatDuplicate_SMF(Mat mat, MatDuplicateOption op, Mat *M)
    153: {
    154:   MatSubMatFreeCtx ctx;

    156:   PetscFunctionBegin;
    157:   PetscCall(MatShellGetContext(mat, &ctx));
    158:   PetscCall(MatCreateSubMatrixFree(ctx->A, ctx->Rows, ctx->Cols, M));
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: PetscErrorCode MatEqual_SMF(Mat A, Mat B, PetscBool *flg)
    163: {
    164:   MatSubMatFreeCtx ctx1, ctx2;
    165:   PetscBool        flg1, flg2, flg3;

    167:   PetscFunctionBegin;
    168:   PetscCall(MatShellGetContext(A, &ctx1));
    169:   PetscCall(MatShellGetContext(B, &ctx2));
    170:   PetscCall(ISEqual(ctx1->Rows, ctx2->Rows, &flg2));
    171:   PetscCall(ISEqual(ctx1->Cols, ctx2->Cols, &flg3));
    172:   if (flg2 == PETSC_FALSE || flg3 == PETSC_FALSE) {
    173:     *flg = PETSC_FALSE;
    174:   } else {
    175:     PetscCall(MatEqual(ctx1->A, ctx2->A, &flg1));
    176:     if (flg1 == PETSC_FALSE) {
    177:       *flg = PETSC_FALSE;
    178:     } else {
    179:       *flg = PETSC_TRUE;
    180:     }
    181:   }
    182:   PetscFunctionReturn(PETSC_SUCCESS);
    183: }

    185: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
    186: {
    187:   MatSubMatFreeCtx ctx;

    189:   PetscFunctionBegin;
    190:   PetscCall(MatShellGetContext(mat, &ctx));
    191:   PetscCall(MatScale(ctx->A, a));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: PetscErrorCode MatTranspose_SMF(Mat mat, Mat *B)
    196: {
    197:   PetscFunctionBegin;
    198:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support for transpose for MatCreateSubMatrixFree() matrix");
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v)
    203: {
    204:   MatSubMatFreeCtx ctx;

    206:   PetscFunctionBegin;
    207:   PetscCall(MatShellGetContext(mat, &ctx));
    208:   PetscCall(MatGetDiagonal(ctx->A, v));
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
    213: {
    214:   MatSubMatFreeCtx ctx;

    216:   PetscFunctionBegin;
    217:   PetscCall(MatShellGetContext(M, &ctx));
    218:   PetscCall(MatGetRowMax(ctx->A, D, NULL));
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }

    222: PetscErrorCode MatCreateSubMatrices_SMF(Mat A, PetscInt n, IS *irow, IS *icol, MatReuse scall, Mat **B)
    223: {
    224:   PetscInt i;

    226:   PetscFunctionBegin;
    227:   if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));

    229:   for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SMF(A, irow[i], icol[i], scall, &(*B)[i]));
    230:   PetscFunctionReturn(PETSC_SUCCESS);
    231: }

    233: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
    234: {
    235:   MatSubMatFreeCtx ctx;

    237:   PetscFunctionBegin;
    238:   PetscCall(MatShellGetContext(mat, &ctx));
    239:   if (newmat) PetscCall(MatDestroy(&*newmat));
    240:   PetscCall(MatCreateSubMatrixFree(ctx->A, isrow, iscol, newmat));
    241:   PetscFunctionReturn(PETSC_SUCCESS);
    242: }

    244: PetscErrorCode MatGetRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
    245: {
    246:   MatSubMatFreeCtx ctx;

    248:   PetscFunctionBegin;
    249:   PetscCall(MatShellGetContext(mat, &ctx));
    250:   PetscCall(MatGetRow(ctx->A, row, ncols, cols, vals));
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: PetscErrorCode MatRestoreRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
    255: {
    256:   MatSubMatFreeCtx ctx;

    258:   PetscFunctionBegin;
    259:   PetscCall(MatShellGetContext(mat, &ctx));
    260:   PetscCall(MatRestoreRow(ctx->A, row, ncols, cols, vals));
    261:   PetscFunctionReturn(PETSC_SUCCESS);
    262: }

    264: PetscErrorCode MatGetColumnVector_SMF(Mat mat, Vec Y, PetscInt col)
    265: {
    266:   MatSubMatFreeCtx ctx;

    268:   PetscFunctionBegin;
    269:   PetscCall(MatShellGetContext(mat, &ctx));
    270:   PetscCall(MatGetColumnVector(ctx->A, Y, col));
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: PetscErrorCode MatNorm_SMF(Mat mat, NormType type, PetscReal *norm)
    275: {
    276:   MatSubMatFreeCtx ctx;

    278:   PetscFunctionBegin;
    279:   PetscCall(MatShellGetContext(mat, &ctx));
    280:   if (type == NORM_FROBENIUS) {
    281:     *norm = 1.0;
    282:   } else if (type == NORM_1 || type == NORM_INFINITY) {
    283:     *norm = 1.0;
    284:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }