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

 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: .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:   PetscInt         mloc, nloc, m, n;

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

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

 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

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

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

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

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

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

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

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

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

116: PetscErrorCode MatDestroy_SMF(Mat mat)
117: {
118:   MatSubMatFreeCtx ctx;

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

130: PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer)
131: {
132:   MatSubMatFreeCtx ctx;

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

140: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
141: {
142:   MatSubMatFreeCtx ctx;

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

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

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

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

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

183: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
184: {
185:   MatSubMatFreeCtx ctx;

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

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

200: PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v)
201: {
202:   MatSubMatFreeCtx ctx;

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

210: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
211: {
212:   MatSubMatFreeCtx ctx;

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

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

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

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

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

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

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

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

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

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

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

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

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

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