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