xref: /petsc/src/tao/matrix/submatfree.c (revision c77c71ff2d9eaa2c74538bf9bf94eff01b512dbf) 
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 .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;
30 
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));
39 
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 
62   PetscFunctionReturn(PETSC_SUCCESS);
63 }
64 
65 PetscErrorCode MatSMFResetRowColumn(Mat mat, IS Rows, IS Cols)
66 {
67   MatSubMatFreeCtx ctx;
68 
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 }
79 
80 PetscErrorCode MatMult_SMF(Mat mat, Vec a, Vec y)
81 {
82   MatSubMatFreeCtx ctx;
83 
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 }
92 
93 PetscErrorCode MatMultTranspose_SMF(Mat mat, Vec a, Vec y)
94 {
95   MatSubMatFreeCtx ctx;
96 
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 }
105 
106 PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D, InsertMode is)
107 {
108   MatSubMatFreeCtx ctx;
109 
110   PetscFunctionBegin;
111   PetscCall(MatShellGetContext(M, &ctx));
112   PetscCall(MatDiagonalSet(ctx->A, D, is));
113   PetscFunctionReturn(PETSC_SUCCESS);
114 }
115 
116 PetscErrorCode MatDestroy_SMF(Mat mat)
117 {
118   MatSubMatFreeCtx ctx;
119 
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 }
129 
130 PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer)
131 {
132   MatSubMatFreeCtx ctx;
133 
134   PetscFunctionBegin;
135   PetscCall(MatShellGetContext(mat, &ctx));
136   PetscCall(MatView(ctx->A, viewer));
137   PetscFunctionReturn(PETSC_SUCCESS);
138 }
139 
140 PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
141 {
142   MatSubMatFreeCtx ctx;
143 
144   PetscFunctionBegin;
145   PetscCall(MatShellGetContext(Y, &ctx));
146   PetscCall(MatShift(ctx->A, a));
147   PetscFunctionReturn(PETSC_SUCCESS);
148 }
149 
150 PetscErrorCode MatDuplicate_SMF(Mat mat, MatDuplicateOption op, Mat *M)
151 {
152   MatSubMatFreeCtx ctx;
153 
154   PetscFunctionBegin;
155   PetscCall(MatShellGetContext(mat, &ctx));
156   PetscCall(MatCreateSubMatrixFree(ctx->A, ctx->Rows, ctx->Cols, M));
157   PetscFunctionReturn(PETSC_SUCCESS);
158 }
159 
160 PetscErrorCode MatEqual_SMF(Mat A, Mat B, PetscBool *flg)
161 {
162   MatSubMatFreeCtx ctx1, ctx2;
163   PetscBool        flg1, flg2, flg3;
164 
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 }
182 
183 PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
184 {
185   MatSubMatFreeCtx ctx;
186 
187   PetscFunctionBegin;
188   PetscCall(MatShellGetContext(mat, &ctx));
189   PetscCall(MatScale(ctx->A, a));
190   PetscFunctionReturn(PETSC_SUCCESS);
191 }
192 
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 }
199 
200 PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v)
201 {
202   MatSubMatFreeCtx ctx;
203 
204   PetscFunctionBegin;
205   PetscCall(MatShellGetContext(mat, &ctx));
206   PetscCall(MatGetDiagonal(ctx->A, v));
207   PetscFunctionReturn(PETSC_SUCCESS);
208 }
209 
210 PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
211 {
212   MatSubMatFreeCtx ctx;
213 
214   PetscFunctionBegin;
215   PetscCall(MatShellGetContext(M, &ctx));
216   PetscCall(MatGetRowMax(ctx->A, D, NULL));
217   PetscFunctionReturn(PETSC_SUCCESS);
218 }
219 
220 PetscErrorCode MatCreateSubMatrices_SMF(Mat A, PetscInt n, IS *irow, IS *icol, MatReuse scall, Mat **B)
221 {
222   PetscInt i;
223 
224   PetscFunctionBegin;
225   if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));
226 
227   for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SMF(A, irow[i], icol[i], scall, &(*B)[i]));
228   PetscFunctionReturn(PETSC_SUCCESS);
229 }
230 
231 PetscErrorCode MatCreateSubMatrix_SMF(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
232 {
233   MatSubMatFreeCtx ctx;
234 
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 }
241 
242 PetscErrorCode MatGetRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
243 {
244   MatSubMatFreeCtx ctx;
245 
246   PetscFunctionBegin;
247   PetscCall(MatShellGetContext(mat, &ctx));
248   PetscCall(MatGetRow(ctx->A, row, ncols, cols, vals));
249   PetscFunctionReturn(PETSC_SUCCESS);
250 }
251 
252 PetscErrorCode MatRestoreRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
253 {
254   MatSubMatFreeCtx ctx;
255 
256   PetscFunctionBegin;
257   PetscCall(MatShellGetContext(mat, &ctx));
258   PetscCall(MatRestoreRow(ctx->A, row, ncols, cols, vals));
259   PetscFunctionReturn(PETSC_SUCCESS);
260 }
261 
262 PetscErrorCode MatGetColumnVector_SMF(Mat mat, Vec Y, PetscInt col)
263 {
264   MatSubMatFreeCtx ctx;
265 
266   PetscFunctionBegin;
267   PetscCall(MatShellGetContext(mat, &ctx));
268   PetscCall(MatGetColumnVector(ctx->A, Y, col));
269   PetscFunctionReturn(PETSC_SUCCESS);
270 }
271 
272 PetscErrorCode MatNorm_SMF(Mat mat, NormType type, PetscReal *norm)
273 {
274   MatSubMatFreeCtx ctx;
275 
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 }
285