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