Actual source code: transm.c
1: #include <../src/mat/impls/shell/shell.h>
3: static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y)
4: {
5: Mat A;
7: PetscFunctionBegin;
8: PetscCall(MatShellGetContext(N, &A));
9: PetscCall(MatMultTranspose(A, x, y));
10: PetscFunctionReturn(PETSC_SUCCESS);
11: }
13: static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y)
14: {
15: Mat A;
17: PetscFunctionBegin;
18: PetscCall(MatShellGetContext(N, &A));
19: PetscCall(MatMult(A, x, y));
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: static PetscErrorCode MatDestroy_Transpose(Mat N)
24: {
25: Mat A;
27: PetscFunctionBegin;
28: PetscCall(MatShellGetContext(N, &A));
29: PetscCall(MatDestroy(&A));
30: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
31: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
32: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
37: {
38: Mat A, C;
40: PetscFunctionBegin;
41: PetscCall(MatShellGetContext(N, &A));
42: PetscCall(MatDuplicate(A, op, &C));
43: PetscCall(MatCreateTranspose(C, m));
44: PetscCall(MatDestroy(&C));
45: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
50: {
51: Mat A;
53: PetscFunctionBegin;
54: PetscCall(MatShellGetContext(mat, &A));
55: *has = PETSC_FALSE;
56: if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
57: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
58: } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
59: PetscCall(MatHasOperation(A, MATOP_MULT, has));
60: } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
65: {
66: Mat A, B, C, Ain, Bin, Cin;
67: PetscBool Aistrans, Bistrans, Cistrans;
68: PetscInt Atrans, Btrans, Ctrans;
69: MatProductType ptype;
71: PetscFunctionBegin;
72: MatCheckProduct(D, 1);
73: A = D->product->A;
74: B = D->product->B;
75: C = D->product->C;
76: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
77: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
78: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
79: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
80: Atrans = 0;
81: Ain = A;
82: while (Aistrans) {
83: Atrans++;
84: PetscCall(MatTransposeGetMat(Ain, &Ain));
85: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
86: }
87: Btrans = 0;
88: Bin = B;
89: while (Bistrans) {
90: Btrans++;
91: PetscCall(MatTransposeGetMat(Bin, &Bin));
92: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
93: }
94: Ctrans = 0;
95: Cin = C;
96: while (Cistrans) {
97: Ctrans++;
98: PetscCall(MatTransposeGetMat(Cin, &Cin));
99: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
100: }
101: Atrans = Atrans % 2;
102: Btrans = Btrans % 2;
103: Ctrans = Ctrans % 2;
104: ptype = D->product->type; /* same product type by default */
105: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
106: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
107: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
109: if (Atrans || Btrans || Ctrans) {
110: ptype = MATPRODUCT_UNSPECIFIED;
111: switch (D->product->type) {
112: case MATPRODUCT_AB:
113: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
114: /* TODO custom implementation ? */
115: } else if (Atrans) { /* At * B */
116: ptype = MATPRODUCT_AtB;
117: } else { /* A * Bt */
118: ptype = MATPRODUCT_ABt;
119: }
120: break;
121: case MATPRODUCT_AtB:
122: if (Atrans && Btrans) { /* A * Bt */
123: ptype = MATPRODUCT_ABt;
124: } else if (Atrans) { /* A * B */
125: ptype = MATPRODUCT_AB;
126: } else { /* At * Bt we do not have support for this */
127: /* TODO custom implementation ? */
128: }
129: break;
130: case MATPRODUCT_ABt:
131: if (Atrans && Btrans) { /* At * B */
132: ptype = MATPRODUCT_AtB;
133: } else if (Atrans) { /* At * Bt we do not have support for this */
134: /* TODO custom implementation ? */
135: } else { /* A * B */
136: ptype = MATPRODUCT_AB;
137: }
138: break;
139: case MATPRODUCT_PtAP:
140: if (Atrans) { /* PtAtP */
141: /* TODO custom implementation ? */
142: } else { /* RARt */
143: ptype = MATPRODUCT_RARt;
144: }
145: break;
146: case MATPRODUCT_RARt:
147: if (Atrans) { /* RAtRt */
148: /* TODO custom implementation ? */
149: } else { /* PtAP */
150: ptype = MATPRODUCT_PtAP;
151: }
152: break;
153: case MATPRODUCT_ABC:
154: /* TODO custom implementation ? */
155: break;
156: default:
157: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
158: }
159: }
160: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
161: PetscCall(MatProductSetType(D, ptype));
162: PetscCall(MatProductSetFromOptions(D));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
167: {
168: Mat A;
170: PetscFunctionBegin;
171: PetscCall(MatShellGetContext(N, &A));
172: PetscCall(MatGetDiagonal(A, v));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
177: {
178: Mat a, b;
180: PetscFunctionBegin;
181: PetscCall(MatShellGetContext(A, &a));
182: PetscCall(MatShellGetContext(B, &b));
183: PetscCall(MatCopy(a, b, str));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
188: {
189: Mat A;
190: PetscScalar vscale = 1.0, vshift = 0.0;
191: PetscBool flg;
193: PetscFunctionBegin;
194: PetscCall(MatShellGetContext(N, &A));
195: PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
196: 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 */
197: PetscCheck(!((Mat_Shell *)N->data)->zrows && !((Mat_Shell *)N->data)->zcols, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat");
198: PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatAXPY() has been called on the input Mat");
199: PetscCheck(!((Mat_Shell *)N->data)->left && !((Mat_Shell *)N->data)->right, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalScale() has been called on the input Mat");
200: PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalSet() has been called on the input Mat");
201: vscale = ((Mat_Shell *)N->data)->vscale;
202: vshift = ((Mat_Shell *)N->data)->vshift;
203: }
204: if (flg) {
205: Mat B;
207: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
208: if (reuse != MAT_INPLACE_MATRIX) {
209: PetscCall(MatConvert(B, newtype, reuse, newmat));
210: PetscCall(MatDestroy(&B));
211: } else {
212: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
213: PetscCall(MatHeaderReplace(N, &B));
214: }
215: } else { /* use basic converter as fallback */
216: flg = (PetscBool)(N->ops->getrow != NULL);
217: PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
218: }
219: if (flg) {
220: PetscCall(MatScale(*newmat, vscale));
221: PetscCall(MatShift(*newmat, vshift));
222: }
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
227: {
228: PetscFunctionBegin;
229: PetscCall(MatShellGetContext(N, M));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*@
234: MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`
236: Logically Collective
238: Input Parameter:
239: . A - the `MATTRANSPOSEVIRTUAL` matrix
241: Output Parameter:
242: . M - the matrix object stored inside `A`
244: Level: intermediate
246: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
247: @*/
248: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
249: {
250: PetscFunctionBegin;
253: PetscAssertPointer(M, 2);
254: PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*MC
259: MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix
261: Level: advanced
263: Developer Notes:
264: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
266: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
268: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
269: `MATNORMALHERMITIAN`, `MATNORMAL`
270: M*/
272: /*@
273: MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'
275: Collective
277: Input Parameter:
278: . A - the (possibly rectangular) matrix
280: Output Parameter:
281: . N - the matrix that represents A'
283: Level: intermediate
285: Note:
286: The transpose A' is NOT actually formed! Rather the new matrix
287: object performs the matrix-vector product by using the `MatMultTranspose()` on
288: the original matrix
290: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
291: `MATNORMALHERMITIAN`
292: @*/
293: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
294: {
295: VecType vtype;
297: PetscFunctionBegin;
298: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
299: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
300: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
301: PetscCall(MatSetType(*N, MATSHELL));
302: PetscCall(MatShellSetContext(*N, A));
303: PetscCall(PetscObjectReference((PetscObject)A));
305: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
306: PetscCall(MatGetVecType(A, &vtype));
307: PetscCall(MatSetVecType(*N, vtype));
308: #if defined(PETSC_HAVE_DEVICE)
309: PetscCall(MatBindToCPU(*N, A->boundtocpu));
310: #endif
311: PetscCall(MatSetUp(*N));
313: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Transpose));
314: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Transpose));
315: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Transpose));
316: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Transpose));
317: PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_Transpose));
318: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Transpose));
319: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Transpose));
320: PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_Transpose));
322: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
323: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
324: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
325: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
326: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
327: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }