1 #include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/
2
MatMult_Transpose(Mat N,Vec x,Vec y)3 static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y)
4 {
5 Mat A;
6
7 PetscFunctionBegin;
8 PetscCall(MatShellGetContext(N, &A));
9 PetscCall(MatMultTranspose(A, x, y));
10 PetscFunctionReturn(PETSC_SUCCESS);
11 }
12
MatMultTranspose_Transpose(Mat N,Vec x,Vec y)13 static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y)
14 {
15 Mat A;
16
17 PetscFunctionBegin;
18 PetscCall(MatShellGetContext(N, &A));
19 PetscCall(MatMult(A, x, y));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
MatSolve_Transpose_LU(Mat N,Vec b,Vec x)23 static PetscErrorCode MatSolve_Transpose_LU(Mat N, Vec b, Vec x)
24 {
25 Mat A;
26
27 PetscFunctionBegin;
28 PetscCall(MatShellGetContext(N, &A));
29 PetscCall(MatSolveTranspose(A, b, x));
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
MatSolveAdd_Transpose_LU(Mat N,Vec b,Vec y,Vec x)33 static PetscErrorCode MatSolveAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x)
34 {
35 Mat A;
36
37 PetscFunctionBegin;
38 PetscCall(MatShellGetContext(N, &A));
39 PetscCall(MatSolveTransposeAdd(A, b, y, x));
40 PetscFunctionReturn(PETSC_SUCCESS);
41 }
42
MatSolveTranspose_Transpose_LU(Mat N,Vec b,Vec x)43 static PetscErrorCode MatSolveTranspose_Transpose_LU(Mat N, Vec b, Vec x)
44 {
45 Mat A;
46
47 PetscFunctionBegin;
48 PetscCall(MatShellGetContext(N, &A));
49 PetscCall(MatSolve(A, b, x));
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52
MatSolveTransposeAdd_Transpose_LU(Mat N,Vec b,Vec y,Vec x)53 static PetscErrorCode MatSolveTransposeAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x)
54 {
55 Mat A;
56
57 PetscFunctionBegin;
58 PetscCall(MatShellGetContext(N, &A));
59 PetscCall(MatSolveAdd(A, b, y, x));
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
MatMatSolve_Transpose_LU(Mat N,Mat B,Mat X)63 static PetscErrorCode MatMatSolve_Transpose_LU(Mat N, Mat B, Mat X)
64 {
65 Mat A;
66
67 PetscFunctionBegin;
68 PetscCall(MatShellGetContext(N, &A));
69 PetscCall(MatMatSolveTranspose(A, B, X));
70 PetscFunctionReturn(PETSC_SUCCESS);
71 }
72
MatMatSolveTranspose_Transpose_LU(Mat N,Mat B,Mat X)73 static PetscErrorCode MatMatSolveTranspose_Transpose_LU(Mat N, Mat B, Mat X)
74 {
75 Mat A;
76
77 PetscFunctionBegin;
78 PetscCall(MatShellGetContext(N, &A));
79 PetscCall(MatMatSolve(A, B, X));
80 PetscFunctionReturn(PETSC_SUCCESS);
81 }
82
MatLUFactor_Transpose(Mat N,IS row,IS col,const MatFactorInfo * minfo)83 static PetscErrorCode MatLUFactor_Transpose(Mat N, IS row, IS col, const MatFactorInfo *minfo)
84 {
85 Mat A;
86
87 PetscFunctionBegin;
88 PetscCall(MatShellGetContext(N, &A));
89 PetscCall(MatLUFactor(A, col, row, minfo));
90 PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_LU));
91 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_LU));
92 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_LU));
93 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_LU));
94 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_LU));
95 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_LU));
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
MatSolve_Transpose_Cholesky(Mat N,Vec b,Vec x)99 static PetscErrorCode MatSolve_Transpose_Cholesky(Mat N, Vec b, Vec x)
100 {
101 Mat A;
102
103 PetscFunctionBegin;
104 PetscCall(MatShellGetContext(N, &A));
105 PetscCall(MatSolveTranspose(A, b, x));
106 PetscFunctionReturn(PETSC_SUCCESS);
107 }
108
MatSolveAdd_Transpose_Cholesky(Mat N,Vec b,Vec y,Vec x)109 static PetscErrorCode MatSolveAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x)
110 {
111 Mat A;
112
113 PetscFunctionBegin;
114 PetscCall(MatShellGetContext(N, &A));
115 PetscCall(MatSolveTransposeAdd(A, b, y, x));
116 PetscFunctionReturn(PETSC_SUCCESS);
117 }
118
MatSolveTranspose_Transpose_Cholesky(Mat N,Vec b,Vec x)119 static PetscErrorCode MatSolveTranspose_Transpose_Cholesky(Mat N, Vec b, Vec x)
120 {
121 Mat A;
122
123 PetscFunctionBegin;
124 PetscCall(MatShellGetContext(N, &A));
125 PetscCall(MatSolve(A, b, x));
126 PetscFunctionReturn(PETSC_SUCCESS);
127 }
128
MatSolveTransposeAdd_Transpose_Cholesky(Mat N,Vec b,Vec y,Vec x)129 static PetscErrorCode MatSolveTransposeAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x)
130 {
131 Mat A;
132
133 PetscFunctionBegin;
134 PetscCall(MatShellGetContext(N, &A));
135 PetscCall(MatSolveAdd(A, b, y, x));
136 PetscFunctionReturn(PETSC_SUCCESS);
137 }
138
MatMatSolve_Transpose_Cholesky(Mat N,Mat B,Mat X)139 static PetscErrorCode MatMatSolve_Transpose_Cholesky(Mat N, Mat B, Mat X)
140 {
141 Mat A;
142
143 PetscFunctionBegin;
144 PetscCall(MatShellGetContext(N, &A));
145 PetscCall(MatMatSolveTranspose(A, B, X));
146 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
MatMatSolveTranspose_Transpose_Cholesky(Mat N,Mat B,Mat X)149 static PetscErrorCode MatMatSolveTranspose_Transpose_Cholesky(Mat N, Mat B, Mat X)
150 {
151 Mat A;
152
153 PetscFunctionBegin;
154 PetscCall(MatShellGetContext(N, &A));
155 PetscCall(MatMatSolve(A, B, X));
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
MatCholeskyFactor_Transpose(Mat N,IS perm,const MatFactorInfo * minfo)159 static PetscErrorCode MatCholeskyFactor_Transpose(Mat N, IS perm, const MatFactorInfo *minfo)
160 {
161 Mat A;
162
163 PetscFunctionBegin;
164 PetscCall(MatShellGetContext(N, &A));
165 PetscCall(MatCholeskyFactor(A, perm, minfo));
166 PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_Cholesky));
167 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_Cholesky));
168 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_Cholesky));
169 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_Cholesky));
170 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_Cholesky));
171 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_Cholesky));
172 PetscFunctionReturn(PETSC_SUCCESS);
173 }
174
MatLUFactorNumeric_Transpose(Mat F,Mat N,const MatFactorInfo * info)175 static PetscErrorCode MatLUFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info)
176 {
177 Mat A, FA;
178
179 PetscFunctionBegin;
180 PetscCall(MatShellGetContext(N, &A));
181 PetscCall(MatShellGetContext(F, &FA));
182 PetscCall(MatLUFactorNumeric(FA, A, info));
183 PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_LU));
184 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_LU));
185 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_LU));
186 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_LU));
187 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_LU));
188 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_LU));
189 PetscFunctionReturn(PETSC_SUCCESS);
190 }
191
MatLUFactorSymbolic_Transpose(Mat F,Mat N,IS row,IS col,const MatFactorInfo * info)192 static PetscErrorCode MatLUFactorSymbolic_Transpose(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info)
193 {
194 Mat A, FA;
195
196 PetscFunctionBegin;
197 PetscCall(MatShellGetContext(N, &A));
198 PetscCall(MatShellGetContext(F, &FA));
199 PetscCall(MatLUFactorSymbolic(FA, A, row, col, info));
200 PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (PetscErrorCodeFn *)MatLUFactorNumeric_Transpose));
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
MatCholeskyFactorNumeric_Transpose(Mat F,Mat N,const MatFactorInfo * info)204 static PetscErrorCode MatCholeskyFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info)
205 {
206 Mat A, FA;
207
208 PetscFunctionBegin;
209 PetscCall(MatShellGetContext(N, &A));
210 PetscCall(MatShellGetContext(F, &FA));
211 PetscCall(MatCholeskyFactorNumeric(FA, A, info));
212 PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_Cholesky));
213 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_Cholesky));
214 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_Cholesky));
215 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_Cholesky));
216 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_Cholesky));
217 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_Cholesky));
218 PetscFunctionReturn(PETSC_SUCCESS);
219 }
220
MatCholeskyFactorSymbolic_Transpose(Mat F,Mat N,IS perm,const MatFactorInfo * info)221 static PetscErrorCode MatCholeskyFactorSymbolic_Transpose(Mat F, Mat N, IS perm, const MatFactorInfo *info)
222 {
223 Mat A, FA;
224
225 PetscFunctionBegin;
226 PetscCall(MatShellGetContext(N, &A));
227 PetscCall(MatShellGetContext(F, &FA));
228 PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info));
229 PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (PetscErrorCodeFn *)MatCholeskyFactorNumeric_Transpose));
230 PetscFunctionReturn(PETSC_SUCCESS);
231 }
232
MatGetFactor_Transpose(Mat N,MatSolverType type,MatFactorType ftype,Mat * F)233 static PetscErrorCode MatGetFactor_Transpose(Mat N, MatSolverType type, MatFactorType ftype, Mat *F)
234 {
235 Mat A, FA;
236
237 PetscFunctionBegin;
238 PetscCall(MatShellGetContext(N, &A));
239 PetscCall(MatGetFactor(A, type, ftype, &FA));
240 PetscCall(MatCreateTranspose(FA, F));
241 if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatLUFactorSymbolic_Transpose));
242 else if (ftype == MAT_FACTOR_CHOLESKY) {
243 PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatCholeskyFactorSymbolic_Transpose));
244 PetscCall(MatPropagateSymmetryOptions(A, FA));
245 } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]);
246 (*F)->factortype = ftype;
247 PetscCall(MatDestroy(&FA));
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
MatDestroy_Transpose(Mat N)251 static PetscErrorCode MatDestroy_Transpose(Mat N)
252 {
253 Mat A;
254
255 PetscFunctionBegin;
256 PetscCall(MatShellGetContext(N, &A));
257 PetscCall(MatDestroy(&A));
258 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
259 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
260 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
261 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL));
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
MatGetInfo_Transpose(Mat N,MatInfoType flag,MatInfo * info)265 static PetscErrorCode MatGetInfo_Transpose(Mat N, MatInfoType flag, MatInfo *info)
266 {
267 Mat A;
268
269 PetscFunctionBegin;
270 PetscCall(MatShellGetContext(N, &A));
271 PetscCall(MatGetInfo(A, flag, info));
272 PetscFunctionReturn(PETSC_SUCCESS);
273 }
274
MatFactorGetSolverType_Transpose(Mat N,MatSolverType * type)275 static PetscErrorCode MatFactorGetSolverType_Transpose(Mat N, MatSolverType *type)
276 {
277 Mat A;
278
279 PetscFunctionBegin;
280 PetscCall(MatShellGetContext(N, &A));
281 PetscCall(MatFactorGetSolverType(A, type));
282 PetscFunctionReturn(PETSC_SUCCESS);
283 }
284
MatDuplicate_Transpose(Mat N,MatDuplicateOption op,Mat * m)285 static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
286 {
287 Mat A, C;
288
289 PetscFunctionBegin;
290 PetscCall(MatShellGetContext(N, &A));
291 PetscCall(MatDuplicate(A, op, &C));
292 PetscCall(MatCreateTranspose(C, m));
293 if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
294 PetscCall(MatDestroy(&C));
295 PetscFunctionReturn(PETSC_SUCCESS);
296 }
297
MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool * has)298 static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
299 {
300 Mat A;
301
302 PetscFunctionBegin;
303 PetscCall(MatShellGetContext(mat, &A));
304 *has = PETSC_FALSE;
305 if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
306 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
307 } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
308 PetscCall(MatHasOperation(A, MATOP_MULT, has));
309 } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
310 PetscFunctionReturn(PETSC_SUCCESS);
311 }
312
313 typedef struct {
314 PetscErrorCode (*numeric)(Mat);
315 PetscCtxDestroyFn *destroy;
316 PetscScalar scale;
317 Mat D;
318 void *data;
319 } MatProductCtx_Transpose;
320
MatProductCtxDestroy_Transpose(PetscCtxRt ptr)321 static PetscErrorCode MatProductCtxDestroy_Transpose(PetscCtxRt ptr)
322 {
323 MatProductCtx_Transpose *data = *(MatProductCtx_Transpose **)ptr;
324 PetscContainer container;
325
326 PetscFunctionBegin;
327 if (data->data) PetscCall((*data->destroy)(&data->data));
328 PetscCall(PetscObjectQuery((PetscObject)data->D, "MatProductCtx_Transpose", (PetscObject *)&container));
329 PetscCall(PetscContainerDestroy(&container));
330 PetscCall(PetscObjectCompose((PetscObject)data->D, "MatProductCtx_Transpose", NULL));
331 PetscCall(PetscFree(data));
332 PetscFunctionReturn(PETSC_SUCCESS);
333 }
334
MatProductNumeric_Transpose(Mat D)335 static PetscErrorCode MatProductNumeric_Transpose(Mat D)
336 {
337 Mat_Product *product;
338 MatProductCtx_Transpose *data;
339 PetscContainer container;
340
341 PetscFunctionBegin;
342 MatCheckProduct(D, 1);
343 PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
344 product = D->product;
345 PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
346 PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_Transpose missing");
347 PetscCall(PetscContainerGetPointer(container, &data));
348 data = (MatProductCtx_Transpose *)product->data;
349 product->data = data->data;
350 PetscCall((*data->numeric)(D));
351 PetscCall(MatScale(D, data->scale));
352 product->data = data;
353 PetscFunctionReturn(PETSC_SUCCESS);
354 }
355
MatProductSymbolic_Transpose(Mat D)356 static PetscErrorCode MatProductSymbolic_Transpose(Mat D)
357 {
358 Mat_Product *product;
359 MatProductCtx_Transpose *data;
360 PetscContainer container;
361
362 PetscFunctionBegin;
363 MatCheckProduct(D, 1);
364 product = D->product;
365 if (D->ops->productsymbolic == MatProductSymbolic_Transpose) {
366 PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
367 PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
368 PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_Transpose missing");
369 PetscCall(PetscContainerGetPointer(container, &data));
370 PetscCall(MatProductSetFromOptions(D));
371 PetscCall(MatProductSymbolic(D));
372 data->numeric = D->ops->productnumeric;
373 data->destroy = product->destroy;
374 data->data = product->data;
375 D->ops->productnumeric = MatProductNumeric_Transpose;
376 product->destroy = MatProductCtxDestroy_Transpose;
377 product->data = data;
378 }
379 PetscFunctionReturn(PETSC_SUCCESS);
380 }
381
MatProductSetFromOptions_Transpose(Mat D)382 static PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
383 {
384 Mat A, B, C, Ain, Bin, Cin;
385 PetscScalar scale = 1.0, vscale;
386 PetscBool Aistrans, Bistrans, Cistrans;
387 PetscInt Atrans, Btrans, Ctrans;
388 PetscContainer container = NULL;
389 MatProductCtx_Transpose *data;
390 MatProductType ptype;
391
392 PetscFunctionBegin;
393 MatCheckProduct(D, 1);
394 A = D->product->A;
395 B = D->product->B;
396 C = D->product->C;
397 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
398 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
399 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
400 PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
401 Atrans = 0;
402 Ain = A;
403 while (Aistrans) {
404 Atrans++;
405 PetscCall(MatShellGetScalingShifts(Ain, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
406 scale *= vscale;
407 PetscCall(MatTransposeGetMat(Ain, &Ain));
408 PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
409 }
410 Btrans = 0;
411 Bin = B;
412 while (Bistrans) {
413 Btrans++;
414 PetscCall(MatShellGetScalingShifts(Bin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
415 scale *= vscale;
416 PetscCall(MatTransposeGetMat(Bin, &Bin));
417 PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
418 }
419 Ctrans = 0;
420 Cin = C;
421 while (Cistrans) {
422 Ctrans++;
423 PetscCall(MatShellGetScalingShifts(Cin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
424 scale *= vscale;
425 PetscCall(MatTransposeGetMat(Cin, &Cin));
426 PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
427 }
428 Atrans = Atrans % 2;
429 Btrans = Btrans % 2;
430 Ctrans = Ctrans % 2;
431 ptype = D->product->type; /* same product type by default */
432 if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
433 if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
434 if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
435
436 if (Atrans || Btrans || Ctrans) {
437 if (scale != 1.0) {
438 PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
439 if (!container) {
440 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)D), &container));
441 PetscCall(PetscNew(&data));
442 PetscCall(PetscContainerSetPointer(container, data));
443 PetscCall(PetscObjectCompose((PetscObject)D, "MatProductCtx_Transpose", (PetscObject)container));
444 } else PetscCall(PetscContainerGetPointer(container, &data));
445 data->scale = scale;
446 data->D = D;
447 }
448 ptype = MATPRODUCT_UNSPECIFIED;
449 switch (D->product->type) {
450 case MATPRODUCT_AB:
451 if (Atrans && Btrans) { /* At * Bt we do not have support for this */
452 /* TODO custom implementation ? */
453 } else if (Atrans) { /* At * B */
454 ptype = MATPRODUCT_AtB;
455 } else { /* A * Bt */
456 ptype = MATPRODUCT_ABt;
457 }
458 break;
459 case MATPRODUCT_AtB:
460 if (Atrans && Btrans) { /* A * Bt */
461 ptype = MATPRODUCT_ABt;
462 } else if (Atrans) { /* A * B */
463 ptype = MATPRODUCT_AB;
464 } else { /* At * Bt we do not have support for this */
465 /* TODO custom implementation ? */
466 }
467 break;
468 case MATPRODUCT_ABt:
469 if (Atrans && Btrans) { /* At * B */
470 ptype = MATPRODUCT_AtB;
471 } else if (Atrans) { /* At * Bt we do not have support for this */
472 /* TODO custom implementation ? */
473 } else { /* A * B */
474 ptype = MATPRODUCT_AB;
475 }
476 break;
477 case MATPRODUCT_PtAP:
478 if (Atrans) { /* PtAtP */
479 /* TODO custom implementation ? */
480 } else { /* RARt */
481 ptype = MATPRODUCT_RARt;
482 }
483 break;
484 case MATPRODUCT_RARt:
485 if (Atrans) { /* RAtRt */
486 /* TODO custom implementation ? */
487 } else { /* PtAP */
488 ptype = MATPRODUCT_PtAP;
489 }
490 break;
491 case MATPRODUCT_ABC:
492 /* TODO custom implementation ? */
493 break;
494 default:
495 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
496 }
497 }
498 PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
499 PetscCall(MatProductSetType(D, ptype));
500 if (container) D->ops->productsymbolic = MatProductSymbolic_Transpose;
501 else PetscCall(MatProductSetFromOptions(D));
502 PetscFunctionReturn(PETSC_SUCCESS);
503 }
504
MatGetDiagonal_Transpose(Mat N,Vec v)505 static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
506 {
507 Mat A;
508
509 PetscFunctionBegin;
510 PetscCall(MatShellGetContext(N, &A));
511 PetscCall(MatGetDiagonal(A, v));
512 PetscFunctionReturn(PETSC_SUCCESS);
513 }
514
MatCopy_Transpose(Mat A,Mat B,MatStructure str)515 static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
516 {
517 Mat a, b;
518
519 PetscFunctionBegin;
520 PetscCall(MatShellGetContext(A, &a));
521 PetscCall(MatShellGetContext(B, &b));
522 PetscCall(MatCopy(a, b, str));
523 PetscFunctionReturn(PETSC_SUCCESS);
524 }
525
MatConvert_Transpose(Mat N,MatType newtype,MatReuse reuse,Mat * newmat)526 static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
527 {
528 Mat A;
529 PetscScalar vscale = 1.0, vshift = 0.0;
530 PetscBool flg;
531
532 PetscFunctionBegin;
533 PetscCall(MatShellGetContext(N, &A));
534 PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
535 if (flg || N->ops->getrow) { /* if this condition is false, MatConvert_Shell() will be called in MatConvert_Basic(), so the following checks are not needed */
536 PetscCall(MatShellGetScalingShifts(N, &vshift, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
537 }
538 if (flg) {
539 Mat B;
540
541 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
542 if (reuse != MAT_INPLACE_MATRIX) {
543 PetscCall(MatConvert(B, newtype, reuse, newmat));
544 PetscCall(MatDestroy(&B));
545 } else {
546 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
547 PetscCall(MatHeaderReplace(N, &B));
548 }
549 } else { /* use basic converter as fallback */
550 flg = (PetscBool)(N->ops->getrow != NULL);
551 PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
552 }
553 if (flg) {
554 PetscCall(MatScale(*newmat, vscale));
555 PetscCall(MatShift(*newmat, vshift));
556 }
557 PetscFunctionReturn(PETSC_SUCCESS);
558 }
559
MatTransposeGetMat_Transpose(Mat N,Mat * M)560 static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
561 {
562 PetscFunctionBegin;
563 PetscCall(MatShellGetContext(N, M));
564 PetscFunctionReturn(PETSC_SUCCESS);
565 }
566
567 /*@
568 MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`
569
570 Logically Collective
571
572 Input Parameter:
573 . A - the `MATTRANSPOSEVIRTUAL` matrix
574
575 Output Parameter:
576 . M - the matrix object stored inside `A`
577
578 Level: intermediate
579
580 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
581 @*/
MatTransposeGetMat(Mat A,Mat * M)582 PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
583 {
584 PetscFunctionBegin;
585 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
586 PetscValidType(A, 1);
587 PetscAssertPointer(M, 2);
588 PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
589 PetscFunctionReturn(PETSC_SUCCESS);
590 }
591
592 /*MC
593 MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix
594
595 Level: advanced
596
597 Developer Notes:
598 This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
599
600 Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
601
602 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
603 `MATNORMALHERMITIAN`, `MATNORMAL`
604 M*/
605
606 /*@
607 MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'
608
609 Collective
610
611 Input Parameter:
612 . A - the (possibly rectangular) matrix
613
614 Output Parameter:
615 . N - the matrix that represents A'
616
617 Level: intermediate
618
619 Note:
620 The transpose A' is NOT actually formed! Rather the new matrix
621 object performs the matrix-vector product by using the `MatMultTranspose()` on
622 the original matrix
623
624 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
625 `MATNORMALHERMITIAN`
626 @*/
MatCreateTranspose(Mat A,Mat * N)627 PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
628 {
629 VecType vtype;
630
631 PetscFunctionBegin;
632 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
633 PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
634 PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
635 PetscCall(MatSetType(*N, MATSHELL));
636 PetscCall(MatShellSetContext(*N, A));
637 PetscCall(PetscObjectReference((PetscObject)A));
638
639 PetscCall(MatSetBlockSizes(*N, A->cmap->bs, A->rmap->bs));
640 PetscCall(MatGetVecType(A, &vtype));
641 PetscCall(MatSetVecType(*N, vtype));
642 #if defined(PETSC_HAVE_DEVICE)
643 PetscCall(MatBindToCPU(*N, A->boundtocpu));
644 #endif
645 PetscCall(MatSetUp(*N));
646
647 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Transpose));
648 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Transpose));
649 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Transpose));
650 PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (PetscErrorCodeFn *)MatLUFactor_Transpose));
651 PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (PetscErrorCodeFn *)MatCholeskyFactor_Transpose));
652 PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (PetscErrorCodeFn *)MatGetFactor_Transpose));
653 PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (PetscErrorCodeFn *)MatGetInfo_Transpose));
654 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_Transpose));
655 PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (PetscErrorCodeFn *)MatHasOperation_Transpose));
656 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Transpose));
657 PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_Transpose));
658 PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (PetscErrorCodeFn *)MatConvert_Transpose));
659
660 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
661 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
662 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_Transpose));
663 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
664 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
665 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
666 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
667 PetscFunctionReturn(PETSC_SUCCESS);
668 }
669