Actual source code: transm.c
1: #include <petsc/private/matimpl.h>
3: typedef struct {
4: Mat A;
5: } Mat_Transpose;
7: static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y)
8: {
9: Mat_Transpose *Na = (Mat_Transpose *)N->data;
11: PetscFunctionBegin;
12: PetscCall(MatMultTranspose(Na->A, x, y));
13: PetscFunctionReturn(PETSC_SUCCESS);
14: }
16: static PetscErrorCode MatMultAdd_Transpose(Mat N, Vec v1, Vec v2, Vec v3)
17: {
18: Mat_Transpose *Na = (Mat_Transpose *)N->data;
20: PetscFunctionBegin;
21: PetscCall(MatMultTransposeAdd(Na->A, v1, v2, v3));
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y)
26: {
27: Mat_Transpose *Na = (Mat_Transpose *)N->data;
29: PetscFunctionBegin;
30: PetscCall(MatMult(Na->A, x, y));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode MatMultTransposeAdd_Transpose(Mat N, Vec v1, Vec v2, Vec v3)
35: {
36: Mat_Transpose *Na = (Mat_Transpose *)N->data;
38: PetscFunctionBegin;
39: PetscCall(MatMultAdd(Na->A, v1, v2, v3));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode MatDestroy_Transpose(Mat N)
44: {
45: Mat_Transpose *Na = (Mat_Transpose *)N->data;
47: PetscFunctionBegin;
48: PetscCall(MatDestroy(&Na->A));
49: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
50: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
51: PetscCall(PetscFree(N->data));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
56: {
57: Mat_Transpose *Na = (Mat_Transpose *)N->data;
59: PetscFunctionBegin;
60: if (op == MAT_COPY_VALUES) {
61: PetscCall(MatTranspose(Na->A, MAT_INITIAL_MATRIX, m));
62: } else if (op == MAT_DO_NOT_COPY_VALUES) {
63: PetscCall(MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m));
64: PetscCall(MatTranspose(*m, MAT_INPLACE_MATRIX, m));
65: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode MatCreateVecs_Transpose(Mat A, Vec *r, Vec *l)
70: {
71: Mat_Transpose *Aa = (Mat_Transpose *)A->data;
73: PetscFunctionBegin;
74: PetscCall(MatCreateVecs(Aa->A, l, r));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode MatAXPY_Transpose(Mat Y, PetscScalar a, Mat X, MatStructure str)
79: {
80: Mat_Transpose *Ya = (Mat_Transpose *)Y->data;
81: Mat_Transpose *Xa = (Mat_Transpose *)X->data;
82: Mat M = Ya->A;
83: Mat N = Xa->A;
85: PetscFunctionBegin;
86: PetscCall(MatAXPY(M, a, N, str));
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
91: {
92: Mat_Transpose *X = (Mat_Transpose *)mat->data;
93: PetscFunctionBegin;
95: *has = PETSC_FALSE;
96: if (op == MATOP_MULT) {
97: PetscCall(MatHasOperation(X->A, MATOP_MULT_TRANSPOSE, has));
98: } else if (op == MATOP_MULT_TRANSPOSE) {
99: PetscCall(MatHasOperation(X->A, MATOP_MULT, has));
100: } else if (op == MATOP_MULT_ADD) {
101: PetscCall(MatHasOperation(X->A, MATOP_MULT_TRANSPOSE_ADD, has));
102: } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
103: PetscCall(MatHasOperation(X->A, MATOP_MULT_ADD, has));
104: } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
109: {
110: Mat A, B, C, Ain, Bin, Cin;
111: PetscBool Aistrans, Bistrans, Cistrans;
112: PetscInt Atrans, Btrans, Ctrans;
113: MatProductType ptype;
115: PetscFunctionBegin;
116: MatCheckProduct(D, 1);
117: A = D->product->A;
118: B = D->product->B;
119: C = D->product->C;
120: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
121: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
122: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
123: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
124: Atrans = 0;
125: Ain = A;
126: while (Aistrans) {
127: Atrans++;
128: PetscCall(MatTransposeGetMat(Ain, &Ain));
129: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
130: }
131: Btrans = 0;
132: Bin = B;
133: while (Bistrans) {
134: Btrans++;
135: PetscCall(MatTransposeGetMat(Bin, &Bin));
136: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
137: }
138: Ctrans = 0;
139: Cin = C;
140: while (Cistrans) {
141: Ctrans++;
142: PetscCall(MatTransposeGetMat(Cin, &Cin));
143: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
144: }
145: Atrans = Atrans % 2;
146: Btrans = Btrans % 2;
147: Ctrans = Ctrans % 2;
148: ptype = D->product->type; /* same product type by default */
149: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
150: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
151: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
153: if (Atrans || Btrans || Ctrans) {
154: ptype = MATPRODUCT_UNSPECIFIED;
155: switch (D->product->type) {
156: case MATPRODUCT_AB:
157: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
158: /* TODO custom implementation ? */
159: } else if (Atrans) { /* At * B */
160: ptype = MATPRODUCT_AtB;
161: } else { /* A * Bt */
162: ptype = MATPRODUCT_ABt;
163: }
164: break;
165: case MATPRODUCT_AtB:
166: if (Atrans && Btrans) { /* A * Bt */
167: ptype = MATPRODUCT_ABt;
168: } else if (Atrans) { /* A * B */
169: ptype = MATPRODUCT_AB;
170: } else { /* At * Bt we do not have support for this */
171: /* TODO custom implementation ? */
172: }
173: break;
174: case MATPRODUCT_ABt:
175: if (Atrans && Btrans) { /* At * B */
176: ptype = MATPRODUCT_AtB;
177: } else if (Atrans) { /* At * Bt we do not have support for this */
178: /* TODO custom implementation ? */
179: } else { /* A * B */
180: ptype = MATPRODUCT_AB;
181: }
182: break;
183: case MATPRODUCT_PtAP:
184: if (Atrans) { /* PtAtP */
185: /* TODO custom implementation ? */
186: } else { /* RARt */
187: ptype = MATPRODUCT_RARt;
188: }
189: break;
190: case MATPRODUCT_RARt:
191: if (Atrans) { /* RAtRt */
192: /* TODO custom implementation ? */
193: } else { /* PtAP */
194: ptype = MATPRODUCT_PtAP;
195: }
196: break;
197: case MATPRODUCT_ABC:
198: /* TODO custom implementation ? */
199: break;
200: default:
201: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
202: }
203: }
204: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
205: PetscCall(MatProductSetType(D, ptype));
206: PetscCall(MatProductSetFromOptions(D));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: static PetscErrorCode MatGetDiagonal_Transpose(Mat A, Vec v)
211: {
212: Mat_Transpose *Aa = (Mat_Transpose *)A->data;
214: PetscFunctionBegin;
215: PetscCall(MatGetDiagonal(Aa->A, v));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: static PetscErrorCode MatConvert_Transpose(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
220: {
221: Mat_Transpose *Aa = (Mat_Transpose *)A->data;
222: PetscBool flg;
224: PetscFunctionBegin;
225: PetscCall(MatHasOperation(Aa->A, MATOP_TRANSPOSE, &flg));
226: if (flg) {
227: Mat B;
229: PetscCall(MatTranspose(Aa->A, MAT_INITIAL_MATRIX, &B));
230: if (reuse != MAT_INPLACE_MATRIX) {
231: PetscCall(MatConvert(B, newtype, reuse, newmat));
232: PetscCall(MatDestroy(&B));
233: } else {
234: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
235: PetscCall(MatHeaderReplace(A, &B));
236: }
237: } else { /* use basic converter as fallback */
238: PetscCall(MatConvert_Basic(A, newtype, reuse, newmat));
239: }
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode MatTransposeGetMat_Transpose(Mat A, Mat *M)
244: {
245: Mat_Transpose *Aa = (Mat_Transpose *)A->data;
247: PetscFunctionBegin;
248: *M = Aa->A;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /*@
253: MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`
255: Logically Collective
257: Input Parameter:
258: . A - the `MATTRANSPOSEVIRTUAL` matrix
260: Output Parameter:
261: . M - the matrix object stored inside `A`
263: Level: intermediate
265: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
266: @*/
267: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
268: {
269: PetscFunctionBegin;
272: PetscAssertPointer(M, 2);
273: PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*MC
278: MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix
280: Level: advanced
282: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
283: `MATNORMALHERMITIAN`, `MATNORMAL`
284: M*/
286: /*@
287: MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'
289: Collective
291: Input Parameter:
292: . A - the (possibly rectangular) matrix
294: Output Parameter:
295: . N - the matrix that represents A'
297: Level: intermediate
299: Note:
300: The transpose A' is NOT actually formed! Rather the new matrix
301: object performs the matrix-vector product by using the `MatMultTranspose()` on
302: the original matrix
304: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
305: `MATNORMALHERMITIAN`
306: @*/
307: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
308: {
309: Mat_Transpose *Na;
310: VecType vtype;
312: PetscFunctionBegin;
313: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
314: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
315: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
316: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
318: PetscCall(PetscNew(&Na));
319: (*N)->data = (void *)Na;
320: PetscCall(PetscObjectReference((PetscObject)A));
321: Na->A = A;
323: (*N)->ops->destroy = MatDestroy_Transpose;
324: (*N)->ops->mult = MatMult_Transpose;
325: (*N)->ops->multadd = MatMultAdd_Transpose;
326: (*N)->ops->multtranspose = MatMultTranspose_Transpose;
327: (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose;
328: (*N)->ops->duplicate = MatDuplicate_Transpose;
329: (*N)->ops->getvecs = MatCreateVecs_Transpose;
330: (*N)->ops->axpy = MatAXPY_Transpose;
331: (*N)->ops->hasoperation = MatHasOperation_Transpose;
332: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
333: (*N)->ops->getdiagonal = MatGetDiagonal_Transpose;
334: (*N)->ops->convert = MatConvert_Transpose;
335: (*N)->assembled = PETSC_TRUE;
337: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
338: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
339: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
340: PetscCall(MatGetVecType(A, &vtype));
341: PetscCall(MatSetVecType(*N, vtype));
342: #if defined(PETSC_HAVE_DEVICE)
343: PetscCall(MatBindToCPU(*N, A->boundtocpu));
344: #endif
345: PetscCall(MatSetUp(*N));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }