1 #include <petsctao.h> /*I "petsctao.h" I*/
2 #include <../src/tao/matrix/submatfree.h>
3
4 /*@
5 MatCreateSubMatrixFree - Creates a reduced matrix by masking a
6 full matrix.
7
8 Collective
9
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
14
15 Output Parameter:
16 . J - New matrix
17
18 Level: developer
19
20 Note:
21 The caller is responsible for destroying the input objects after matrix J has been destroyed.
22
23 Developer Note:
24 This should be moved/supported in `Mat`
25
26 .seealso: `MatCreate()`
27 @*/
MatCreateSubMatrixFree(Mat mat,IS Rows,IS Cols,Mat * J)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;
33
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));
42
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 }
66
MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols)67 PetscErrorCode MatSMFResetRowColumn(Mat mat, IS Rows, IS Cols)
68 {
69 MatSubMatFreeCtx ctx;
70
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 }
81
MatMult_SMF(Mat mat,Vec a,Vec y)82 PetscErrorCode MatMult_SMF(Mat mat, Vec a, Vec y)
83 {
84 MatSubMatFreeCtx ctx;
85
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 }
94
MatMultTranspose_SMF(Mat mat,Vec a,Vec y)95 PetscErrorCode MatMultTranspose_SMF(Mat mat, Vec a, Vec y)
96 {
97 MatSubMatFreeCtx ctx;
98
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 }
107
MatDiagonalSet_SMF(Mat M,Vec D,InsertMode is)108 PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D, InsertMode is)
109 {
110 MatSubMatFreeCtx ctx;
111
112 PetscFunctionBegin;
113 PetscCall(MatShellGetContext(M, &ctx));
114 PetscCall(MatDiagonalSet(ctx->A, D, is));
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
MatDestroy_SMF(Mat mat)118 PetscErrorCode MatDestroy_SMF(Mat mat)
119 {
120 MatSubMatFreeCtx ctx;
121
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 }
131
MatView_SMF(Mat mat,PetscViewer viewer)132 PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer)
133 {
134 MatSubMatFreeCtx ctx;
135
136 PetscFunctionBegin;
137 PetscCall(MatShellGetContext(mat, &ctx));
138 PetscCall(MatView(ctx->A, viewer));
139 PetscFunctionReturn(PETSC_SUCCESS);
140 }
141
MatShift_SMF(Mat Y,PetscReal a)142 PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
143 {
144 MatSubMatFreeCtx ctx;
145
146 PetscFunctionBegin;
147 PetscCall(MatShellGetContext(Y, &ctx));
148 PetscCall(MatShift(ctx->A, a));
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat * M)152 PetscErrorCode MatDuplicate_SMF(Mat mat, MatDuplicateOption op, Mat *M)
153 {
154 MatSubMatFreeCtx ctx;
155
156 PetscFunctionBegin;
157 PetscCall(MatShellGetContext(mat, &ctx));
158 PetscCall(MatCreateSubMatrixFree(ctx->A, ctx->Rows, ctx->Cols, M));
159 PetscFunctionReturn(PETSC_SUCCESS);
160 }
161
MatEqual_SMF(Mat A,Mat B,PetscBool * flg)162 PetscErrorCode MatEqual_SMF(Mat A, Mat B, PetscBool *flg)
163 {
164 MatSubMatFreeCtx ctx1, ctx2;
165 PetscBool flg1, flg2, flg3;
166
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 }
184
MatScale_SMF(Mat mat,PetscReal a)185 PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
186 {
187 MatSubMatFreeCtx ctx;
188
189 PetscFunctionBegin;
190 PetscCall(MatShellGetContext(mat, &ctx));
191 PetscCall(MatScale(ctx->A, a));
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
MatTranspose_SMF(Mat mat,Mat * B)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 }
201
MatGetDiagonal_SMF(Mat mat,Vec v)202 PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v)
203 {
204 MatSubMatFreeCtx ctx;
205
206 PetscFunctionBegin;
207 PetscCall(MatShellGetContext(mat, &ctx));
208 PetscCall(MatGetDiagonal(ctx->A, v));
209 PetscFunctionReturn(PETSC_SUCCESS);
210 }
211
MatGetRowMax_SMF(Mat M,Vec D)212 PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
213 {
214 MatSubMatFreeCtx ctx;
215
216 PetscFunctionBegin;
217 PetscCall(MatShellGetContext(M, &ctx));
218 PetscCall(MatGetRowMax(ctx->A, D, NULL));
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
MatCreateSubMatrices_SMF(Mat A,PetscInt n,IS * irow,IS * icol,MatReuse scall,Mat ** B)222 PetscErrorCode MatCreateSubMatrices_SMF(Mat A, PetscInt n, IS *irow, IS *icol, MatReuse scall, Mat **B)
223 {
224 PetscInt i;
225
226 PetscFunctionBegin;
227 if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));
228
229 for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SMF(A, irow[i], icol[i], scall, &(*B)[i]));
230 PetscFunctionReturn(PETSC_SUCCESS);
231 }
232
MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,Mat * newmat)233 PetscErrorCode MatCreateSubMatrix_SMF(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
234 {
235 MatSubMatFreeCtx ctx;
236
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 }
243
MatGetRow_SMF(Mat mat,PetscInt row,PetscInt * ncols,const PetscInt ** cols,const PetscScalar ** vals)244 PetscErrorCode MatGetRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
245 {
246 MatSubMatFreeCtx ctx;
247
248 PetscFunctionBegin;
249 PetscCall(MatShellGetContext(mat, &ctx));
250 PetscCall(MatGetRow(ctx->A, row, ncols, cols, vals));
251 PetscFunctionReturn(PETSC_SUCCESS);
252 }
253
MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt * ncols,const PetscInt ** cols,const PetscScalar ** vals)254 PetscErrorCode MatRestoreRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
255 {
256 MatSubMatFreeCtx ctx;
257
258 PetscFunctionBegin;
259 PetscCall(MatShellGetContext(mat, &ctx));
260 PetscCall(MatRestoreRow(ctx->A, row, ncols, cols, vals));
261 PetscFunctionReturn(PETSC_SUCCESS);
262 }
263
MatGetColumnVector_SMF(Mat mat,Vec Y,PetscInt col)264 PetscErrorCode MatGetColumnVector_SMF(Mat mat, Vec Y, PetscInt col)
265 {
266 MatSubMatFreeCtx ctx;
267
268 PetscFunctionBegin;
269 PetscCall(MatShellGetContext(mat, &ctx));
270 PetscCall(MatGetColumnVector(ctx->A, Y, col));
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
MatNorm_SMF(Mat mat,NormType type,PetscReal * norm)274 PetscErrorCode MatNorm_SMF(Mat mat, NormType type, PetscReal *norm)
275 {
276 MatSubMatFreeCtx ctx;
277
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 }
287