Actual source code: htransm.c
1: #include <petsc/private/matimpl.h>
3: typedef struct {
4: Mat A;
5: } Mat_HT;
7: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_HermitianTranspose(Mat D)
8: {
9: Mat A, B, C, Ain, Bin, Cin;
10: PetscBool Aistrans, Bistrans, Cistrans;
11: PetscInt Atrans, Btrans, Ctrans;
12: MatProductType ptype;
14: PetscFunctionBegin;
15: MatCheckProduct(D, 1);
16: A = D->product->A;
17: B = D->product->B;
18: C = D->product->C;
19: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
20: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
21: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
22: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
23: Atrans = 0;
24: Ain = A;
25: while (Aistrans) {
26: Atrans++;
27: PetscCall(MatHermitianTransposeGetMat(Ain, &Ain));
28: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
29: }
30: Btrans = 0;
31: Bin = B;
32: while (Bistrans) {
33: Btrans++;
34: PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
35: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
36: }
37: Ctrans = 0;
38: Cin = C;
39: while (Cistrans) {
40: Ctrans++;
41: PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
42: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
43: }
44: Atrans = Atrans % 2;
45: Btrans = Btrans % 2;
46: Ctrans = Ctrans % 2;
47: ptype = D->product->type; /* same product type by default */
48: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
49: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
50: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
52: if (Atrans || Btrans || Ctrans) {
53: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for complex Hermitian transpose matrices");
54: ptype = MATPRODUCT_UNSPECIFIED;
55: switch (D->product->type) {
56: case MATPRODUCT_AB:
57: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
58: /* TODO custom implementation ? */
59: } else if (Atrans) { /* At * B */
60: ptype = MATPRODUCT_AtB;
61: } else { /* A * Bt */
62: ptype = MATPRODUCT_ABt;
63: }
64: break;
65: case MATPRODUCT_AtB:
66: if (Atrans && Btrans) { /* A * Bt */
67: ptype = MATPRODUCT_ABt;
68: } else if (Atrans) { /* A * B */
69: ptype = MATPRODUCT_AB;
70: } else { /* At * Bt we do not have support for this */
71: /* TODO custom implementation ? */
72: }
73: break;
74: case MATPRODUCT_ABt:
75: if (Atrans && Btrans) { /* At * B */
76: ptype = MATPRODUCT_AtB;
77: } else if (Atrans) { /* At * Bt we do not have support for this */
78: /* TODO custom implementation ? */
79: } else { /* A * B */
80: ptype = MATPRODUCT_AB;
81: }
82: break;
83: case MATPRODUCT_PtAP:
84: if (Atrans) { /* PtAtP */
85: /* TODO custom implementation ? */
86: } else { /* RARt */
87: ptype = MATPRODUCT_RARt;
88: }
89: break;
90: case MATPRODUCT_RARt:
91: if (Atrans) { /* RAtRt */
92: /* TODO custom implementation ? */
93: } else { /* PtAP */
94: ptype = MATPRODUCT_PtAP;
95: }
96: break;
97: case MATPRODUCT_ABC:
98: /* TODO custom implementation ? */
99: break;
100: default:
101: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
102: }
103: }
104: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
105: PetscCall(MatProductSetType(D, ptype));
106: PetscCall(MatProductSetFromOptions(D));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
109: static PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
110: {
111: Mat_HT *Na = (Mat_HT *)N->data;
113: PetscFunctionBegin;
114: PetscCall(MatMultHermitianTranspose(Na->A, x, y));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode MatMultAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
119: {
120: Mat_HT *Na = (Mat_HT *)N->data;
122: PetscFunctionBegin;
123: PetscCall(MatMultHermitianTransposeAdd(Na->A, v1, v2, v3));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
128: {
129: Mat_HT *Na = (Mat_HT *)N->data;
131: PetscFunctionBegin;
132: PetscCall(MatMult(Na->A, x, y));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: static PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
137: {
138: Mat_HT *Na = (Mat_HT *)N->data;
140: PetscFunctionBegin;
141: PetscCall(MatMultAdd(Na->A, v1, v2, v3));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode MatDestroy_HT(Mat N)
146: {
147: Mat_HT *Na = (Mat_HT *)N->data;
149: PetscFunctionBegin;
150: PetscCall(MatDestroy(&Na->A));
151: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
152: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
153: #if !defined(PETSC_USE_COMPLEX)
154: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
155: #endif
156: PetscCall(PetscFree(N->data));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
161: {
162: Mat_HT *Na = (Mat_HT *)N->data;
164: PetscFunctionBegin;
165: if (op == MAT_COPY_VALUES) {
166: PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, m));
167: } else if (op == MAT_DO_NOT_COPY_VALUES) {
168: PetscCall(MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m));
169: PetscCall(MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m));
170: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l)
175: {
176: Mat_HT *Na = (Mat_HT *)N->data;
178: PetscFunctionBegin;
179: PetscCall(MatCreateVecs(Na->A, l, r));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: static PetscErrorCode MatAXPY_HT(Mat Y, PetscScalar a, Mat X, MatStructure str)
184: {
185: Mat_HT *Ya = (Mat_HT *)Y->data;
186: Mat_HT *Xa = (Mat_HT *)X->data;
187: Mat M = Ya->A;
188: Mat N = Xa->A;
190: PetscFunctionBegin;
191: PetscCall(MatAXPY(M, a, N, str));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
196: {
197: Mat_HT *Na = (Mat_HT *)N->data;
199: PetscFunctionBegin;
200: *M = Na->A;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@
205: MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`
207: Logically Collective
209: Input Parameter:
210: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix
212: Output Parameter:
213: . M - the matrix object stored inside A
215: Level: intermediate
217: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
218: @*/
219: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
220: {
221: PetscFunctionBegin;
224: PetscAssertPointer(M, 2);
225: PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);
231: static PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v)
232: {
233: Mat_HT *Na = (Mat_HT *)A->data;
235: PetscFunctionBegin;
236: PetscCall(MatGetDiagonal(Na->A, v));
237: PetscCall(VecConjugate(v));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
242: {
243: Mat_HT *Na = (Mat_HT *)A->data;
244: PetscBool flg;
246: PetscFunctionBegin;
247: PetscCall(MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg));
248: if (flg) {
249: Mat B;
251: PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B));
252: if (reuse != MAT_INPLACE_MATRIX) {
253: PetscCall(MatConvert(B, newtype, reuse, newmat));
254: PetscCall(MatDestroy(&B));
255: } else {
256: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
257: PetscCall(MatHeaderReplace(A, &B));
258: }
259: } else { /* use basic converter as fallback */
260: PetscCall(MatConvert_Basic(A, newtype, reuse, newmat));
261: }
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*MC
266: MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix
268: Level: advanced
270: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`
271: M*/
273: /*@
274: MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'*
276: Collective
278: Input Parameter:
279: . A - the (possibly rectangular) matrix
281: Output Parameter:
282: . N - the matrix that represents A'*
284: Level: intermediate
286: Note:
287: The Hermitian transpose A' is NOT actually formed! Rather the new matrix
288: object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
289: the original matrix
291: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
292: `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
293: @*/
294: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
295: {
296: Mat_HT *Na;
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(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
305: PetscCall(PetscNew(&Na));
306: (*N)->data = (void *)Na;
307: PetscCall(PetscObjectReference((PetscObject)A));
308: Na->A = A;
310: (*N)->ops->destroy = MatDestroy_HT;
311: (*N)->ops->mult = MatMult_HT;
312: (*N)->ops->multadd = MatMultAdd_HT;
313: (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT;
314: (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
315: #if !defined(PETSC_USE_COMPLEX)
316: (*N)->ops->multtranspose = MatMultHermitianTranspose_HT;
317: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT;
318: #endif
319: (*N)->ops->duplicate = MatDuplicate_HT;
320: (*N)->ops->getvecs = MatCreateVecs_HT;
321: (*N)->ops->axpy = MatAXPY_HT;
322: #if !defined(PETSC_USE_COMPLEX)
323: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
324: #endif
325: (*N)->ops->getdiagonal = MatGetDiagonal_HT;
326: (*N)->ops->convert = MatConvert_HT;
327: (*N)->assembled = PETSC_TRUE;
329: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
330: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HermitianTranspose));
331: #if !defined(PETSC_USE_COMPLEX)
332: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
333: #endif
334: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
335: PetscCall(MatGetVecType(A, &vtype));
336: PetscCall(MatSetVecType(*N, vtype));
337: #if defined(PETSC_HAVE_DEVICE)
338: PetscCall(MatBindToCPU(*N, A->boundtocpu));
339: #endif
340: PetscCall(MatSetUp(*N));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }