Actual source code: htransm.c
1: #include <../src/mat/impls/shell/shell.h>
3: typedef struct {
4: PetscErrorCode (*numeric)(Mat);
5: PetscErrorCode (*destroy)(void *);
6: Mat B;
7: PetscScalar scale;
8: PetscBool conjugate;
9: PetscContainer container;
10: void *stash;
11: } MatProductData;
13: static PetscErrorCode DestroyMatProductData(void *ptr)
14: {
15: MatProductData *data = (MatProductData *)ptr;
17: PetscFunctionBegin;
18: if (data->stash) PetscCall((*data->destroy)(data->stash));
19: if (data->conjugate) PetscCall(MatDestroy(&data->B));
20: PetscCall(PetscContainerDestroy(&data->container));
21: PetscCall(PetscFree(data));
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: static PetscErrorCode MatProductNumeric_HT(Mat D)
26: {
27: Mat_Product *product;
28: Mat B;
29: MatProductData *data;
30: PetscContainer container;
32: PetscFunctionBegin;
33: MatCheckProduct(D, 1);
34: PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
35: product = D->product;
36: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductData", (PetscObject *)&container));
37: PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductData missing");
38: PetscCall(PetscContainerGetPointer(container, (void **)&data));
39: B = product->B;
40: data = (MatProductData *)product->data;
41: if (data->conjugate) {
42: PetscCall(MatCopy(product->B, data->B, SAME_NONZERO_PATTERN));
43: PetscCall(MatConjugate(data->B));
44: product->B = data->B;
45: }
46: product->data = data->stash;
47: PetscCall((*data->numeric)(D));
48: if (data->conjugate) {
49: PetscCall(MatConjugate(D));
50: product->B = B;
51: }
52: PetscCall(MatScale(D, data->scale));
53: product->data = data;
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode MatProductSymbolic_HT(Mat D)
58: {
59: Mat_Product *product;
60: Mat B;
61: MatProductData *data;
62: PetscContainer container;
64: PetscFunctionBegin;
65: MatCheckProduct(D, 1);
66: product = D->product;
67: B = product->B;
68: if (D->ops->productsymbolic == MatProductSymbolic_HT) {
69: PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
70: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductData", (PetscObject *)&container));
71: PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductData missing");
72: PetscCall(PetscContainerGetPointer(container, (void **)&data));
73: PetscCall(MatProductSetFromOptions(D));
74: if (data->conjugate) {
75: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &data->B));
76: product->B = data->B;
77: }
78: PetscCall(MatProductSymbolic(D));
79: data->numeric = D->ops->productnumeric;
80: data->destroy = product->destroy;
81: data->stash = product->data;
82: D->ops->productnumeric = MatProductNumeric_HT;
83: product->destroy = DestroyMatProductData;
84: if (data->conjugate) product->B = B;
85: product->data = data;
86: }
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: static PetscErrorCode MatProductSetFromOptions_HT(Mat D)
91: {
92: Mat A, B, C, Ain, Bin, Cin;
93: PetscScalar scale = 1.0, vscale;
94: PetscBool Aistrans, Bistrans, Cistrans, conjugate = PETSC_FALSE;
95: PetscInt Atrans, Btrans, Ctrans;
96: PetscContainer container = NULL;
97: MatProductData *data;
98: MatProductType ptype;
100: PetscFunctionBegin;
101: MatCheckProduct(D, 1);
102: A = D->product->A;
103: B = D->product->B;
104: C = D->product->C;
105: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
106: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
107: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
108: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
109: Atrans = 0;
110: Ain = A;
111: while (Aistrans) {
112: Atrans++;
113: PetscCall(MatShellGetScalingShifts(Ain, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &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));
114: conjugate = (PetscBool)!conjugate;
115: scale *= vscale;
116: PetscCall(MatHermitianTransposeGetMat(Ain, &Ain));
117: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
118: }
119: Btrans = 0;
120: Bin = B;
121: while (Bistrans) {
122: Btrans++;
123: PetscCall(MatShellGetScalingShifts(Bin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &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));
124: scale *= vscale;
125: PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
126: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
127: }
128: Ctrans = 0;
129: Cin = C;
130: while (Cistrans) {
131: Ctrans++;
132: PetscCall(MatShellGetScalingShifts(Cin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &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));
133: scale *= vscale;
134: PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
135: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
136: }
137: Atrans = Atrans % 2;
138: Btrans = Btrans % 2;
139: Ctrans = Ctrans % 2;
140: ptype = D->product->type; /* same product type by default */
141: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
142: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
143: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
145: if (Atrans || Btrans || Ctrans) {
146: PetscCheck(!PetscDefined(USE_COMPLEX) || (!Btrans && !Ctrans), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for complex Hermitian transpose matrices");
147: if ((PetscDefined(USE_COMPLEX) && Atrans) || scale != 1.0) {
148: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductData", (PetscObject *)&container));
149: if (!container) {
150: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)D), &container));
151: PetscCall(PetscNew(&data));
152: data->scale = scale;
153: data->conjugate = (PetscBool)Atrans;
154: data->container = container;
155: PetscCall(PetscContainerSetPointer(container, data));
156: PetscCall(PetscObjectCompose((PetscObject)D, "MatProductData", (PetscObject)container));
157: }
158: }
159: ptype = MATPRODUCT_UNSPECIFIED;
160: switch (D->product->type) {
161: case MATPRODUCT_AB:
162: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
163: /* TODO custom implementation ? */
164: } else if (Atrans) { /* At * B */
165: ptype = MATPRODUCT_AtB;
166: } else { /* A * Bt */
167: ptype = MATPRODUCT_ABt;
168: }
169: break;
170: case MATPRODUCT_AtB:
171: if (Atrans && Btrans) { /* A * Bt */
172: ptype = MATPRODUCT_ABt;
173: } else if (Atrans) { /* A * B */
174: ptype = MATPRODUCT_AB;
175: } else { /* At * Bt we do not have support for this */
176: /* TODO custom implementation ? */
177: }
178: break;
179: case MATPRODUCT_ABt:
180: if (Atrans && Btrans) { /* At * B */
181: ptype = MATPRODUCT_AtB;
182: } else if (Atrans) { /* At * Bt we do not have support for this */
183: /* TODO custom implementation ? */
184: } else { /* A * B */
185: ptype = MATPRODUCT_AB;
186: }
187: break;
188: case MATPRODUCT_PtAP:
189: if (Atrans) { /* PtAtP */
190: /* TODO custom implementation ? */
191: } else { /* RARt */
192: ptype = MATPRODUCT_RARt;
193: }
194: break;
195: case MATPRODUCT_RARt:
196: if (Atrans) { /* RAtRt */
197: /* TODO custom implementation ? */
198: } else { /* PtAP */
199: ptype = MATPRODUCT_PtAP;
200: }
201: break;
202: case MATPRODUCT_ABC:
203: /* TODO custom implementation ? */
204: break;
205: default:
206: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
207: }
208: }
209: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
210: PetscCall(MatProductSetType(D, ptype));
211: if (container) D->ops->productsymbolic = MatProductSymbolic_HT;
212: else PetscCall(MatProductSetFromOptions(D));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: static PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
217: {
218: Mat A;
220: PetscFunctionBegin;
221: PetscCall(MatShellGetContext(N, &A));
222: PetscCall(MatMultHermitianTranspose(A, x, y));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
227: {
228: Mat A;
230: PetscFunctionBegin;
231: PetscCall(MatShellGetContext(N, &A));
232: PetscCall(MatMult(A, x, y));
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: static PetscErrorCode MatSolve_HT_LU(Mat N, Vec b, Vec x)
237: {
238: Mat A;
239: Vec w;
241: PetscFunctionBegin;
242: PetscCall(MatShellGetContext(N, &A));
243: PetscCall(VecDuplicate(b, &w));
244: PetscCall(VecCopy(b, w));
245: PetscCall(VecConjugate(w));
246: PetscCall(MatSolveTranspose(A, w, x));
247: PetscCall(VecConjugate(x));
248: PetscCall(VecDestroy(&w));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: static PetscErrorCode MatSolveAdd_HT_LU(Mat N, Vec b, Vec y, Vec x)
253: {
254: Mat A;
255: Vec v, w;
257: PetscFunctionBegin;
258: PetscCall(MatShellGetContext(N, &A));
259: PetscCall(VecDuplicate(b, &v));
260: PetscCall(VecDuplicate(b, &w));
261: PetscCall(VecCopy(y, v));
262: PetscCall(VecCopy(b, w));
263: PetscCall(VecConjugate(v));
264: PetscCall(VecConjugate(w));
265: PetscCall(MatSolveTransposeAdd(A, w, v, x));
266: PetscCall(VecConjugate(x));
267: PetscCall(VecDestroy(&v));
268: PetscCall(VecDestroy(&w));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode MatMatSolve_HT_LU(Mat N, Mat B, Mat X)
273: {
274: Mat A, W;
276: PetscFunctionBegin;
277: PetscCall(MatShellGetContext(N, &A));
278: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &W));
279: PetscCall(MatConjugate(W));
280: PetscCall(MatMatSolveTranspose(A, W, X));
281: PetscCall(MatConjugate(X));
282: PetscCall(MatDestroy(&W));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: static PetscErrorCode MatLUFactor_HT(Mat N, IS row, IS col, const MatFactorInfo *minfo)
287: {
288: Mat A;
290: PetscFunctionBegin;
291: PetscCall(MatShellGetContext(N, &A));
292: PetscCall(MatLUFactor(A, col, row, minfo));
293: PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (void (*)(void))MatSolve_HT_LU));
294: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_LU));
295: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_LU));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode MatSolve_HT_Cholesky(Mat N, Vec b, Vec x)
300: {
301: Mat A;
303: PetscFunctionBegin;
304: PetscCall(MatShellGetContext(N, &A));
305: PetscCall(MatSolve(A, b, x));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode MatSolveAdd_HT_Cholesky(Mat N, Vec b, Vec y, Vec x)
310: {
311: Mat A;
312: Vec v, w;
314: PetscFunctionBegin;
315: PetscCall(MatShellGetContext(N, &A));
316: PetscCall(VecDuplicate(b, &v));
317: PetscCall(VecDuplicate(b, &w));
318: PetscCall(VecCopy(y, v));
319: PetscCall(VecCopy(b, w));
320: PetscCall(VecConjugate(v));
321: PetscCall(VecConjugate(w));
322: PetscCall(MatSolveTransposeAdd(A, w, v, x));
323: PetscCall(VecConjugate(x));
324: PetscCall(VecDestroy(&v));
325: PetscCall(VecDestroy(&w));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: static PetscErrorCode MatMatSolve_HT_Cholesky(Mat N, Mat B, Mat X)
330: {
331: Mat A, W;
333: PetscFunctionBegin;
334: PetscCall(MatShellGetContext(N, &A));
335: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &W));
336: PetscCall(MatConjugate(W));
337: PetscCall(MatMatSolveTranspose(A, W, X));
338: PetscCall(MatConjugate(X));
339: PetscCall(MatDestroy(&W));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: static PetscErrorCode MatCholeskyFactor_HT(Mat N, IS perm, const MatFactorInfo *minfo)
344: {
345: Mat A;
347: PetscFunctionBegin;
348: PetscCall(MatShellGetContext(N, &A));
349: PetscCheck(!PetscDefined(USE_COMPLEX) || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cholesky supported only if original matrix is Hermitian");
350: PetscCall(MatCholeskyFactor(A, perm, minfo));
351: PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (void (*)(void))MatSolve_HT_Cholesky));
352: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_Cholesky));
353: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_Cholesky));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode MatLUFactorNumeric_HT(Mat F, Mat N, const MatFactorInfo *info)
358: {
359: Mat A, FA;
361: PetscFunctionBegin;
362: PetscCall(MatShellGetContext(N, &A));
363: PetscCall(MatShellGetContext(F, &FA));
364: PetscCall(MatLUFactorNumeric(FA, A, info));
365: PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (void (*)(void))MatSolve_HT_LU));
366: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_LU));
367: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_LU));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode MatLUFactorSymbolic_HT(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info)
372: {
373: Mat A, FA;
375: PetscFunctionBegin;
376: PetscCall(MatShellGetContext(N, &A));
377: PetscCall(MatShellGetContext(F, &FA));
378: PetscCall(MatLUFactorSymbolic(FA, A, row, col, info));
379: PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (void (*)(void))MatLUFactorNumeric_HT));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode MatCholeskyFactorNumeric_HT(Mat F, Mat N, const MatFactorInfo *info)
384: {
385: Mat A, FA;
387: PetscFunctionBegin;
388: PetscCall(MatShellGetContext(N, &A));
389: PetscCall(MatShellGetContext(F, &FA));
390: PetscCall(MatCholeskyFactorNumeric(FA, A, info));
391: PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (void (*)(void))MatSolve_HT_Cholesky));
392: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_HT_Cholesky));
393: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_HT_Cholesky));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static PetscErrorCode MatCholeskyFactorSymbolic_HT(Mat F, Mat N, IS perm, const MatFactorInfo *info)
398: {
399: Mat A, FA;
401: PetscFunctionBegin;
402: PetscCall(MatShellGetContext(N, &A));
403: PetscCall(MatShellGetContext(F, &FA));
404: PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info));
405: PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (void (*)(void))MatCholeskyFactorNumeric_HT));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode MatGetFactor_HT(Mat N, MatSolverType type, MatFactorType ftype, Mat *F)
410: {
411: Mat A, FA;
413: PetscFunctionBegin;
414: PetscCall(MatShellGetContext(N, &A));
415: PetscCall(MatGetFactor(A, type, ftype, &FA));
416: PetscCall(MatCreateTranspose(FA, F));
417: if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (void (*)(void))MatLUFactorSymbolic_HT));
418: else if (ftype == MAT_FACTOR_CHOLESKY) {
419: PetscCheck(!PetscDefined(USE_COMPLEX) || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cholesky supported only if original matrix is Hermitian");
420: PetscCall(MatPropagateSymmetryOptions(A, FA));
421: PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (void (*)(void))MatCholeskyFactorSymbolic_HT));
422: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]);
423: (*F)->factortype = ftype;
424: PetscCall(MatDestroy(&FA));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode MatDestroy_HT(Mat N)
429: {
430: Mat A;
432: PetscFunctionBegin;
433: PetscCall(MatShellGetContext(N, &A));
434: PetscCall(MatDestroy(&A));
435: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
436: #if !defined(PETSC_USE_COMPLEX)
437: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
438: #endif
439: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
440: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
441: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL));
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode MatGetInfo_HT(Mat N, MatInfoType flag, MatInfo *info)
446: {
447: Mat A;
449: PetscFunctionBegin;
450: PetscCall(MatShellGetContext(N, &A));
451: PetscCall(MatGetInfo(A, flag, info));
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: static PetscErrorCode MatFactorGetSolverType_HT(Mat N, MatSolverType *type)
456: {
457: Mat A;
459: PetscFunctionBegin;
460: PetscCall(MatShellGetContext(N, &A));
461: PetscCall(MatFactorGetSolverType(A, type));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
466: {
467: Mat A, C;
469: PetscFunctionBegin;
470: PetscCall(MatShellGetContext(N, &A));
471: PetscCall(MatDuplicate(A, op, &C));
472: PetscCall(MatCreateHermitianTranspose(C, m));
473: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
474: PetscCall(MatDestroy(&C));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has)
479: {
480: Mat A;
482: PetscFunctionBegin;
483: PetscCall(MatShellGetContext(mat, &A));
484: *has = PETSC_FALSE;
485: if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
486: PetscCall(MatHasOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, has));
487: if (!*has) PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
488: } else if (op == MATOP_MULT_HERMITIAN_TRANSPOSE || op == MATOP_MULT_HERMITIAN_TRANS_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
489: PetscCall(MatHasOperation(A, MATOP_MULT, has));
490: } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
495: {
496: PetscFunctionBegin;
497: PetscCall(MatShellGetContext(N, M));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: /*@
502: MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`
504: Logically Collective
506: Input Parameter:
507: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix
509: Output Parameter:
510: . M - the matrix object stored inside A
512: Level: intermediate
514: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
515: @*/
516: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
517: {
518: PetscFunctionBegin;
521: PetscAssertPointer(M, 2);
522: PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v)
527: {
528: Mat A;
530: PetscFunctionBegin;
531: PetscCall(MatShellGetContext(N, &A));
532: PetscCall(MatGetDiagonal(A, v));
533: PetscCall(VecConjugate(v));
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: static PetscErrorCode MatCopy_HT(Mat A, Mat B, MatStructure str)
538: {
539: Mat a, b;
541: PetscFunctionBegin;
542: PetscCall(MatShellGetContext(A, &a));
543: PetscCall(MatShellGetContext(B, &b));
544: PetscCall(MatCopy(a, b, str));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
549: {
550: Mat A;
551: PetscScalar vscale = 1.0, vshift = 0.0;
552: PetscBool flg;
554: PetscFunctionBegin;
555: PetscCall(MatShellGetContext(N, &A));
556: PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg));
557: 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 */
558: 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));
559: }
560: if (flg) {
561: Mat B;
563: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B));
564: if (reuse != MAT_INPLACE_MATRIX) {
565: PetscCall(MatConvert(B, newtype, reuse, newmat));
566: PetscCall(MatDestroy(&B));
567: } else {
568: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
569: PetscCall(MatHeaderReplace(N, &B));
570: }
571: } else { /* use basic converter as fallback */
572: flg = (PetscBool)(N->ops->getrow != NULL);
573: PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
574: }
575: if (flg) {
576: PetscCall(MatScale(*newmat, vscale));
577: PetscCall(MatShift(*newmat, vshift));
578: }
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*MC
583: MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix
585: Level: advanced
587: Developer Notes:
588: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
590: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
592: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`
593: M*/
595: /*@
596: MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'*
598: Collective
600: Input Parameter:
601: . A - the (possibly rectangular) matrix
603: Output Parameter:
604: . N - the matrix that represents A'*
606: Level: intermediate
608: Note:
609: The Hermitian transpose A' is NOT actually formed! Rather the new matrix
610: object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
611: the original matrix
613: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
614: `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
615: @*/
616: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
617: {
618: VecType vtype;
620: PetscFunctionBegin;
621: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
622: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
623: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
624: PetscCall(MatSetType(*N, MATSHELL));
625: PetscCall(MatShellSetContext(*N, A));
626: PetscCall(PetscObjectReference((PetscObject)A));
628: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
629: PetscCall(MatGetVecType(A, &vtype));
630: PetscCall(MatSetVecType(*N, vtype));
631: #if defined(PETSC_HAVE_DEVICE)
632: PetscCall(MatBindToCPU(*N, A->boundtocpu));
633: #endif
634: PetscCall(MatSetUp(*N));
636: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT));
637: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT));
638: PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
639: #if !defined(PETSC_USE_COMPLEX)
640: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
641: #endif
642: PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (void (*)(void))MatLUFactor_HT));
643: PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (void (*)(void))MatCholeskyFactor_HT));
644: PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (void (*)(void))MatGetFactor_HT));
645: PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (void (*)(void))MatGetInfo_HT));
646: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT));
647: PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_HT));
648: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT));
649: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_HT));
650: PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT));
652: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
653: #if !defined(PETSC_USE_COMPLEX)
654: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
655: #endif
656: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT));
657: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_HT));
658: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
659: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
660: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
661: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }