Actual source code: transm.c
1: #include <../src/mat/impls/shell/shell.h>
3: static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y)
4: {
5: Mat A;
7: PetscFunctionBegin;
8: PetscCall(MatShellGetContext(N, &A));
9: PetscCall(MatMultTranspose(A, x, y));
10: PetscFunctionReturn(PETSC_SUCCESS);
11: }
13: static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y)
14: {
15: Mat A;
17: PetscFunctionBegin;
18: PetscCall(MatShellGetContext(N, &A));
19: PetscCall(MatMult(A, x, y));
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: static PetscErrorCode MatSolve_Transpose_LU(Mat N, Vec b, Vec x)
24: {
25: Mat A;
27: PetscFunctionBegin;
28: PetscCall(MatShellGetContext(N, &A));
29: PetscCall(MatSolveTranspose(A, b, x));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: static PetscErrorCode MatSolveAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x)
34: {
35: Mat A;
37: PetscFunctionBegin;
38: PetscCall(MatShellGetContext(N, &A));
39: PetscCall(MatSolveTransposeAdd(A, b, y, x));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode MatSolveTranspose_Transpose_LU(Mat N, Vec b, Vec x)
44: {
45: Mat A;
47: PetscFunctionBegin;
48: PetscCall(MatShellGetContext(N, &A));
49: PetscCall(MatSolve(A, b, x));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode MatSolveTransposeAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x)
54: {
55: Mat A;
57: PetscFunctionBegin;
58: PetscCall(MatShellGetContext(N, &A));
59: PetscCall(MatSolveAdd(A, b, y, x));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode MatMatSolve_Transpose_LU(Mat N, Mat B, Mat X)
64: {
65: Mat A;
67: PetscFunctionBegin;
68: PetscCall(MatShellGetContext(N, &A));
69: PetscCall(MatMatSolveTranspose(A, B, X));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode MatMatSolveTranspose_Transpose_LU(Mat N, Mat B, Mat X)
74: {
75: Mat A;
77: PetscFunctionBegin;
78: PetscCall(MatShellGetContext(N, &A));
79: PetscCall(MatMatSolve(A, B, X));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode MatLUFactor_Transpose(Mat N, IS row, IS col, const MatFactorInfo *minfo)
84: {
85: Mat A;
87: PetscFunctionBegin;
88: PetscCall(MatShellGetContext(N, &A));
89: PetscCall(MatLUFactor(A, col, row, minfo));
90: PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_LU));
91: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_LU));
92: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_LU));
93: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_LU));
94: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_LU));
95: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_LU));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode MatSolve_Transpose_Cholesky(Mat N, Vec b, Vec x)
100: {
101: Mat A;
103: PetscFunctionBegin;
104: PetscCall(MatShellGetContext(N, &A));
105: PetscCall(MatSolveTranspose(A, b, x));
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: static PetscErrorCode MatSolveAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x)
110: {
111: Mat A;
113: PetscFunctionBegin;
114: PetscCall(MatShellGetContext(N, &A));
115: PetscCall(MatSolveTransposeAdd(A, b, y, x));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode MatSolveTranspose_Transpose_Cholesky(Mat N, Vec b, Vec x)
120: {
121: Mat A;
123: PetscFunctionBegin;
124: PetscCall(MatShellGetContext(N, &A));
125: PetscCall(MatSolve(A, b, x));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode MatSolveTransposeAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x)
130: {
131: Mat A;
133: PetscFunctionBegin;
134: PetscCall(MatShellGetContext(N, &A));
135: PetscCall(MatSolveAdd(A, b, y, x));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: static PetscErrorCode MatMatSolve_Transpose_Cholesky(Mat N, Mat B, Mat X)
140: {
141: Mat A;
143: PetscFunctionBegin;
144: PetscCall(MatShellGetContext(N, &A));
145: PetscCall(MatMatSolveTranspose(A, B, X));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode MatMatSolveTranspose_Transpose_Cholesky(Mat N, Mat B, Mat X)
150: {
151: Mat A;
153: PetscFunctionBegin;
154: PetscCall(MatShellGetContext(N, &A));
155: PetscCall(MatMatSolve(A, B, X));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode MatCholeskyFactor_Transpose(Mat N, IS perm, const MatFactorInfo *minfo)
160: {
161: Mat A;
163: PetscFunctionBegin;
164: PetscCall(MatShellGetContext(N, &A));
165: PetscCall(MatCholeskyFactor(A, perm, minfo));
166: PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_Cholesky));
167: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_Cholesky));
168: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_Cholesky));
169: PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_Cholesky));
170: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_Cholesky));
171: PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_Cholesky));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode MatLUFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info)
176: {
177: Mat A, FA;
179: PetscFunctionBegin;
180: PetscCall(MatShellGetContext(N, &A));
181: PetscCall(MatShellGetContext(F, &FA));
182: PetscCall(MatLUFactorNumeric(FA, A, info));
183: PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_LU));
184: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_LU));
185: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_LU));
186: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_LU));
187: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_LU));
188: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_LU));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode MatLUFactorSymbolic_Transpose(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info)
193: {
194: Mat A, FA;
196: PetscFunctionBegin;
197: PetscCall(MatShellGetContext(N, &A));
198: PetscCall(MatShellGetContext(F, &FA));
199: PetscCall(MatLUFactorSymbolic(FA, A, row, col, info));
200: PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (PetscErrorCodeFn *)MatLUFactorNumeric_Transpose));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode MatCholeskyFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info)
205: {
206: Mat A, FA;
208: PetscFunctionBegin;
209: PetscCall(MatShellGetContext(N, &A));
210: PetscCall(MatShellGetContext(F, &FA));
211: PetscCall(MatCholeskyFactorNumeric(FA, A, info));
212: PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (PetscErrorCodeFn *)MatSolve_Transpose_Cholesky));
213: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (PetscErrorCodeFn *)MatSolveAdd_Transpose_Cholesky));
214: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatSolveTranspose_Transpose_Cholesky));
215: PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (PetscErrorCodeFn *)MatSolveTransposeAdd_Transpose_Cholesky));
216: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (PetscErrorCodeFn *)MatMatSolve_Transpose_Cholesky));
217: PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (PetscErrorCodeFn *)MatMatSolveTranspose_Transpose_Cholesky));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode MatCholeskyFactorSymbolic_Transpose(Mat F, Mat N, IS perm, const MatFactorInfo *info)
222: {
223: Mat A, FA;
225: PetscFunctionBegin;
226: PetscCall(MatShellGetContext(N, &A));
227: PetscCall(MatShellGetContext(F, &FA));
228: PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info));
229: PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (PetscErrorCodeFn *)MatCholeskyFactorNumeric_Transpose));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode MatGetFactor_Transpose(Mat N, MatSolverType type, MatFactorType ftype, Mat *F)
234: {
235: Mat A, FA;
237: PetscFunctionBegin;
238: PetscCall(MatShellGetContext(N, &A));
239: PetscCall(MatGetFactor(A, type, ftype, &FA));
240: PetscCall(MatCreateTranspose(FA, F));
241: if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatLUFactorSymbolic_Transpose));
242: else if (ftype == MAT_FACTOR_CHOLESKY) {
243: PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (PetscErrorCodeFn *)MatCholeskyFactorSymbolic_Transpose));
244: PetscCall(MatPropagateSymmetryOptions(A, FA));
245: } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]);
246: (*F)->factortype = ftype;
247: PetscCall(MatDestroy(&FA));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode MatDestroy_Transpose(Mat N)
252: {
253: Mat A;
255: PetscFunctionBegin;
256: PetscCall(MatShellGetContext(N, &A));
257: PetscCall(MatDestroy(&A));
258: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
259: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
260: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
261: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL));
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: static PetscErrorCode MatGetInfo_Transpose(Mat N, MatInfoType flag, MatInfo *info)
266: {
267: Mat A;
269: PetscFunctionBegin;
270: PetscCall(MatShellGetContext(N, &A));
271: PetscCall(MatGetInfo(A, flag, info));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode MatFactorGetSolverType_Transpose(Mat N, MatSolverType *type)
276: {
277: Mat A;
279: PetscFunctionBegin;
280: PetscCall(MatShellGetContext(N, &A));
281: PetscCall(MatFactorGetSolverType(A, type));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
286: {
287: Mat A, C;
289: PetscFunctionBegin;
290: PetscCall(MatShellGetContext(N, &A));
291: PetscCall(MatDuplicate(A, op, &C));
292: PetscCall(MatCreateTranspose(C, m));
293: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
294: PetscCall(MatDestroy(&C));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
299: {
300: Mat A;
302: PetscFunctionBegin;
303: PetscCall(MatShellGetContext(mat, &A));
304: *has = PETSC_FALSE;
305: if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
306: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
307: } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
308: PetscCall(MatHasOperation(A, MATOP_MULT, has));
309: } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: typedef struct {
314: PetscErrorCode (*numeric)(Mat);
315: PetscCtxDestroyFn *destroy;
316: PetscScalar scale;
317: Mat D;
318: void *data;
319: } MatProductCtx_Transpose;
321: static PetscErrorCode MatProductCtxDestroy_Transpose(PetscCtxRt ptr)
322: {
323: MatProductCtx_Transpose *data = *(MatProductCtx_Transpose **)ptr;
324: PetscContainer container;
326: PetscFunctionBegin;
327: if (data->data) PetscCall((*data->destroy)(&data->data));
328: PetscCall(PetscObjectQuery((PetscObject)data->D, "MatProductCtx_Transpose", (PetscObject *)&container));
329: PetscCall(PetscContainerDestroy(&container));
330: PetscCall(PetscObjectCompose((PetscObject)data->D, "MatProductCtx_Transpose", NULL));
331: PetscCall(PetscFree(data));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: static PetscErrorCode MatProductNumeric_Transpose(Mat D)
336: {
337: Mat_Product *product;
338: MatProductCtx_Transpose *data;
339: PetscContainer container;
341: PetscFunctionBegin;
342: MatCheckProduct(D, 1);
343: PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
344: product = D->product;
345: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
346: PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_Transpose missing");
347: PetscCall(PetscContainerGetPointer(container, &data));
348: data = (MatProductCtx_Transpose *)product->data;
349: product->data = data->data;
350: PetscCall((*data->numeric)(D));
351: PetscCall(MatScale(D, data->scale));
352: product->data = data;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode MatProductSymbolic_Transpose(Mat D)
357: {
358: Mat_Product *product;
359: MatProductCtx_Transpose *data;
360: PetscContainer container;
362: PetscFunctionBegin;
363: MatCheckProduct(D, 1);
364: product = D->product;
365: if (D->ops->productsymbolic == MatProductSymbolic_Transpose) {
366: PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
367: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
368: PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_Transpose missing");
369: PetscCall(PetscContainerGetPointer(container, &data));
370: PetscCall(MatProductSetFromOptions(D));
371: PetscCall(MatProductSymbolic(D));
372: data->numeric = D->ops->productnumeric;
373: data->destroy = product->destroy;
374: data->data = product->data;
375: D->ops->productnumeric = MatProductNumeric_Transpose;
376: product->destroy = MatProductCtxDestroy_Transpose;
377: product->data = data;
378: }
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
383: {
384: Mat A, B, C, Ain, Bin, Cin;
385: PetscScalar scale = 1.0, vscale;
386: PetscBool Aistrans, Bistrans, Cistrans;
387: PetscInt Atrans, Btrans, Ctrans;
388: PetscContainer container = NULL;
389: MatProductCtx_Transpose *data;
390: MatProductType ptype;
392: PetscFunctionBegin;
393: MatCheckProduct(D, 1);
394: A = D->product->A;
395: B = D->product->B;
396: C = D->product->C;
397: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
398: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
399: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
400: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
401: Atrans = 0;
402: Ain = A;
403: while (Aistrans) {
404: Atrans++;
405: 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));
406: scale *= vscale;
407: PetscCall(MatTransposeGetMat(Ain, &Ain));
408: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
409: }
410: Btrans = 0;
411: Bin = B;
412: while (Bistrans) {
413: Btrans++;
414: 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));
415: scale *= vscale;
416: PetscCall(MatTransposeGetMat(Bin, &Bin));
417: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
418: }
419: Ctrans = 0;
420: Cin = C;
421: while (Cistrans) {
422: Ctrans++;
423: 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));
424: scale *= vscale;
425: PetscCall(MatTransposeGetMat(Cin, &Cin));
426: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
427: }
428: Atrans = Atrans % 2;
429: Btrans = Btrans % 2;
430: Ctrans = Ctrans % 2;
431: ptype = D->product->type; /* same product type by default */
432: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
433: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
434: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
436: if (Atrans || Btrans || Ctrans) {
437: if (scale != 1.0) {
438: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
439: if (!container) {
440: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)D), &container));
441: PetscCall(PetscNew(&data));
442: PetscCall(PetscContainerSetPointer(container, data));
443: PetscCall(PetscObjectCompose((PetscObject)D, "MatProductCtx_Transpose", (PetscObject)container));
444: } else PetscCall(PetscContainerGetPointer(container, &data));
445: data->scale = scale;
446: data->D = D;
447: }
448: ptype = MATPRODUCT_UNSPECIFIED;
449: switch (D->product->type) {
450: case MATPRODUCT_AB:
451: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
452: /* TODO custom implementation ? */
453: } else if (Atrans) { /* At * B */
454: ptype = MATPRODUCT_AtB;
455: } else { /* A * Bt */
456: ptype = MATPRODUCT_ABt;
457: }
458: break;
459: case MATPRODUCT_AtB:
460: if (Atrans && Btrans) { /* A * Bt */
461: ptype = MATPRODUCT_ABt;
462: } else if (Atrans) { /* A * B */
463: ptype = MATPRODUCT_AB;
464: } else { /* At * Bt we do not have support for this */
465: /* TODO custom implementation ? */
466: }
467: break;
468: case MATPRODUCT_ABt:
469: if (Atrans && Btrans) { /* At * B */
470: ptype = MATPRODUCT_AtB;
471: } else if (Atrans) { /* At * Bt we do not have support for this */
472: /* TODO custom implementation ? */
473: } else { /* A * B */
474: ptype = MATPRODUCT_AB;
475: }
476: break;
477: case MATPRODUCT_PtAP:
478: if (Atrans) { /* PtAtP */
479: /* TODO custom implementation ? */
480: } else { /* RARt */
481: ptype = MATPRODUCT_RARt;
482: }
483: break;
484: case MATPRODUCT_RARt:
485: if (Atrans) { /* RAtRt */
486: /* TODO custom implementation ? */
487: } else { /* PtAP */
488: ptype = MATPRODUCT_PtAP;
489: }
490: break;
491: case MATPRODUCT_ABC:
492: /* TODO custom implementation ? */
493: break;
494: default:
495: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
496: }
497: }
498: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
499: PetscCall(MatProductSetType(D, ptype));
500: if (container) D->ops->productsymbolic = MatProductSymbolic_Transpose;
501: else PetscCall(MatProductSetFromOptions(D));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
506: {
507: Mat A;
509: PetscFunctionBegin;
510: PetscCall(MatShellGetContext(N, &A));
511: PetscCall(MatGetDiagonal(A, v));
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
516: {
517: Mat a, b;
519: PetscFunctionBegin;
520: PetscCall(MatShellGetContext(A, &a));
521: PetscCall(MatShellGetContext(B, &b));
522: PetscCall(MatCopy(a, b, str));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
527: {
528: Mat A;
529: PetscScalar vscale = 1.0, vshift = 0.0;
530: PetscBool flg;
532: PetscFunctionBegin;
533: PetscCall(MatShellGetContext(N, &A));
534: PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
535: 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 */
536: 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));
537: }
538: if (flg) {
539: Mat B;
541: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
542: if (reuse != MAT_INPLACE_MATRIX) {
543: PetscCall(MatConvert(B, newtype, reuse, newmat));
544: PetscCall(MatDestroy(&B));
545: } else {
546: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
547: PetscCall(MatHeaderReplace(N, &B));
548: }
549: } else { /* use basic converter as fallback */
550: flg = (PetscBool)(N->ops->getrow != NULL);
551: PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
552: }
553: if (flg) {
554: PetscCall(MatScale(*newmat, vscale));
555: PetscCall(MatShift(*newmat, vshift));
556: }
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
561: {
562: PetscFunctionBegin;
563: PetscCall(MatShellGetContext(N, M));
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*@
568: MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`
570: Logically Collective
572: Input Parameter:
573: . A - the `MATTRANSPOSEVIRTUAL` matrix
575: Output Parameter:
576: . M - the matrix object stored inside `A`
578: Level: intermediate
580: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
581: @*/
582: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
583: {
584: PetscFunctionBegin;
587: PetscAssertPointer(M, 2);
588: PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*MC
593: MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix
595: Level: advanced
597: Developer Notes:
598: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
600: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
602: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
603: `MATNORMALHERMITIAN`, `MATNORMAL`
604: M*/
606: /*@
607: MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'
609: Collective
611: Input Parameter:
612: . A - the (possibly rectangular) matrix
614: Output Parameter:
615: . N - the matrix that represents A'
617: Level: intermediate
619: Note:
620: The transpose A' is NOT actually formed! Rather the new matrix
621: object performs the matrix-vector product by using the `MatMultTranspose()` on
622: the original matrix
624: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
625: `MATNORMALHERMITIAN`
626: @*/
627: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
628: {
629: VecType vtype;
631: PetscFunctionBegin;
632: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
633: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
634: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
635: PetscCall(MatSetType(*N, MATSHELL));
636: PetscCall(MatShellSetContext(*N, A));
637: PetscCall(PetscObjectReference((PetscObject)A));
639: PetscCall(MatSetBlockSizes(*N, A->cmap->bs, A->rmap->bs));
640: PetscCall(MatGetVecType(A, &vtype));
641: PetscCall(MatSetVecType(*N, vtype));
642: #if defined(PETSC_HAVE_DEVICE)
643: PetscCall(MatBindToCPU(*N, A->boundtocpu));
644: #endif
645: PetscCall(MatSetUp(*N));
647: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Transpose));
648: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Transpose));
649: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Transpose));
650: PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (PetscErrorCodeFn *)MatLUFactor_Transpose));
651: PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (PetscErrorCodeFn *)MatCholeskyFactor_Transpose));
652: PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (PetscErrorCodeFn *)MatGetFactor_Transpose));
653: PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (PetscErrorCodeFn *)MatGetInfo_Transpose));
654: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_Transpose));
655: PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (PetscErrorCodeFn *)MatHasOperation_Transpose));
656: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Transpose));
657: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_Transpose));
658: PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (PetscErrorCodeFn *)MatConvert_Transpose));
660: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
661: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
662: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_Transpose));
663: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
664: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
665: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
666: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }