Actual source code: htransm.c
1: #include <../src/mat/impls/shell/shell.h>
3: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_HT(Mat D)
4: {
5: Mat A, B, C, Ain, Bin, Cin;
6: PetscBool Aistrans, Bistrans, Cistrans;
7: PetscInt Atrans, Btrans, Ctrans;
8: MatProductType ptype;
10: PetscFunctionBegin;
11: MatCheckProduct(D, 1);
12: A = D->product->A;
13: B = D->product->B;
14: C = D->product->C;
15: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
16: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
17: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
18: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
19: Atrans = 0;
20: Ain = A;
21: while (Aistrans) {
22: Atrans++;
23: PetscCall(MatHermitianTransposeGetMat(Ain, &Ain));
24: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
25: }
26: Btrans = 0;
27: Bin = B;
28: while (Bistrans) {
29: Btrans++;
30: PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
31: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
32: }
33: Ctrans = 0;
34: Cin = C;
35: while (Cistrans) {
36: Ctrans++;
37: PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
38: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
39: }
40: Atrans = Atrans % 2;
41: Btrans = Btrans % 2;
42: Ctrans = Ctrans % 2;
43: ptype = D->product->type; /* same product type by default */
44: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
45: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
46: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
48: if (Atrans || Btrans || Ctrans) {
49: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for complex Hermitian transpose matrices");
50: ptype = MATPRODUCT_UNSPECIFIED;
51: switch (D->product->type) {
52: case MATPRODUCT_AB:
53: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
54: /* TODO custom implementation ? */
55: } else if (Atrans) { /* At * B */
56: ptype = MATPRODUCT_AtB;
57: } else { /* A * Bt */
58: ptype = MATPRODUCT_ABt;
59: }
60: break;
61: case MATPRODUCT_AtB:
62: if (Atrans && Btrans) { /* A * Bt */
63: ptype = MATPRODUCT_ABt;
64: } else if (Atrans) { /* A * B */
65: ptype = MATPRODUCT_AB;
66: } else { /* At * Bt we do not have support for this */
67: /* TODO custom implementation ? */
68: }
69: break;
70: case MATPRODUCT_ABt:
71: if (Atrans && Btrans) { /* At * B */
72: ptype = MATPRODUCT_AtB;
73: } else if (Atrans) { /* At * Bt we do not have support for this */
74: /* TODO custom implementation ? */
75: } else { /* A * B */
76: ptype = MATPRODUCT_AB;
77: }
78: break;
79: case MATPRODUCT_PtAP:
80: if (Atrans) { /* PtAtP */
81: /* TODO custom implementation ? */
82: } else { /* RARt */
83: ptype = MATPRODUCT_RARt;
84: }
85: break;
86: case MATPRODUCT_RARt:
87: if (Atrans) { /* RAtRt */
88: /* TODO custom implementation ? */
89: } else { /* PtAP */
90: ptype = MATPRODUCT_PtAP;
91: }
92: break;
93: case MATPRODUCT_ABC:
94: /* TODO custom implementation ? */
95: break;
96: default:
97: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
98: }
99: }
100: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
101: PetscCall(MatProductSetType(D, ptype));
102: PetscCall(MatProductSetFromOptions(D));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
105: static PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
106: {
107: Mat A;
109: PetscFunctionBegin;
110: PetscCall(MatShellGetContext(N, &A));
111: PetscCall(MatMultHermitianTranspose(A, x, y));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
116: {
117: Mat A;
119: PetscFunctionBegin;
120: PetscCall(MatShellGetContext(N, &A));
121: PetscCall(MatMult(A, x, y));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
124: static PetscErrorCode MatDestroy_HT(Mat N)
125: {
126: Mat A;
128: PetscFunctionBegin;
129: PetscCall(MatShellGetContext(N, &A));
130: PetscCall(MatDestroy(&A));
131: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
132: #if !defined(PETSC_USE_COMPLEX)
133: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
134: #endif
135: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
136: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
141: {
142: Mat A, C;
144: PetscFunctionBegin;
145: PetscCall(MatShellGetContext(N, &A));
146: PetscCall(MatDuplicate(A, op, &C));
147: PetscCall(MatCreateHermitianTranspose(C, m));
148: PetscCall(MatDestroy(&C));
149: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has)
154: {
155: Mat A;
157: PetscFunctionBegin;
158: PetscCall(MatShellGetContext(mat, &A));
159: *has = PETSC_FALSE;
160: if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
161: PetscCall(MatHasOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, has));
162: if (!*has) PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
163: } else if (op == MATOP_MULT_HERMITIAN_TRANSPOSE || op == MATOP_MULT_HERMITIAN_TRANS_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
164: PetscCall(MatHasOperation(A, MATOP_MULT, has));
165: } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
170: {
171: PetscFunctionBegin;
172: PetscCall(MatShellGetContext(N, M));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*@
177: MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`
179: Logically Collective
181: Input Parameter:
182: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix
184: Output Parameter:
185: . M - the matrix object stored inside A
187: Level: intermediate
189: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
190: @*/
191: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
192: {
193: PetscFunctionBegin;
196: PetscAssertPointer(M, 2);
197: PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v)
202: {
203: Mat A;
205: PetscFunctionBegin;
206: PetscCall(MatShellGetContext(N, &A));
207: PetscCall(MatGetDiagonal(A, v));
208: PetscCall(VecConjugate(v));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode MatCopy_HT(Mat A, Mat B, MatStructure str)
213: {
214: Mat a, b;
216: PetscFunctionBegin;
217: PetscCall(MatShellGetContext(A, &a));
218: PetscCall(MatShellGetContext(B, &b));
219: PetscCall(MatCopy(a, b, str));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
224: {
225: Mat A;
226: PetscBool flg;
228: PetscFunctionBegin;
229: PetscCall(MatShellGetContext(N, &A));
230: PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg));
231: if (flg) {
232: Mat B;
234: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B));
235: if (reuse != MAT_INPLACE_MATRIX) {
236: PetscCall(MatConvert(B, newtype, reuse, newmat));
237: PetscCall(MatDestroy(&B));
238: } else {
239: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
240: PetscCall(MatHeaderReplace(N, &B));
241: }
242: } else { /* use basic converter as fallback */
243: PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*MC
249: MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix
251: Level: advanced
253: Developer Notes:
254: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
256: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
258: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`
259: M*/
261: /*@
262: MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'*
264: Collective
266: Input Parameter:
267: . A - the (possibly rectangular) matrix
269: Output Parameter:
270: . N - the matrix that represents A'*
272: Level: intermediate
274: Note:
275: The Hermitian transpose A' is NOT actually formed! Rather the new matrix
276: object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
277: the original matrix
279: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
280: `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
281: @*/
282: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
283: {
284: VecType vtype;
286: PetscFunctionBegin;
287: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
288: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
289: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
290: PetscCall(MatSetType(*N, MATSHELL));
291: PetscCall(MatShellSetContext(*N, A));
292: PetscCall(PetscObjectReference((PetscObject)A));
294: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
295: PetscCall(MatGetVecType(A, &vtype));
296: PetscCall(MatSetVecType(*N, vtype));
297: #if defined(PETSC_HAVE_DEVICE)
298: PetscCall(MatBindToCPU(*N, A->boundtocpu));
299: #endif
300: PetscCall(MatSetUp(*N));
302: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT));
303: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT));
304: PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
305: #if !defined(PETSC_USE_COMPLEX)
306: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
307: #endif
308: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT));
309: PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_HT));
310: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT));
311: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_HT));
312: PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT));
314: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
315: #if !defined(PETSC_USE_COMPLEX)
316: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
317: #endif
318: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT));
319: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
320: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
321: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
322: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }