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: MatCheckProduct(D, 1);
16: A = D->product->A;
17: B = D->product->B;
18: C = D->product->C;
19: PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans);
20: PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans);
21: PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans);
23: Atrans = 0;
24: Ain = A;
25: while (Aistrans) {
26: Atrans++;
27: MatHermitianTransposeGetMat(Ain, &Ain);
28: PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans);
29: }
30: Btrans = 0;
31: Bin = B;
32: while (Bistrans) {
33: Btrans++;
34: MatHermitianTransposeGetMat(Bin, &Bin);
35: PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans);
36: }
37: Ctrans = 0;
38: Cin = C;
39: while (Cistrans) {
40: Ctrans++;
41: MatHermitianTransposeGetMat(Cin, &Cin);
42: 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) {
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: MatProductReplaceMats(Ain, Bin, Cin, D);
105: MatProductSetType(D, ptype);
106: MatProductSetFromOptions(D);
107: return 0;
108: }
109: PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
110: {
111: Mat_HT *Na = (Mat_HT *)N->data;
113: MatMultHermitianTranspose(Na->A, x, y);
114: return 0;
115: }
117: PetscErrorCode MatMultAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
118: {
119: Mat_HT *Na = (Mat_HT *)N->data;
121: MatMultHermitianTransposeAdd(Na->A, v1, v2, v3);
122: return 0;
123: }
125: PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
126: {
127: Mat_HT *Na = (Mat_HT *)N->data;
129: MatMult(Na->A, x, y);
130: return 0;
131: }
133: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
134: {
135: Mat_HT *Na = (Mat_HT *)N->data;
137: MatMultAdd(Na->A, v1, v2, v3);
138: return 0;
139: }
141: PetscErrorCode MatDestroy_HT(Mat N)
142: {
143: Mat_HT *Na = (Mat_HT *)N->data;
145: MatDestroy(&Na->A);
146: PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL);
147: PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL);
148: #if !defined(PETSC_USE_COMPLEX)
149: PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL);
150: #endif
151: PetscFree(N->data);
152: return 0;
153: }
155: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
156: {
157: Mat_HT *Na = (Mat_HT *)N->data;
159: if (op == MAT_COPY_VALUES) {
160: MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, m);
161: } else if (op == MAT_DO_NOT_COPY_VALUES) {
162: MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m);
163: MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m);
164: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
165: return 0;
166: }
168: PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l)
169: {
170: Mat_HT *Na = (Mat_HT *)N->data;
172: MatCreateVecs(Na->A, l, r);
173: return 0;
174: }
176: PetscErrorCode MatAXPY_HT(Mat Y, PetscScalar a, Mat X, MatStructure str)
177: {
178: Mat_HT *Ya = (Mat_HT *)Y->data;
179: Mat_HT *Xa = (Mat_HT *)X->data;
180: Mat M = Ya->A;
181: Mat N = Xa->A;
183: MatAXPY(M, a, N, str);
184: return 0;
185: }
187: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
188: {
189: Mat_HT *Na = (Mat_HT *)N->data;
191: *M = Na->A;
192: return 0;
193: }
195: /*@
196: MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`
198: Logically collective
200: Input Parameter:
201: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix
203: Output Parameter:
204: . M - the matrix object stored inside A
206: Level: intermediate
208: .seealso: `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
209: @*/
210: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
211: {
215: PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
216: return 0;
217: }
219: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);
221: PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v)
222: {
223: Mat_HT *Na = (Mat_HT *)A->data;
225: MatGetDiagonal(Na->A, v);
226: VecConjugate(v);
227: return 0;
228: }
230: PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
231: {
232: Mat_HT *Na = (Mat_HT *)A->data;
233: PetscBool flg;
235: MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg);
236: if (flg) {
237: Mat B;
239: MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B);
240: if (reuse != MAT_INPLACE_MATRIX) {
241: MatConvert(B, newtype, reuse, newmat);
242: MatDestroy(&B);
243: } else {
244: MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B);
245: MatHeaderReplace(A, &B);
246: }
247: } else { /* use basic converter as fallback */
248: MatConvert_Basic(A, newtype, reuse, newmat);
249: }
250: return 0;
251: }
253: /*MC
254: MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix
256: Level: advanced
258: .seealso: `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: `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
280: `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`
281: @*/
282: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
283: {
284: PetscInt m, n;
285: Mat_HT *Na;
286: VecType vtype;
288: MatGetLocalSize(A, &m, &n);
289: MatCreate(PetscObjectComm((PetscObject)A), N);
290: MatSetSizes(*N, n, m, PETSC_DECIDE, PETSC_DECIDE);
291: PetscLayoutSetUp((*N)->rmap);
292: PetscLayoutSetUp((*N)->cmap);
293: PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL);
295: PetscNew(&Na);
296: (*N)->data = (void *)Na;
297: PetscObjectReference((PetscObject)A);
298: Na->A = A;
300: (*N)->ops->destroy = MatDestroy_HT;
301: (*N)->ops->mult = MatMult_HT;
302: (*N)->ops->multadd = MatMultAdd_HT;
303: (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT;
304: (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
305: #if !defined(PETSC_USE_COMPLEX)
306: (*N)->ops->multtranspose = MatMultHermitianTranspose_HT;
307: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT;
308: #endif
309: (*N)->ops->duplicate = MatDuplicate_HT;
310: (*N)->ops->getvecs = MatCreateVecs_HT;
311: (*N)->ops->axpy = MatAXPY_HT;
312: #if !defined(PETSC_USE_COMPLEX)
313: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
314: #endif
315: (*N)->ops->getdiagonal = MatGetDiagonal_HT;
316: (*N)->ops->convert = MatConvert_HT;
317: (*N)->assembled = PETSC_TRUE;
319: PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT);
320: PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HermitianTranspose);
321: #if !defined(PETSC_USE_COMPLEX)
322: PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT);
323: #endif
324: MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
325: MatGetVecType(A, &vtype);
326: MatSetVecType(*N, vtype);
327: #if defined(PETSC_HAVE_DEVICE)
328: MatBindToCPU(*N, A->boundtocpu);
329: #endif
330: MatSetUp(*N);
331: return 0;
332: }