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: PetscContainer container;
318: void *data;
319: } MatProductCtx_Transpose;
321: static PetscErrorCode MatProductCtxDestroy_Transpose(void **ptr)
322: {
323: MatProductCtx_Transpose *data = *(MatProductCtx_Transpose **)ptr;
325: PetscFunctionBegin;
326: if (data->data) PetscCall((*data->destroy)(&data->data));
327: PetscCall(PetscContainerDestroy(&data->container));
328: PetscCall(PetscFree(data));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode MatProductNumeric_Transpose(Mat D)
333: {
334: Mat_Product *product;
335: MatProductCtx_Transpose *data;
336: PetscContainer container;
338: PetscFunctionBegin;
339: MatCheckProduct(D, 1);
340: PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
341: product = D->product;
342: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
343: PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_Transpose missing");
344: PetscCall(PetscContainerGetPointer(container, (void **)&data));
345: data = (MatProductCtx_Transpose *)product->data;
346: product->data = data->data;
347: PetscCall((*data->numeric)(D));
348: PetscCall(MatScale(D, data->scale));
349: product->data = data;
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: static PetscErrorCode MatProductSymbolic_Transpose(Mat D)
354: {
355: Mat_Product *product;
356: MatProductCtx_Transpose *data;
357: PetscContainer container;
359: PetscFunctionBegin;
360: MatCheckProduct(D, 1);
361: product = D->product;
362: if (D->ops->productsymbolic == MatProductSymbolic_Transpose) {
363: PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
364: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
365: PetscCheck(container, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatProductCtx_Transpose missing");
366: PetscCall(PetscContainerGetPointer(container, (void **)&data));
367: PetscCall(MatProductSetFromOptions(D));
368: PetscCall(MatProductSymbolic(D));
369: data->numeric = D->ops->productnumeric;
370: data->destroy = product->destroy;
371: data->data = product->data;
372: D->ops->productnumeric = MatProductNumeric_Transpose;
373: product->destroy = MatProductCtxDestroy_Transpose;
374: product->data = data;
375: }
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: static PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
380: {
381: Mat A, B, C, Ain, Bin, Cin;
382: PetscScalar scale = 1.0, vscale;
383: PetscBool Aistrans, Bistrans, Cistrans;
384: PetscInt Atrans, Btrans, Ctrans;
385: PetscContainer container = NULL;
386: MatProductCtx_Transpose *data;
387: MatProductType ptype;
389: PetscFunctionBegin;
390: MatCheckProduct(D, 1);
391: A = D->product->A;
392: B = D->product->B;
393: C = D->product->C;
394: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
395: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
396: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
397: PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
398: Atrans = 0;
399: Ain = A;
400: while (Aistrans) {
401: Atrans++;
402: 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));
403: scale *= vscale;
404: PetscCall(MatTransposeGetMat(Ain, &Ain));
405: PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
406: }
407: Btrans = 0;
408: Bin = B;
409: while (Bistrans) {
410: Btrans++;
411: 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));
412: scale *= vscale;
413: PetscCall(MatTransposeGetMat(Bin, &Bin));
414: PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
415: }
416: Ctrans = 0;
417: Cin = C;
418: while (Cistrans) {
419: Ctrans++;
420: 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));
421: scale *= vscale;
422: PetscCall(MatTransposeGetMat(Cin, &Cin));
423: PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
424: }
425: Atrans = Atrans % 2;
426: Btrans = Btrans % 2;
427: Ctrans = Ctrans % 2;
428: ptype = D->product->type; /* same product type by default */
429: if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
430: if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
431: if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;
433: if (Atrans || Btrans || Ctrans) {
434: if (scale != 1.0) {
435: PetscCall(PetscObjectQuery((PetscObject)D, "MatProductCtx_Transpose", (PetscObject *)&container));
436: if (!container) {
437: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)D), &container));
438: PetscCall(PetscNew(&data));
439: PetscCall(PetscContainerSetPointer(container, data));
440: PetscCall(PetscObjectCompose((PetscObject)D, "MatProductCtx_Transpose", (PetscObject)container));
441: } else PetscCall(PetscContainerGetPointer(container, (void **)&data));
442: data->scale = scale;
443: data->container = container;
444: }
445: ptype = MATPRODUCT_UNSPECIFIED;
446: switch (D->product->type) {
447: case MATPRODUCT_AB:
448: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
449: /* TODO custom implementation ? */
450: } else if (Atrans) { /* At * B */
451: ptype = MATPRODUCT_AtB;
452: } else { /* A * Bt */
453: ptype = MATPRODUCT_ABt;
454: }
455: break;
456: case MATPRODUCT_AtB:
457: if (Atrans && Btrans) { /* A * Bt */
458: ptype = MATPRODUCT_ABt;
459: } else if (Atrans) { /* A * B */
460: ptype = MATPRODUCT_AB;
461: } else { /* At * Bt we do not have support for this */
462: /* TODO custom implementation ? */
463: }
464: break;
465: case MATPRODUCT_ABt:
466: if (Atrans && Btrans) { /* At * B */
467: ptype = MATPRODUCT_AtB;
468: } else if (Atrans) { /* At * Bt we do not have support for this */
469: /* TODO custom implementation ? */
470: } else { /* A * B */
471: ptype = MATPRODUCT_AB;
472: }
473: break;
474: case MATPRODUCT_PtAP:
475: if (Atrans) { /* PtAtP */
476: /* TODO custom implementation ? */
477: } else { /* RARt */
478: ptype = MATPRODUCT_RARt;
479: }
480: break;
481: case MATPRODUCT_RARt:
482: if (Atrans) { /* RAtRt */
483: /* TODO custom implementation ? */
484: } else { /* PtAP */
485: ptype = MATPRODUCT_PtAP;
486: }
487: break;
488: case MATPRODUCT_ABC:
489: /* TODO custom implementation ? */
490: break;
491: default:
492: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
493: }
494: }
495: PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
496: PetscCall(MatProductSetType(D, ptype));
497: if (container) D->ops->productsymbolic = MatProductSymbolic_Transpose;
498: else PetscCall(MatProductSetFromOptions(D));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
503: {
504: Mat A;
506: PetscFunctionBegin;
507: PetscCall(MatShellGetContext(N, &A));
508: PetscCall(MatGetDiagonal(A, v));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
513: {
514: Mat a, b;
516: PetscFunctionBegin;
517: PetscCall(MatShellGetContext(A, &a));
518: PetscCall(MatShellGetContext(B, &b));
519: PetscCall(MatCopy(a, b, str));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
524: {
525: Mat A;
526: PetscScalar vscale = 1.0, vshift = 0.0;
527: PetscBool flg;
529: PetscFunctionBegin;
530: PetscCall(MatShellGetContext(N, &A));
531: PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
532: 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 */
533: 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));
534: }
535: if (flg) {
536: Mat B;
538: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
539: if (reuse != MAT_INPLACE_MATRIX) {
540: PetscCall(MatConvert(B, newtype, reuse, newmat));
541: PetscCall(MatDestroy(&B));
542: } else {
543: PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
544: PetscCall(MatHeaderReplace(N, &B));
545: }
546: } else { /* use basic converter as fallback */
547: flg = (PetscBool)(N->ops->getrow != NULL);
548: PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
549: }
550: if (flg) {
551: PetscCall(MatScale(*newmat, vscale));
552: PetscCall(MatShift(*newmat, vshift));
553: }
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
558: {
559: PetscFunctionBegin;
560: PetscCall(MatShellGetContext(N, M));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: /*@
565: MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`
567: Logically Collective
569: Input Parameter:
570: . A - the `MATTRANSPOSEVIRTUAL` matrix
572: Output Parameter:
573: . M - the matrix object stored inside `A`
575: Level: intermediate
577: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
578: @*/
579: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
580: {
581: PetscFunctionBegin;
584: PetscAssertPointer(M, 2);
585: PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*MC
590: MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix
592: Level: advanced
594: Developer Notes:
595: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
597: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
599: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
600: `MATNORMALHERMITIAN`, `MATNORMAL`
601: M*/
603: /*@
604: MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'
606: Collective
608: Input Parameter:
609: . A - the (possibly rectangular) matrix
611: Output Parameter:
612: . N - the matrix that represents A'
614: Level: intermediate
616: Note:
617: The transpose A' is NOT actually formed! Rather the new matrix
618: object performs the matrix-vector product by using the `MatMultTranspose()` on
619: the original matrix
621: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
622: `MATNORMALHERMITIAN`
623: @*/
624: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
625: {
626: VecType vtype;
628: PetscFunctionBegin;
629: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
630: PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
631: PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
632: PetscCall(MatSetType(*N, MATSHELL));
633: PetscCall(MatShellSetContext(*N, A));
634: PetscCall(PetscObjectReference((PetscObject)A));
636: PetscCall(MatSetBlockSizes(*N, A->cmap->bs, A->rmap->bs));
637: PetscCall(MatGetVecType(A, &vtype));
638: PetscCall(MatSetVecType(*N, vtype));
639: #if defined(PETSC_HAVE_DEVICE)
640: PetscCall(MatBindToCPU(*N, A->boundtocpu));
641: #endif
642: PetscCall(MatSetUp(*N));
644: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Transpose));
645: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Transpose));
646: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Transpose));
647: PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (PetscErrorCodeFn *)MatLUFactor_Transpose));
648: PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (PetscErrorCodeFn *)MatCholeskyFactor_Transpose));
649: PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (PetscErrorCodeFn *)MatGetFactor_Transpose));
650: PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (PetscErrorCodeFn *)MatGetInfo_Transpose));
651: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (PetscErrorCodeFn *)MatDuplicate_Transpose));
652: PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (PetscErrorCodeFn *)MatHasOperation_Transpose));
653: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_Transpose));
654: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_Transpose));
655: PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (PetscErrorCodeFn *)MatConvert_Transpose));
657: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
658: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
659: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_Transpose));
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, MATTRANSPOSEVIRTUAL));
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }