xref: /petsc/src/tao/matrix/submatfree.c (revision 2d30e087755efd99e28fdfe792ffbeb2ee1ea928) 
1 #include <petsctao.h>                     /*I "petsctao.h" I*/
2 #include <../src/tao/matrix/submatfree.h> /*I "submatfree.h" I*/
3 
4 /*@C
5   MatCreateSubMatrixFree - Creates a reduced matrix by masking a
6   full matrix.
7 
8    Collective on matrix
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 Parameters:
16 .  J - New matrix
17 
18    Notes:
19    The caller is responsible for destroying the input objects after matrix J has been destroyed.
20 
21    Level: developer
22 
23 .seealso: `MatCreate()`
24 @*/
25 PetscErrorCode MatCreateSubMatrixFree(Mat mat, IS Rows, IS Cols, Mat *J) {
26   MPI_Comm         comm = PetscObjectComm((PetscObject)mat);
27   MatSubMatFreeCtx ctx;
28   PetscInt         mloc, nloc, m, n;
29 
30   PetscFunctionBegin;
31   PetscCall(PetscNew(&ctx));
32   ctx->A = mat;
33   PetscCall(MatGetSize(mat, &m, &n));
34   PetscCall(MatGetLocalSize(mat, &mloc, &nloc));
35   PetscCall(MatCreateVecs(mat, NULL, &ctx->VC));
36   ctx->VR = ctx->VC;
37   PetscCall(PetscObjectReference((PetscObject)mat));
38 
39   ctx->Rows = Rows;
40   ctx->Cols = Cols;
41   PetscCall(PetscObjectReference((PetscObject)Rows));
42   PetscCall(PetscObjectReference((PetscObject)Cols));
43   PetscCall(MatCreateShell(comm, mloc, nloc, m, n, ctx, J));
44   PetscCall(MatShellSetManageScalingShifts(*J));
45   PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))MatMult_SMF));
46   PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))MatDestroy_SMF));
47   PetscCall(MatShellSetOperation(*J, MATOP_VIEW, (void (*)(void))MatView_SMF));
48   PetscCall(MatShellSetOperation(*J, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_SMF));
49   PetscCall(MatShellSetOperation(*J, MATOP_DIAGONAL_SET, (void (*)(void))MatDiagonalSet_SMF));
50   PetscCall(MatShellSetOperation(*J, MATOP_SHIFT, (void (*)(void))MatShift_SMF));
51   PetscCall(MatShellSetOperation(*J, MATOP_EQUAL, (void (*)(void))MatEqual_SMF));
52   PetscCall(MatShellSetOperation(*J, MATOP_SCALE, (void (*)(void))MatScale_SMF));
53   PetscCall(MatShellSetOperation(*J, MATOP_TRANSPOSE, (void (*)(void))MatTranspose_SMF));
54   PetscCall(MatShellSetOperation(*J, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_SMF));
55   PetscCall(MatShellSetOperation(*J, MATOP_CREATE_SUBMATRICES, (void (*)(void))MatCreateSubMatrices_SMF));
56   PetscCall(MatShellSetOperation(*J, MATOP_NORM, (void (*)(void))MatNorm_SMF));
57   PetscCall(MatShellSetOperation(*J, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_SMF));
58   PetscCall(MatShellSetOperation(*J, MATOP_CREATE_SUBMATRIX, (void (*)(void))MatCreateSubMatrix_SMF));
59   PetscCall(MatShellSetOperation(*J, MATOP_GET_ROW_MAX, (void (*)(void))MatDuplicate_SMF));
60 
61   PetscCall(PetscLogObjectParent((PetscObject)mat, (PetscObject)(*J)));
62   PetscFunctionReturn(0);
63 }
64 
65 PetscErrorCode MatSMFResetRowColumn(Mat mat, IS Rows, IS Cols) {
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(0);
77 }
78 
79 PetscErrorCode MatMult_SMF(Mat mat, Vec a, Vec y) {
80   MatSubMatFreeCtx ctx;
81 
82   PetscFunctionBegin;
83   PetscCall(MatShellGetContext(mat, &ctx));
84   PetscCall(VecCopy(a, ctx->VR));
85   PetscCall(VecISSet(ctx->VR, ctx->Cols, 0.0));
86   PetscCall(MatMult(ctx->A, ctx->VR, y));
87   PetscCall(VecISSet(y, ctx->Rows, 0.0));
88   PetscFunctionReturn(0);
89 }
90 
91 PetscErrorCode MatMultTranspose_SMF(Mat mat, Vec a, Vec y) {
92   MatSubMatFreeCtx ctx;
93 
94   PetscFunctionBegin;
95   PetscCall(MatShellGetContext(mat, &ctx));
96   PetscCall(VecCopy(a, ctx->VC));
97   PetscCall(VecISSet(ctx->VC, ctx->Rows, 0.0));
98   PetscCall(MatMultTranspose(ctx->A, ctx->VC, y));
99   PetscCall(VecISSet(y, ctx->Cols, 0.0));
100   PetscFunctionReturn(0);
101 }
102 
103 PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D, InsertMode is) {
104   MatSubMatFreeCtx ctx;
105 
106   PetscFunctionBegin;
107   PetscCall(MatShellGetContext(M, &ctx));
108   PetscCall(MatDiagonalSet(ctx->A, D, is));
109   PetscFunctionReturn(0);
110 }
111 
112 PetscErrorCode MatDestroy_SMF(Mat mat) {
113   MatSubMatFreeCtx ctx;
114 
115   PetscFunctionBegin;
116   PetscCall(MatShellGetContext(mat, &ctx));
117   PetscCall(MatDestroy(&ctx->A));
118   PetscCall(ISDestroy(&ctx->Rows));
119   PetscCall(ISDestroy(&ctx->Cols));
120   PetscCall(VecDestroy(&ctx->VC));
121   PetscCall(PetscFree(ctx));
122   PetscFunctionReturn(0);
123 }
124 
125 PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer) {
126   MatSubMatFreeCtx ctx;
127 
128   PetscFunctionBegin;
129   PetscCall(MatShellGetContext(mat, &ctx));
130   PetscCall(MatView(ctx->A, viewer));
131   PetscFunctionReturn(0);
132 }
133 
134 PetscErrorCode MatShift_SMF(Mat Y, PetscReal a) {
135   MatSubMatFreeCtx ctx;
136 
137   PetscFunctionBegin;
138   PetscCall(MatShellGetContext(Y, &ctx));
139   PetscCall(MatShift(ctx->A, a));
140   PetscFunctionReturn(0);
141 }
142 
143 PetscErrorCode MatDuplicate_SMF(Mat mat, MatDuplicateOption op, Mat *M) {
144   MatSubMatFreeCtx ctx;
145 
146   PetscFunctionBegin;
147   PetscCall(MatShellGetContext(mat, &ctx));
148   PetscCall(MatCreateSubMatrixFree(ctx->A, ctx->Rows, ctx->Cols, M));
149   PetscFunctionReturn(0);
150 }
151 
152 PetscErrorCode MatEqual_SMF(Mat A, Mat B, PetscBool *flg) {
153   MatSubMatFreeCtx ctx1, ctx2;
154   PetscBool        flg1, flg2, flg3;
155 
156   PetscFunctionBegin;
157   PetscCall(MatShellGetContext(A, &ctx1));
158   PetscCall(MatShellGetContext(B, &ctx2));
159   PetscCall(ISEqual(ctx1->Rows, ctx2->Rows, &flg2));
160   PetscCall(ISEqual(ctx1->Cols, ctx2->Cols, &flg3));
161   if (flg2 == PETSC_FALSE || flg3 == PETSC_FALSE) {
162     *flg = PETSC_FALSE;
163   } else {
164     PetscCall(MatEqual(ctx1->A, ctx2->A, &flg1));
165     if (flg1 == PETSC_FALSE) {
166       *flg = PETSC_FALSE;
167     } else {
168       *flg = PETSC_TRUE;
169     }
170   }
171   PetscFunctionReturn(0);
172 }
173 
174 PetscErrorCode MatScale_SMF(Mat mat, PetscReal a) {
175   MatSubMatFreeCtx ctx;
176 
177   PetscFunctionBegin;
178   PetscCall(MatShellGetContext(mat, &ctx));
179   PetscCall(MatScale(ctx->A, a));
180   PetscFunctionReturn(0);
181 }
182 
183 PetscErrorCode MatTranspose_SMF(Mat mat, Mat *B) {
184   PetscFunctionBegin;
185   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support for transpose for MatCreateSubMatrixFree() matrix");
186   PetscFunctionReturn(0);
187 }
188 
189 PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v) {
190   MatSubMatFreeCtx ctx;
191 
192   PetscFunctionBegin;
193   PetscCall(MatShellGetContext(mat, &ctx));
194   PetscCall(MatGetDiagonal(ctx->A, v));
195   PetscFunctionReturn(0);
196 }
197 
198 PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D) {
199   MatSubMatFreeCtx ctx;
200 
201   PetscFunctionBegin;
202   PetscCall(MatShellGetContext(M, &ctx));
203   PetscCall(MatGetRowMax(ctx->A, D, NULL));
204   PetscFunctionReturn(0);
205 }
206 
207 PetscErrorCode MatCreateSubMatrices_SMF(Mat A, PetscInt n, IS *irow, IS *icol, MatReuse scall, Mat **B) {
208   PetscInt i;
209 
210   PetscFunctionBegin;
211   if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));
212 
213   for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SMF(A, irow[i], icol[i], scall, &(*B)[i]));
214   PetscFunctionReturn(0);
215 }
216 
217 PetscErrorCode MatCreateSubMatrix_SMF(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat) {
218   MatSubMatFreeCtx ctx;
219 
220   PetscFunctionBegin;
221   PetscCall(MatShellGetContext(mat, &ctx));
222   if (newmat) PetscCall(MatDestroy(&*newmat));
223   PetscCall(MatCreateSubMatrixFree(ctx->A, isrow, iscol, newmat));
224   PetscFunctionReturn(0);
225 }
226 
227 PetscErrorCode MatGetRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals) {
228   MatSubMatFreeCtx ctx;
229 
230   PetscFunctionBegin;
231   PetscCall(MatShellGetContext(mat, &ctx));
232   PetscCall(MatGetRow(ctx->A, row, ncols, cols, vals));
233   PetscFunctionReturn(0);
234 }
235 
236 PetscErrorCode MatRestoreRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals) {
237   MatSubMatFreeCtx ctx;
238 
239   PetscFunctionBegin;
240   PetscCall(MatShellGetContext(mat, &ctx));
241   PetscCall(MatRestoreRow(ctx->A, row, ncols, cols, vals));
242   PetscFunctionReturn(0);
243 }
244 
245 PetscErrorCode MatGetColumnVector_SMF(Mat mat, Vec Y, PetscInt col) {
246   MatSubMatFreeCtx ctx;
247 
248   PetscFunctionBegin;
249   PetscCall(MatShellGetContext(mat, &ctx));
250   PetscCall(MatGetColumnVector(ctx->A, Y, col));
251   PetscFunctionReturn(0);
252 }
253 
254 PetscErrorCode MatNorm_SMF(Mat mat, NormType type, PetscReal *norm) {
255   MatSubMatFreeCtx ctx;
256 
257   PetscFunctionBegin;
258   PetscCall(MatShellGetContext(mat, &ctx));
259   if (type == NORM_FROBENIUS) {
260     *norm = 1.0;
261   } else if (type == NORM_1 || type == NORM_INFINITY) {
262     *norm = 1.0;
263   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
264   PetscFunctionReturn(0);
265 }
266