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));
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: PetscErrorCode MatSMFResetRowColumn(Mat mat, IS Rows, IS Cols)
65: {
66: MatSubMatFreeCtx ctx;
68: PetscFunctionBegin;
69: PetscCall(MatShellGetContext(mat, &ctx));
70: PetscCall(ISDestroy(&ctx->Rows));
71: PetscCall(ISDestroy(&ctx->Cols));
72: PetscCall(PetscObjectReference((PetscObject)Rows));
73: PetscCall(PetscObjectReference((PetscObject)Cols));
74: ctx->Cols = Cols;
75: ctx->Rows = Rows;
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: PetscErrorCode MatMult_SMF(Mat mat, Vec a, Vec y)
80: {
81: MatSubMatFreeCtx ctx;
83: PetscFunctionBegin;
84: PetscCall(MatShellGetContext(mat, &ctx));
85: PetscCall(VecCopy(a, ctx->VR));
86: PetscCall(VecISSet(ctx->VR, ctx->Cols, 0.0));
87: PetscCall(MatMult(ctx->A, ctx->VR, y));
88: PetscCall(VecISSet(y, ctx->Rows, 0.0));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: PetscErrorCode MatMultTranspose_SMF(Mat mat, Vec a, Vec y)
93: {
94: MatSubMatFreeCtx ctx;
96: PetscFunctionBegin;
97: PetscCall(MatShellGetContext(mat, &ctx));
98: PetscCall(VecCopy(a, ctx->VC));
99: PetscCall(VecISSet(ctx->VC, ctx->Rows, 0.0));
100: PetscCall(MatMultTranspose(ctx->A, ctx->VC, y));
101: PetscCall(VecISSet(y, ctx->Cols, 0.0));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D, InsertMode is)
106: {
107: MatSubMatFreeCtx ctx;
109: PetscFunctionBegin;
110: PetscCall(MatShellGetContext(M, &ctx));
111: PetscCall(MatDiagonalSet(ctx->A, D, is));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: PetscErrorCode MatDestroy_SMF(Mat mat)
116: {
117: MatSubMatFreeCtx ctx;
119: PetscFunctionBegin;
120: PetscCall(MatShellGetContext(mat, &ctx));
121: PetscCall(MatDestroy(&ctx->A));
122: PetscCall(ISDestroy(&ctx->Rows));
123: PetscCall(ISDestroy(&ctx->Cols));
124: PetscCall(VecDestroy(&ctx->VC));
125: PetscCall(PetscFree(ctx));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer)
130: {
131: MatSubMatFreeCtx ctx;
133: PetscFunctionBegin;
134: PetscCall(MatShellGetContext(mat, &ctx));
135: PetscCall(MatView(ctx->A, viewer));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
140: {
141: MatSubMatFreeCtx ctx;
143: PetscFunctionBegin;
144: PetscCall(MatShellGetContext(Y, &ctx));
145: PetscCall(MatShift(ctx->A, a));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: PetscErrorCode MatDuplicate_SMF(Mat mat, MatDuplicateOption op, Mat *M)
150: {
151: MatSubMatFreeCtx ctx;
153: PetscFunctionBegin;
154: PetscCall(MatShellGetContext(mat, &ctx));
155: PetscCall(MatCreateSubMatrixFree(ctx->A, ctx->Rows, ctx->Cols, M));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: PetscErrorCode MatEqual_SMF(Mat A, Mat B, PetscBool *flg)
160: {
161: MatSubMatFreeCtx ctx1, ctx2;
162: PetscBool flg1, flg2, flg3;
164: PetscFunctionBegin;
165: PetscCall(MatShellGetContext(A, &ctx1));
166: PetscCall(MatShellGetContext(B, &ctx2));
167: PetscCall(ISEqual(ctx1->Rows, ctx2->Rows, &flg2));
168: PetscCall(ISEqual(ctx1->Cols, ctx2->Cols, &flg3));
169: if (flg2 == PETSC_FALSE || flg3 == PETSC_FALSE) {
170: *flg = PETSC_FALSE;
171: } else {
172: PetscCall(MatEqual(ctx1->A, ctx2->A, &flg1));
173: if (flg1 == PETSC_FALSE) {
174: *flg = PETSC_FALSE;
175: } else {
176: *flg = PETSC_TRUE;
177: }
178: }
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
183: {
184: MatSubMatFreeCtx ctx;
186: PetscFunctionBegin;
187: PetscCall(MatShellGetContext(mat, &ctx));
188: PetscCall(MatScale(ctx->A, a));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: PetscErrorCode MatTranspose_SMF(Mat mat, Mat *B)
193: {
194: PetscFunctionBegin;
195: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support for transpose for MatCreateSubMatrixFree() matrix");
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v)
200: {
201: MatSubMatFreeCtx ctx;
203: PetscFunctionBegin;
204: PetscCall(MatShellGetContext(mat, &ctx));
205: PetscCall(MatGetDiagonal(ctx->A, v));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
210: {
211: MatSubMatFreeCtx ctx;
213: PetscFunctionBegin;
214: PetscCall(MatShellGetContext(M, &ctx));
215: PetscCall(MatGetRowMax(ctx->A, D, NULL));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: PetscErrorCode MatCreateSubMatrices_SMF(Mat A, PetscInt n, IS *irow, IS *icol, MatReuse scall, Mat **B)
220: {
221: PetscInt i;
223: PetscFunctionBegin;
224: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));
226: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SMF(A, irow[i], icol[i], scall, &(*B)[i]));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
231: {
232: MatSubMatFreeCtx ctx;
234: PetscFunctionBegin;
235: PetscCall(MatShellGetContext(mat, &ctx));
236: if (newmat) PetscCall(MatDestroy(&*newmat));
237: PetscCall(MatCreateSubMatrixFree(ctx->A, isrow, iscol, newmat));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: PetscErrorCode MatGetRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
242: {
243: MatSubMatFreeCtx ctx;
245: PetscFunctionBegin;
246: PetscCall(MatShellGetContext(mat, &ctx));
247: PetscCall(MatGetRow(ctx->A, row, ncols, cols, vals));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: PetscErrorCode MatRestoreRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
252: {
253: MatSubMatFreeCtx ctx;
255: PetscFunctionBegin;
256: PetscCall(MatShellGetContext(mat, &ctx));
257: PetscCall(MatRestoreRow(ctx->A, row, ncols, cols, vals));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: PetscErrorCode MatGetColumnVector_SMF(Mat mat, Vec Y, PetscInt col)
262: {
263: MatSubMatFreeCtx ctx;
265: PetscFunctionBegin;
266: PetscCall(MatShellGetContext(mat, &ctx));
267: PetscCall(MatGetColumnVector(ctx->A, Y, col));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: PetscErrorCode MatNorm_SMF(Mat mat, NormType type, PetscReal *norm)
272: {
273: MatSubMatFreeCtx ctx;
275: PetscFunctionBegin;
276: PetscCall(MatShellGetContext(mat, &ctx));
277: if (type == NORM_FROBENIUS) {
278: *norm = 1.0;
279: } else if (type == NORM_1 || type == NORM_INFINITY) {
280: *norm = 1.0;
281: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }