xref: /petsc/src/tao/matrix/submatfree.c (revision f13dfd9ea68e0ddeee984e65c377a1819eab8a8a) 
1 #include <petsctao.h> /*I "petsctao.h" I*/
2 #include <../src/tao/matrix/submatfree.h>
3 
4 /*@C
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 @*/
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, (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 }
66 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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