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(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: }
106: static PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
107: {
108: Mat A;
110: PetscFunctionBegin;
111: PetscCall(MatShellGetContext(N, &A));
112: PetscCall(MatMultHermitianTranspose(A, x, y));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
117: {
118: Mat A;
120: PetscFunctionBegin;
121: PetscCall(MatShellGetContext(N, &A));
122: PetscCall(MatMult(A, x, y));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode MatSolve_HT_LU(Mat N, Vec b, Vec x)
127: {
128: Mat A;
129: Vec w;
131: PetscFunctionBegin;
132: PetscCall(MatShellGetContext(N, &A));
133: PetscCall(VecDuplicate(b, &w));
134: PetscCall(VecCopy(b, w));
135: PetscCall(VecConjugate(w));
136: PetscCall(MatSolveTranspose(A, w, x));
137: PetscCall(VecConjugate(x));
138: PetscCall(VecDestroy(&w));
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: static PetscErrorCode MatSolveAdd_HT_LU(Mat N, Vec b, Vec y, Vec x)
143: {
144: Mat A;
145: Vec v, w;
147: PetscFunctionBegin;
148: PetscCall(MatShellGetContext(N, &A));
149: PetscCall(VecDuplicate(b, &v));
150: PetscCall(VecDuplicate(b, &w));
151: PetscCall(VecCopy(y, v));
152: PetscCall(VecCopy(b, w));
153: PetscCall(VecConjugate(v));
154: PetscCall(VecConjugate(w));
155: PetscCall(MatSolveTransposeAdd(A, w, v, x));
156: PetscCall(VecConjugate(x));
157: PetscCall(VecDestroy(&v));
158: PetscCall(VecDestroy(&w));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode MatMatSolve_HT_LU(Mat N, Mat B, Mat X)
163: {
164: Mat A, W;
166: PetscFunctionBegin;
167: PetscCall(MatShellGetContext(N, &A));
168: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &W));
169: PetscCall(MatConjugate(W));
170: PetscCall(MatMatSolveTranspose(A, W, X));
171: PetscCall(MatConjugate(X));
172: PetscCall(MatDestroy(&W));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode MatLUFactor_HT(Mat N, IS row, IS col, const MatFactorInfo *minfo)
177: {
178: Mat A;
180: PetscFunctionBegin;
181: PetscCall(MatShellGetContext(N, &A));
182: PetscCall(MatLUFactor(A, col, row, minfo));
183: PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (void (*)(void))MatSolve_HT_LU));
184: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_LU));
185: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_LU));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode MatSolve_HT_Cholesky(Mat N, Vec b, Vec x)
190: {
191: Mat A;
193: PetscFunctionBegin;
194: PetscCall(MatShellGetContext(N, &A));
195: PetscCall(MatSolve(A, b, x));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode MatSolveAdd_HT_Cholesky(Mat N, Vec b, Vec y, Vec x)
200: {
201: Mat A;
202: Vec v, w;
204: PetscFunctionBegin;
205: PetscCall(MatShellGetContext(N, &A));
206: PetscCall(VecDuplicate(b, &v));
207: PetscCall(VecDuplicate(b, &w));
208: PetscCall(VecCopy(y, v));
209: PetscCall(VecCopy(b, w));
210: PetscCall(VecConjugate(v));
211: PetscCall(VecConjugate(w));
212: PetscCall(MatSolveTransposeAdd(A, w, v, x));
213: PetscCall(VecConjugate(x));
214: PetscCall(VecDestroy(&v));
215: PetscCall(VecDestroy(&w));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: static PetscErrorCode MatMatSolve_HT_Cholesky(Mat N, Mat B, Mat X)
220: {
221: Mat A, W;
223: PetscFunctionBegin;
224: PetscCall(MatShellGetContext(N, &A));
225: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &W));
226: PetscCall(MatConjugate(W));
227: PetscCall(MatMatSolveTranspose(A, W, X));
228: PetscCall(MatConjugate(X));
229: PetscCall(MatDestroy(&W));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode MatCholeskyFactor_HT(Mat N, IS perm, const MatFactorInfo *minfo)
234: {
235: Mat A;
237: PetscFunctionBegin;
238: PetscCall(MatShellGetContext(N, &A));
239: PetscCheck(!PetscDefined(USE_COMPLEX) || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cholesky supported only if original matrix is Hermitian");
240: PetscCall(MatCholeskyFactor(A, perm, minfo));
241: PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (void (*)(void))MatSolve_HT_Cholesky));
242: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_Cholesky));
243: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_Cholesky));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: static PetscErrorCode MatLUFactorNumeric_HT(Mat F, Mat N, const MatFactorInfo *info)
248: {
249: Mat A, FA;
251: PetscFunctionBegin;
252: PetscCall(MatShellGetContext(N, &A));
253: PetscCall(MatShellGetContext(F, &FA));
254: PetscCall(MatLUFactorNumeric(FA, A, info));
255: PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (void (*)(void))MatSolve_HT_LU));
256: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_LU));
257: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_LU));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode MatLUFactorSymbolic_HT(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info)
262: {
263: Mat A, FA;
265: PetscFunctionBegin;
266: PetscCall(MatShellGetContext(N, &A));
267: PetscCall(MatShellGetContext(F, &FA));
268: PetscCall(MatLUFactorSymbolic(FA, A, row, col, info));
269: PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (void (*)(void))MatLUFactorNumeric_HT));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode MatCholeskyFactorNumeric_HT(Mat F, Mat N, const MatFactorInfo *info)
274: {
275: Mat A, FA;
277: PetscFunctionBegin;
278: PetscCall(MatShellGetContext(N, &A));
279: PetscCall(MatShellGetContext(F, &FA));
280: PetscCall(MatCholeskyFactorNumeric(FA, A, info));
281: PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (void (*)(void))MatSolve_HT_Cholesky));
282: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_Cholesky));
283: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_Cholesky));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode MatCholeskyFactorSymbolic_HT(Mat F, Mat N, IS perm, const MatFactorInfo *info)
288: {
289: Mat A, FA;
291: PetscFunctionBegin;
292: PetscCall(MatShellGetContext(N, &A));
293: PetscCall(MatShellGetContext(F, &FA));
294: PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info));
295: PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (void (*)(void))MatCholeskyFactorNumeric_HT));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode MatGetFactor_HT(Mat N, MatSolverType type, MatFactorType ftype, Mat *F)
300: {
301: Mat A, FA;
303: PetscFunctionBegin;
304: PetscCall(MatShellGetContext(N, &A));
305: PetscCall(MatGetFactor(A, type, ftype, &FA));
306: PetscCall(MatCreateTranspose(FA, F));
307: if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (void (*)(void))MatLUFactorSymbolic_HT));
308: else if (ftype == MAT_FACTOR_CHOLESKY) {
309: PetscCheck(!PetscDefined(USE_COMPLEX) || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cholesky supported only if original matrix is Hermitian");
310: PetscCall(MatPropagateSymmetryOptions(A, FA));
311: PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (void (*)(void))MatCholeskyFactorSymbolic_HT));
312: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]);
313: (*F)->factortype = ftype;
314: PetscCall(MatDestroy(&FA));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode MatDestroy_HT(Mat N)
319: {
320: Mat A;
322: PetscFunctionBegin;
323: PetscCall(MatShellGetContext(N, &A));
324: PetscCall(MatDestroy(&A));
325: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
326: #if !defined(PETSC_USE_COMPLEX)
327: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
328: #endif
329: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
330: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
331: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: static PetscErrorCode MatGetInfo_HT(Mat N, MatInfoType flag, MatInfo *info)
336: {
337: Mat A;
339: PetscFunctionBegin;
340: PetscCall(MatShellGetContext(N, &A));
341: PetscCall(MatGetInfo(A, flag, info));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: static PetscErrorCode MatFactorGetSolverType_HT(Mat N, MatSolverType *type)
346: {
347: Mat A;
349: PetscFunctionBegin;
350: PetscCall(MatShellGetContext(N, &A));
351: PetscCall(MatFactorGetSolverType(A, type));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
356: {
357: Mat A, C;
359: PetscFunctionBegin;
360: PetscCall(MatShellGetContext(N, &A));
361: PetscCall(MatDuplicate(A, op, &C));
362: PetscCall(MatCreateHermitianTranspose(C, m));
363: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
364: PetscCall(MatDestroy(&C));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has)
369: {
370: Mat A;
372: PetscFunctionBegin;
373: PetscCall(MatShellGetContext(mat, &A));
374: *has = PETSC_FALSE;
375: if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
376: PetscCall(MatHasOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, has));
377: if (!*has) PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
378: } else if (op == MATOP_MULT_HERMITIAN_TRANSPOSE || op == MATOP_MULT_HERMITIAN_TRANS_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
379: PetscCall(MatHasOperation(A, MATOP_MULT, has));
380: } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
385: {
386: PetscFunctionBegin;
387: PetscCall(MatShellGetContext(N, M));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`
394: Logically Collective
396: Input Parameter:
397: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix
399: Output Parameter:
400: . M - the matrix object stored inside A
402: Level: intermediate
404: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
405: @*/
406: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
407: {
408: PetscFunctionBegin;
411: PetscAssertPointer(M, 2);
412: PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v)
417: {
418: Mat A;
420: PetscFunctionBegin;
421: PetscCall(MatShellGetContext(N, &A));
422: PetscCall(MatGetDiagonal(A, v));
423: PetscCall(VecConjugate(v));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: static PetscErrorCode MatCopy_HT(Mat A, Mat B, MatStructure str)
428: {
429: Mat a, b;
431: PetscFunctionBegin;
432: PetscCall(MatShellGetContext(A, &a));
433: PetscCall(MatShellGetContext(B, &b));
434: PetscCall(MatCopy(a, b, str));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
439: {
440: Mat A;
441: PetscScalar vscale = 1.0, vshift = 0.0;
442: PetscBool flg;
444: PetscFunctionBegin;
445: PetscCall(MatShellGetContext(N, &A));
446: PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg));
447: 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 */
448: PetscCheck(!((Mat_Shell *)N->data)->zrows && !((Mat_Shell *)N->data)->zcols, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat");
449: PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatAXPY() has been called on the input Mat");
450: PetscCheck(!((Mat_Shell *)N->data)->left && !((Mat_Shell *)N->data)->right, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalScale() has been called on the input Mat");
451: PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalSet() has been called on the input Mat");
452: vscale = ((Mat_Shell *)N->data)->vscale;
453: vshift = ((Mat_Shell *)N->data)->vshift;
454: }
455: if (flg) {
456: Mat B;
458: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B));
459: if (reuse != MAT_INPLACE_MATRIX) {
460: PetscCall(MatConvert(B, newtype, reuse, newmat));
461: PetscCall(MatDestroy(&B));
462: } else {
463: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
464: PetscCall(MatHeaderReplace(N, &B));
465: }
466: } else { /* use basic converter as fallback */
467: flg = (PetscBool)(N->ops->getrow != NULL);
468: PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
469: }
470: if (flg) {
471: PetscCall(MatScale(*newmat, vscale));
472: PetscCall(MatShift(*newmat, vshift));
473: }
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: /*MC
478: MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix
480: Level: advanced
482: Developer Notes:
483: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
485: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
487: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`
488: M*/
490: /*@
491: MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'*
493: Collective
495: Input Parameter:
496: . A - the (possibly rectangular) matrix
498: Output Parameter:
499: . N - the matrix that represents A'*
501: Level: intermediate
503: Note:
504: The Hermitian transpose A' is NOT actually formed! Rather the new matrix
505: object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
506: the original matrix
508: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
509: `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
510: @*/
511: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
512: {
513: VecType vtype;
515: PetscFunctionBegin;
516: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
517: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
518: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
519: PetscCall(MatSetType(*N, MATSHELL));
520: PetscCall(MatShellSetContext(*N, A));
521: PetscCall(PetscObjectReference((PetscObject)A));
523: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
524: PetscCall(MatGetVecType(A, &vtype));
525: PetscCall(MatSetVecType(*N, vtype));
526: #if defined(PETSC_HAVE_DEVICE)
527: PetscCall(MatBindToCPU(*N, A->boundtocpu));
528: #endif
529: PetscCall(MatSetUp(*N));
531: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT));
532: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT));
533: PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
534: #if !defined(PETSC_USE_COMPLEX)
535: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
536: #endif
537: PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (void (*)(void))MatLUFactor_HT));
538: PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (void (*)(void))MatCholeskyFactor_HT));
539: PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (void (*)(void))MatGetFactor_HT));
540: PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (void (*)(void))MatGetInfo_HT));
541: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT));
542: PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_HT));
543: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT));
544: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_HT));
545: PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT));
547: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
548: #if !defined(PETSC_USE_COMPLEX)
549: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
550: #endif
551: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT));
552: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_HT));
553: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
554: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
555: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
556: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }