Actual source code: htransm.c
2: #include <petsc/private/matimpl.h>
4: typedef struct {
5: Mat A;
6: } Mat_HT;
8: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_HermitianTranspose(Mat D)
9: {
10: Mat A, B, C, Ain, Bin, Cin;
11: PetscBool Aistrans, Bistrans, Cistrans;
12: PetscInt Atrans, Btrans, Ctrans;
13: MatProductType ptype;
15: PetscFunctionBegin;
16: MatCheckProduct(D, 1);
17: A = D->product->A;
18: B = D->product->B;
19: C = D->product->C;
20: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
21: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
22: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
23: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
24: Atrans = 0;
25: Ain = A;
26: while (Aistrans) {
27: Atrans++;
28: PetscCall(MatHermitianTransposeGetMat(Ain, &Ain));
29: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
30: }
31: Btrans = 0;
32: Bin = B;
33: while (Bistrans) {
34: Btrans++;
35: PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
36: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
37: }
38: Ctrans = 0;
39: Cin = C;
40: while (Cistrans) {
41: Ctrans++;
42: PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
43: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
44: }
45: Atrans = Atrans % 2;
46: Btrans = Btrans % 2;
47: Ctrans = Ctrans % 2;
48: ptype = D->product->type; /* same product type by default */
49: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
50: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
51: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
53: if (Atrans || Btrans || Ctrans) {
54: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for complex Hermitian transpose matrices");
55: ptype = MATPRODUCT_UNSPECIFIED;
56: switch (D->product->type) {
57: case MATPRODUCT_AB:
58: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
59: /* TODO custom implementation ? */
60: } else if (Atrans) { /* At * B */
61: ptype = MATPRODUCT_AtB;
62: } else { /* A * Bt */
63: ptype = MATPRODUCT_ABt;
64: }
65: break;
66: case MATPRODUCT_AtB:
67: if (Atrans && Btrans) { /* A * Bt */
68: ptype = MATPRODUCT_ABt;
69: } else if (Atrans) { /* A * B */
70: ptype = MATPRODUCT_AB;
71: } else { /* At * Bt we do not have support for this */
72: /* TODO custom implementation ? */
73: }
74: break;
75: case MATPRODUCT_ABt:
76: if (Atrans && Btrans) { /* At * B */
77: ptype = MATPRODUCT_AtB;
78: } else if (Atrans) { /* At * Bt we do not have support for this */
79: /* TODO custom implementation ? */
80: } else { /* A * B */
81: ptype = MATPRODUCT_AB;
82: }
83: break;
84: case MATPRODUCT_PtAP:
85: if (Atrans) { /* PtAtP */
86: /* TODO custom implementation ? */
87: } else { /* RARt */
88: ptype = MATPRODUCT_RARt;
89: }
90: break;
91: case MATPRODUCT_RARt:
92: if (Atrans) { /* RAtRt */
93: /* TODO custom implementation ? */
94: } else { /* PtAP */
95: ptype = MATPRODUCT_PtAP;
96: }
97: break;
98: case MATPRODUCT_ABC:
99: /* TODO custom implementation ? */
100: break;
101: default:
102: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
103: }
104: }
105: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
106: PetscCall(MatProductSetType(D, ptype));
107: PetscCall(MatProductSetFromOptions(D));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
110: PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
111: {
112: Mat_HT *Na = (Mat_HT *)N->data;
114: PetscFunctionBegin;
115: PetscCall(MatMultHermitianTranspose(Na->A, x, y));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: PetscErrorCode MatMultAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
120: {
121: Mat_HT *Na = (Mat_HT *)N->data;
123: PetscFunctionBegin;
124: PetscCall(MatMultHermitianTransposeAdd(Na->A, v1, v2, v3));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
129: {
130: Mat_HT *Na = (Mat_HT *)N->data;
132: PetscFunctionBegin;
133: PetscCall(MatMult(Na->A, x, y));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
138: {
139: Mat_HT *Na = (Mat_HT *)N->data;
141: PetscFunctionBegin;
142: PetscCall(MatMultAdd(Na->A, v1, v2, v3));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: PetscErrorCode MatDestroy_HT(Mat N)
147: {
148: Mat_HT *Na = (Mat_HT *)N->data;
150: PetscFunctionBegin;
151: PetscCall(MatDestroy(&Na->A));
152: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
153: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
154: #if !defined(PETSC_USE_COMPLEX)
155: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
156: #endif
157: PetscCall(PetscFree(N->data));
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
162: {
163: Mat_HT *Na = (Mat_HT *)N->data;
165: PetscFunctionBegin;
166: if (op == MAT_COPY_VALUES) {
167: PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, m));
168: } else if (op == MAT_DO_NOT_COPY_VALUES) {
169: PetscCall(MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m));
170: PetscCall(MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m));
171: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l)
176: {
177: Mat_HT *Na = (Mat_HT *)N->data;
179: PetscFunctionBegin;
180: PetscCall(MatCreateVecs(Na->A, l, r));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: PetscErrorCode MatAXPY_HT(Mat Y, PetscScalar a, Mat X, MatStructure str)
185: {
186: Mat_HT *Ya = (Mat_HT *)Y->data;
187: Mat_HT *Xa = (Mat_HT *)X->data;
188: Mat M = Ya->A;
189: Mat N = Xa->A;
191: PetscFunctionBegin;
192: PetscCall(MatAXPY(M, a, N, str));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
197: {
198: Mat_HT *Na = (Mat_HT *)N->data;
200: PetscFunctionBegin;
201: *M = Na->A;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /*@
206: MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`
208: Logically Collective
210: Input Parameter:
211: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix
213: Output Parameter:
214: . M - the matrix object stored inside A
216: Level: intermediate
218: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
219: @*/
220: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
221: {
222: PetscFunctionBegin;
226: PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);
232: PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v)
233: {
234: Mat_HT *Na = (Mat_HT *)A->data;
236: PetscFunctionBegin;
237: PetscCall(MatGetDiagonal(Na->A, v));
238: PetscCall(VecConjugate(v));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
243: {
244: Mat_HT *Na = (Mat_HT *)A->data;
245: PetscBool flg;
247: PetscFunctionBegin;
248: PetscCall(MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg));
249: if (flg) {
250: Mat B;
252: PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B));
253: if (reuse != MAT_INPLACE_MATRIX) {
254: PetscCall(MatConvert(B, newtype, reuse, newmat));
255: PetscCall(MatDestroy(&B));
256: } else {
257: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
258: PetscCall(MatHeaderReplace(A, &B));
259: }
260: } else { /* use basic converter as fallback */
261: PetscCall(MatConvert_Basic(A, newtype, reuse, newmat));
262: }
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*MC
267: MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix
269: Level: advanced
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: PetscInt m, n;
298: Mat_HT *Na;
299: VecType vtype;
301: PetscFunctionBegin;
302: PetscCall(MatGetLocalSize(A, &m, &n));
303: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
304: PetscCall(MatSetSizes(*N, n, m, PETSC_DECIDE, PETSC_DECIDE));
305: PetscCall(PetscLayoutSetUp((*N)->rmap));
306: PetscCall(PetscLayoutSetUp((*N)->cmap));
307: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
309: PetscCall(PetscNew(&Na));
310: (*N)->data = (void *)Na;
311: PetscCall(PetscObjectReference((PetscObject)A));
312: Na->A = A;
314: (*N)->ops->destroy = MatDestroy_HT;
315: (*N)->ops->mult = MatMult_HT;
316: (*N)->ops->multadd = MatMultAdd_HT;
317: (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT;
318: (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
319: #if !defined(PETSC_USE_COMPLEX)
320: (*N)->ops->multtranspose = MatMultHermitianTranspose_HT;
321: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT;
322: #endif
323: (*N)->ops->duplicate = MatDuplicate_HT;
324: (*N)->ops->getvecs = MatCreateVecs_HT;
325: (*N)->ops->axpy = MatAXPY_HT;
326: #if !defined(PETSC_USE_COMPLEX)
327: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
328: #endif
329: (*N)->ops->getdiagonal = MatGetDiagonal_HT;
330: (*N)->ops->convert = MatConvert_HT;
331: (*N)->assembled = PETSC_TRUE;
333: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
334: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HermitianTranspose));
335: #if !defined(PETSC_USE_COMPLEX)
336: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
337: #endif
338: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
339: PetscCall(MatGetVecType(A, &vtype));
340: PetscCall(MatSetVecType(*N, vtype));
341: #if defined(PETSC_HAVE_DEVICE)
342: PetscCall(MatBindToCPU(*N, A->boundtocpu));
343: #endif
344: PetscCall(MatSetUp(*N));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }