Actual source code: transm.c
2: #include <petsc/private/matimpl.h>
4: typedef struct {
5: Mat A;
6: } Mat_Transpose;
8: PetscErrorCode MatMult_Transpose(Mat N,Vec x,Vec y)
9: {
10: Mat_Transpose *Na = (Mat_Transpose*)N->data;
12: MatMultTranspose(Na->A,x,y);
13: return 0;
14: }
16: PetscErrorCode MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
17: {
18: Mat_Transpose *Na = (Mat_Transpose*)N->data;
20: MatMultTransposeAdd(Na->A,v1,v2,v3);
21: return 0;
22: }
24: PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y)
25: {
26: Mat_Transpose *Na = (Mat_Transpose*)N->data;
28: MatMult(Na->A,x,y);
29: return 0;
30: }
32: PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
33: {
34: Mat_Transpose *Na = (Mat_Transpose*)N->data;
36: MatMultAdd(Na->A,v1,v2,v3);
37: return 0;
38: }
40: PetscErrorCode MatDestroy_Transpose(Mat N)
41: {
42: Mat_Transpose *Na = (Mat_Transpose*)N->data;
44: MatDestroy(&Na->A);
45: PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
46: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
47: PetscFree(N->data);
48: return 0;
49: }
51: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
52: {
53: Mat_Transpose *Na = (Mat_Transpose*)N->data;
55: if (op == MAT_COPY_VALUES) {
56: MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);
57: } else if (op == MAT_DO_NOT_COPY_VALUES) {
58: MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
59: MatTranspose(*m,MAT_INPLACE_MATRIX,m);
60: } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
61: return 0;
62: }
64: PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
65: {
66: Mat_Transpose *Aa = (Mat_Transpose*)A->data;
68: MatCreateVecs(Aa->A,l,r);
69: return 0;
70: }
72: PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)
73: {
74: Mat_Transpose *Ya = (Mat_Transpose*)Y->data;
75: Mat_Transpose *Xa = (Mat_Transpose*)X->data;
76: Mat M = Ya->A;
77: Mat N = Xa->A;
79: MatAXPY(M,a,N,str);
80: return 0;
81: }
83: PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
84: {
85: Mat_Transpose *X = (Mat_Transpose*)mat->data;
87: *has = PETSC_FALSE;
88: if (op == MATOP_MULT) {
89: MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);
90: } else if (op == MATOP_MULT_TRANSPOSE) {
91: MatHasOperation(X->A,MATOP_MULT,has);
92: } else if (op == MATOP_MULT_ADD) {
93: MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);
94: } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
95: MatHasOperation(X->A,MATOP_MULT_ADD,has);
96: } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
97: return 0;
98: }
100: /* used by hermitian transpose */
101: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
102: {
103: Mat A,B,C,Ain,Bin,Cin;
104: PetscBool Aistrans,Bistrans,Cistrans;
105: PetscInt Atrans,Btrans,Ctrans;
106: MatProductType ptype;
108: MatCheckProduct(D,1);
109: A = D->product->A;
110: B = D->product->B;
111: C = D->product->C;
112: PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);
113: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);
114: PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);
116: Atrans = 0;
117: Ain = A;
118: while (Aistrans) {
119: Atrans++;
120: MatTransposeGetMat(Ain,&Ain);
121: PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);
122: }
123: Btrans = 0;
124: Bin = B;
125: while (Bistrans) {
126: Btrans++;
127: MatTransposeGetMat(Bin,&Bin);
128: PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);
129: }
130: Ctrans = 0;
131: Cin = C;
132: while (Cistrans) {
133: Ctrans++;
134: MatTransposeGetMat(Cin,&Cin);
135: PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&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) Atrans = 0;
142: if (Bin->symmetric) Btrans = 0;
143: if (Cin && Cin->symmetric) Ctrans = 0;
145: if (Atrans || Btrans || Ctrans) {
146: ptype = MATPRODUCT_UNSPECIFIED;
147: switch (D->product->type) {
148: case MATPRODUCT_AB:
149: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
150: /* TODO custom implementation ? */
151: } else if (Atrans) { /* At * B */
152: ptype = MATPRODUCT_AtB;
153: } else { /* A * Bt */
154: ptype = MATPRODUCT_ABt;
155: }
156: break;
157: case MATPRODUCT_AtB:
158: if (Atrans && Btrans) { /* A * Bt */
159: ptype = MATPRODUCT_ABt;
160: } else if (Atrans) { /* A * B */
161: ptype = MATPRODUCT_AB;
162: } else { /* At * Bt we do not have support for this */
163: /* TODO custom implementation ? */
164: }
165: break;
166: case MATPRODUCT_ABt:
167: if (Atrans && Btrans) { /* At * B */
168: ptype = MATPRODUCT_AtB;
169: } else if (Atrans) { /* At * Bt we do not have support for this */
170: /* TODO custom implementation ? */
171: } else { /* A * B */
172: ptype = MATPRODUCT_AB;
173: }
174: break;
175: case MATPRODUCT_PtAP:
176: if (Atrans) { /* PtAtP */
177: /* TODO custom implementation ? */
178: } else { /* RARt */
179: ptype = MATPRODUCT_RARt;
180: }
181: break;
182: case MATPRODUCT_RARt:
183: if (Atrans) { /* RAtRt */
184: /* TODO custom implementation ? */
185: } else { /* PtAP */
186: ptype = MATPRODUCT_PtAP;
187: }
188: break;
189: case MATPRODUCT_ABC:
190: /* TODO custom implementation ? */
191: break;
192: default: SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[D->product->type]);
193: }
194: }
195: MatProductReplaceMats(Ain,Bin,Cin,D);
196: MatProductSetType(D,ptype);
197: MatProductSetFromOptions(D);
198: return 0;
199: }
201: PetscErrorCode MatGetDiagonal_Transpose(Mat A,Vec v)
202: {
203: Mat_Transpose *Aa = (Mat_Transpose*)A->data;
205: MatGetDiagonal(Aa->A,v);
206: return 0;
207: }
209: PetscErrorCode MatConvert_Transpose(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
210: {
211: Mat_Transpose *Aa = (Mat_Transpose*)A->data;
212: PetscBool flg;
214: MatHasOperation(Aa->A,MATOP_TRANSPOSE,&flg);
215: if (flg) {
216: Mat B;
218: MatTranspose(Aa->A,MAT_INITIAL_MATRIX,&B);
219: if (reuse != MAT_INPLACE_MATRIX) {
220: MatConvert(B,newtype,reuse,newmat);
221: MatDestroy(&B);
222: } else {
223: MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
224: MatHeaderReplace(A,&B);
225: }
226: } else { /* use basic converter as fallback */
227: MatConvert_Basic(A,newtype,reuse,newmat);
228: }
229: return 0;
230: }
232: PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
233: {
234: Mat_Transpose *Aa = (Mat_Transpose*)A->data;
236: *M = Aa->A;
237: return 0;
238: }
240: /*@
241: MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT
243: Logically collective on Mat
245: Input Parameter:
246: . A - the MATTRANSPOSE matrix
248: Output Parameter:
249: . M - the matrix object stored inside A
251: Level: intermediate
253: .seealso: MatCreateTranspose()
255: @*/
256: PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
257: {
261: PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));
262: return 0;
263: }
265: /*@
266: MatCreateTranspose - Creates a new matrix object that behaves like A'
268: Collective on Mat
270: Input Parameter:
271: . A - the (possibly rectangular) matrix
273: Output Parameter:
274: . N - the matrix that represents A'
276: Level: intermediate
278: Notes:
279: The transpose A' is NOT actually formed! Rather the new matrix
280: object performs the matrix-vector product by using the MatMultTranspose() on
281: the original matrix
283: .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()
285: @*/
286: PetscErrorCode MatCreateTranspose(Mat A,Mat *N)
287: {
288: PetscInt m,n;
289: Mat_Transpose *Na;
290: VecType vtype;
292: MatGetLocalSize(A,&m,&n);
293: MatCreate(PetscObjectComm((PetscObject)A),N);
294: MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
295: PetscLayoutSetUp((*N)->rmap);
296: PetscLayoutSetUp((*N)->cmap);
297: PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);
299: PetscNewLog(*N,&Na);
300: (*N)->data = (void*) Na;
301: PetscObjectReference((PetscObject)A);
302: Na->A = A;
304: (*N)->ops->destroy = MatDestroy_Transpose;
305: (*N)->ops->mult = MatMult_Transpose;
306: (*N)->ops->multadd = MatMultAdd_Transpose;
307: (*N)->ops->multtranspose = MatMultTranspose_Transpose;
308: (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose;
309: (*N)->ops->duplicate = MatDuplicate_Transpose;
310: (*N)->ops->getvecs = MatCreateVecs_Transpose;
311: (*N)->ops->axpy = MatAXPY_Transpose;
312: (*N)->ops->hasoperation = MatHasOperation_Transpose;
313: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
314: (*N)->ops->getdiagonal = MatGetDiagonal_Transpose;
315: (*N)->ops->convert = MatConvert_Transpose;
316: (*N)->assembled = PETSC_TRUE;
318: PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);
319: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
320: MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
321: MatGetVecType(A,&vtype);
322: MatSetVecType(*N,vtype);
323: #if defined(PETSC_HAVE_DEVICE)
324: MatBindToCPU(*N,A->boundtocpu);
325: #endif
326: MatSetUp(*N);
327: return 0;
328: }