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: }