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, (void (*)(void))MatMult_SMF));
50: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))MatDestroy_SMF));
51: PetscCall(MatShellSetOperation(*J, MATOP_VIEW, (void (*)(void))MatView_SMF));
52: PetscCall(MatShellSetOperation(*J, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_SMF));
53: PetscCall(MatShellSetOperation(*J, MATOP_DIAGONAL_SET, (void (*)(void))MatDiagonalSet_SMF));
54: PetscCall(MatShellSetOperation(*J, MATOP_SHIFT, (void (*)(void))MatShift_SMF));
55: PetscCall(MatShellSetOperation(*J, MATOP_EQUAL, (void (*)(void))MatEqual_SMF));
56: PetscCall(MatShellSetOperation(*J, MATOP_SCALE, (void (*)(void))MatScale_SMF));
57: PetscCall(MatShellSetOperation(*J, MATOP_TRANSPOSE, (void (*)(void))MatTranspose_SMF));
58: PetscCall(MatShellSetOperation(*J, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_SMF));
59: PetscCall(MatShellSetOperation(*J, MATOP_CREATE_SUBMATRICES, (void (*)(void))MatCreateSubMatrices_SMF));
60: PetscCall(MatShellSetOperation(*J, MATOP_NORM, (void (*)(void))MatNorm_SMF));
61: PetscCall(MatShellSetOperation(*J, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_SMF));
62: PetscCall(MatShellSetOperation(*J, MATOP_CREATE_SUBMATRIX, (void (*)(void))MatCreateSubMatrix_SMF));
63: PetscCall(MatShellSetOperation(*J, MATOP_GET_ROW_MAX, (void (*)(void))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: }