xref: /petsc/src/mat/impls/transpose/transm.c (revision bfe80ac4a46d58cb7760074b25f5e81b2f541d8a) !
1 #include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/
2 
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 
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 
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 
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 
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 
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 
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 
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 
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, (void (*)(void))MatSolve_Transpose_LU));
91   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_LU));
92   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_LU));
93   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_LU));
94   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_LU));
95   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_LU));
96   PetscFunctionReturn(PETSC_SUCCESS);
97 }
98 
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 
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 
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 
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 
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 
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 
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, (void (*)(void))MatSolve_Transpose_Cholesky));
167   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_Cholesky));
168   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_Cholesky));
169   PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_Cholesky));
170   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_Cholesky));
171   PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_Cholesky));
172   PetscFunctionReturn(PETSC_SUCCESS);
173 }
174 
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, (void (*)(void))MatSolve_Transpose_LU));
184   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_LU));
185   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_LU));
186   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_LU));
187   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_LU));
188   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_LU));
189   PetscFunctionReturn(PETSC_SUCCESS);
190 }
191 
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, (void (*)(void))MatLUFactorNumeric_Transpose));
201   PetscFunctionReturn(PETSC_SUCCESS);
202 }
203 
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, (void (*)(void))MatSolve_Transpose_Cholesky));
213   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_Cholesky));
214   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_Cholesky));
215   PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_Cholesky));
216   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_Cholesky));
217   PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_Cholesky));
218   PetscFunctionReturn(PETSC_SUCCESS);
219 }
220 
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, (void (*)(void))MatCholeskyFactorNumeric_Transpose));
230   PetscFunctionReturn(PETSC_SUCCESS);
231 }
232 
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, (void (*)(void))MatLUFactorSymbolic_Transpose));
242   else if (ftype == MAT_FACTOR_CHOLESKY) {
243     PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (void (*)(void))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 
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 
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 
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 
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 
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 static PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
314 {
315   Mat            A, B, C, Ain, Bin, Cin;
316   PetscBool      Aistrans, Bistrans, Cistrans;
317   PetscInt       Atrans, Btrans, Ctrans;
318   MatProductType ptype;
319 
320   PetscFunctionBegin;
321   MatCheckProduct(D, 1);
322   A = D->product->A;
323   B = D->product->B;
324   C = D->product->C;
325   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
326   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
327   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
328   PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
329   Atrans = 0;
330   Ain    = A;
331   while (Aistrans) {
332     Atrans++;
333     PetscCall(MatTransposeGetMat(Ain, &Ain));
334     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
335   }
336   Btrans = 0;
337   Bin    = B;
338   while (Bistrans) {
339     Btrans++;
340     PetscCall(MatTransposeGetMat(Bin, &Bin));
341     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
342   }
343   Ctrans = 0;
344   Cin    = C;
345   while (Cistrans) {
346     Ctrans++;
347     PetscCall(MatTransposeGetMat(Cin, &Cin));
348     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
349   }
350   Atrans = Atrans % 2;
351   Btrans = Btrans % 2;
352   Ctrans = Ctrans % 2;
353   ptype  = D->product->type; /* same product type by default */
354   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
355   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
356   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
357 
358   if (Atrans || Btrans || Ctrans) {
359     ptype = MATPRODUCT_UNSPECIFIED;
360     switch (D->product->type) {
361     case MATPRODUCT_AB:
362       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
363         /* TODO custom implementation ? */
364       } else if (Atrans) { /* At * B */
365         ptype = MATPRODUCT_AtB;
366       } else { /* A * Bt */
367         ptype = MATPRODUCT_ABt;
368       }
369       break;
370     case MATPRODUCT_AtB:
371       if (Atrans && Btrans) { /* A * Bt */
372         ptype = MATPRODUCT_ABt;
373       } else if (Atrans) { /* A * B */
374         ptype = MATPRODUCT_AB;
375       } else { /* At * Bt we do not have support for this */
376         /* TODO custom implementation ? */
377       }
378       break;
379     case MATPRODUCT_ABt:
380       if (Atrans && Btrans) { /* At * B */
381         ptype = MATPRODUCT_AtB;
382       } else if (Atrans) { /* At * Bt we do not have support for this */
383         /* TODO custom implementation ? */
384       } else { /* A * B */
385         ptype = MATPRODUCT_AB;
386       }
387       break;
388     case MATPRODUCT_PtAP:
389       if (Atrans) { /* PtAtP */
390         /* TODO custom implementation ? */
391       } else { /* RARt */
392         ptype = MATPRODUCT_RARt;
393       }
394       break;
395     case MATPRODUCT_RARt:
396       if (Atrans) { /* RAtRt */
397         /* TODO custom implementation ? */
398       } else { /* PtAP */
399         ptype = MATPRODUCT_PtAP;
400       }
401       break;
402     case MATPRODUCT_ABC:
403       /* TODO custom implementation ? */
404       break;
405     default:
406       SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
407     }
408   }
409   PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
410   PetscCall(MatProductSetType(D, ptype));
411   PetscCall(MatProductSetFromOptions(D));
412   PetscFunctionReturn(PETSC_SUCCESS);
413 }
414 
415 static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
416 {
417   Mat A;
418 
419   PetscFunctionBegin;
420   PetscCall(MatShellGetContext(N, &A));
421   PetscCall(MatGetDiagonal(A, v));
422   PetscFunctionReturn(PETSC_SUCCESS);
423 }
424 
425 static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
426 {
427   Mat a, b;
428 
429   PetscFunctionBegin;
430   PetscCall(MatShellGetContext(A, &a));
431   PetscCall(MatShellGetContext(B, &b));
432   PetscCall(MatCopy(a, b, str));
433   PetscFunctionReturn(PETSC_SUCCESS);
434 }
435 
436 static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
437 {
438   Mat         A;
439   PetscScalar vscale = 1.0, vshift = 0.0;
440   PetscBool   flg;
441 
442   PetscFunctionBegin;
443   PetscCall(MatShellGetContext(N, &A));
444   PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
445   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 */
446     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));
447   }
448   if (flg) {
449     Mat B;
450 
451     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
452     if (reuse != MAT_INPLACE_MATRIX) {
453       PetscCall(MatConvert(B, newtype, reuse, newmat));
454       PetscCall(MatDestroy(&B));
455     } else {
456       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
457       PetscCall(MatHeaderReplace(N, &B));
458     }
459   } else { /* use basic converter as fallback */
460     flg = (PetscBool)(N->ops->getrow != NULL);
461     PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
462   }
463   if (flg) {
464     PetscCall(MatScale(*newmat, vscale));
465     PetscCall(MatShift(*newmat, vshift));
466   }
467   PetscFunctionReturn(PETSC_SUCCESS);
468 }
469 
470 static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
471 {
472   PetscFunctionBegin;
473   PetscCall(MatShellGetContext(N, M));
474   PetscFunctionReturn(PETSC_SUCCESS);
475 }
476 
477 /*@
478   MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`
479 
480   Logically Collective
481 
482   Input Parameter:
483 . A - the `MATTRANSPOSEVIRTUAL` matrix
484 
485   Output Parameter:
486 . M - the matrix object stored inside `A`
487 
488   Level: intermediate
489 
490 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
491 @*/
492 PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
493 {
494   PetscFunctionBegin;
495   PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
496   PetscValidType(A, 1);
497   PetscAssertPointer(M, 2);
498   PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
499   PetscFunctionReturn(PETSC_SUCCESS);
500 }
501 
502 /*MC
503    MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix
504 
505   Level: advanced
506 
507   Developer Notes:
508   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
509 
510   Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
511 
512 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
513           `MATNORMALHERMITIAN`, `MATNORMAL`
514 M*/
515 
516 /*@
517   MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'
518 
519   Collective
520 
521   Input Parameter:
522 . A - the (possibly rectangular) matrix
523 
524   Output Parameter:
525 . N - the matrix that represents A'
526 
527   Level: intermediate
528 
529   Note:
530   The transpose A' is NOT actually formed! Rather the new matrix
531   object performs the matrix-vector product by using the `MatMultTranspose()` on
532   the original matrix
533 
534 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
535           `MATNORMALHERMITIAN`
536 @*/
537 PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
538 {
539   VecType vtype;
540 
541   PetscFunctionBegin;
542   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
543   PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
544   PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
545   PetscCall(MatSetType(*N, MATSHELL));
546   PetscCall(MatShellSetContext(*N, A));
547   PetscCall(PetscObjectReference((PetscObject)A));
548 
549   PetscCall(MatSetBlockSizes(*N, A->cmap->bs, A->rmap->bs));
550   PetscCall(MatGetVecType(A, &vtype));
551   PetscCall(MatSetVecType(*N, vtype));
552 #if defined(PETSC_HAVE_DEVICE)
553   PetscCall(MatBindToCPU(*N, A->boundtocpu));
554 #endif
555   PetscCall(MatSetUp(*N));
556 
557   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Transpose));
558   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Transpose));
559   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Transpose));
560   PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (void (*)(void))MatLUFactor_Transpose));
561   PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (void (*)(void))MatCholeskyFactor_Transpose));
562   PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (void (*)(void))MatGetFactor_Transpose));
563   PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (void (*)(void))MatGetInfo_Transpose));
564   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Transpose));
565   PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_Transpose));
566   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Transpose));
567   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Transpose));
568   PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_Transpose));
569 
570   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
571   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
572   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_Transpose));
573   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
574   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
575   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
576   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
577   PetscFunctionReturn(PETSC_SUCCESS);
578 }
579