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