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