Actual source code: transm.c
petsc-3.14.6 2021-03-30
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;
14: MatMultTranspose(Na->A,x,y);
15: return(0);
16: }
18: PetscErrorCode MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
19: {
20: Mat_Transpose *Na = (Mat_Transpose*)N->data;
24: MatMultTransposeAdd(Na->A,v1,v2,v3);
25: return(0);
26: }
28: PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y)
29: {
30: Mat_Transpose *Na = (Mat_Transpose*)N->data;
34: MatMult(Na->A,x,y);
35: return(0);
36: }
38: PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
39: {
40: Mat_Transpose *Na = (Mat_Transpose*)N->data;
44: MatMultAdd(Na->A,v1,v2,v3);
45: return(0);
46: }
48: PetscErrorCode MatDestroy_Transpose(Mat N)
49: {
50: Mat_Transpose *Na = (Mat_Transpose*)N->data;
54: MatDestroy(&Na->A);
55: PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
56: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
57: PetscFree(N->data);
58: return(0);
59: }
61: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
62: {
63: Mat_Transpose *Na = (Mat_Transpose*)N->data;
67: if (op == MAT_COPY_VALUES) {
68: MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);
69: } else if (op == MAT_DO_NOT_COPY_VALUES) {
70: MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
71: MatTranspose(*m,MAT_INPLACE_MATRIX,m);
72: } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
73: return(0);
74: }
76: PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
77: {
78: Mat_Transpose *Aa = (Mat_Transpose*)A->data;
82: MatCreateVecs(Aa->A,l,r);
83: return(0);
84: }
86: PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)
87: {
88: Mat_Transpose *Ya = (Mat_Transpose*)Y->data;
89: Mat_Transpose *Xa = (Mat_Transpose*)X->data;
90: Mat M = Ya->A;
91: Mat N = Xa->A;
95: MatAXPY(M,a,N,str);
96: return(0);
97: }
99: PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
100: {
101: Mat_Transpose *X = (Mat_Transpose*)mat->data;
105: *has = PETSC_FALSE;
106: if (op == MATOP_MULT) {
107: MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);
108: } else if (op == MATOP_MULT_TRANSPOSE) {
109: MatHasOperation(X->A,MATOP_MULT,has);
110: } else if (op == MATOP_MULT_ADD) {
111: MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);
112: } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
113: MatHasOperation(X->A,MATOP_MULT_ADD,has);
114: } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
115: return(0);
116: }
118: /* used by hermitian transpose */
119: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
120: {
121: Mat A,B,C,Ain,Bin,Cin;
122: PetscBool Aistrans,Bistrans,Cistrans;
123: PetscInt Atrans,Btrans,Ctrans;
124: MatProductType ptype;
128: MatCheckProduct(D,1);
129: A = D->product->A;
130: B = D->product->B;
131: C = D->product->C;
132: PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);
133: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);
134: PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);
135: if (!Aistrans && !Bistrans && !Cistrans) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"This should not happen");
136: Atrans = 0;
137: Ain = A;
138: while (Aistrans) {
139: Atrans++;
140: MatTransposeGetMat(Ain,&Ain);
141: PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);
142: }
143: Btrans = 0;
144: Bin = B;
145: while (Bistrans) {
146: Btrans++;
147: MatTransposeGetMat(Bin,&Bin);
148: PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);
149: }
150: Ctrans = 0;
151: Cin = C;
152: while (Cistrans) {
153: Ctrans++;
154: MatTransposeGetMat(Cin,&Cin);
155: PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&Cistrans);
156: }
157: Atrans = Atrans%2;
158: Btrans = Btrans%2;
159: Ctrans = Ctrans%2;
160: ptype = D->product->type; /* same product type by default */
161: if (Ain->symmetric) Atrans = 0;
162: if (Bin->symmetric) Btrans = 0;
163: if (Cin && Cin->symmetric) Ctrans = 0;
165: if (Atrans || Btrans || Ctrans) {
166: ptype = MATPRODUCT_UNSPECIFIED;
167: switch (D->product->type) {
168: case MATPRODUCT_AB:
169: if (Atrans && Btrans) { /* At * Bt we do not have support for this */
170: /* TODO custom implementation ? */
171: } else if (Atrans) { /* At * B */
172: ptype = MATPRODUCT_AtB;
173: } else { /* A * Bt */
174: ptype = MATPRODUCT_ABt;
175: }
176: break;
177: case MATPRODUCT_AtB:
178: if (Atrans && Btrans) { /* A * Bt */
179: ptype = MATPRODUCT_ABt;
180: } else if (Atrans) { /* A * B */
181: ptype = MATPRODUCT_AB;
182: } else { /* At * Bt we do not have support for this */
183: /* TODO custom implementation ? */
184: }
185: break;
186: case MATPRODUCT_ABt:
187: if (Atrans && Btrans) { /* At * B */
188: ptype = MATPRODUCT_AtB;
189: } else if (Atrans) { /* At * Bt we do not have support for this */
190: /* TODO custom implementation ? */
191: } else { /* A * B */
192: ptype = MATPRODUCT_AB;
193: }
194: break;
195: case MATPRODUCT_PtAP:
196: if (Atrans) { /* PtAtP */
197: /* TODO custom implementation ? */
198: } else { /* RARt */
199: ptype = MATPRODUCT_RARt;
200: }
201: break;
202: case MATPRODUCT_RARt:
203: if (Atrans) { /* RAtRt */
204: /* TODO custom implementation ? */
205: } else { /* PtAP */
206: ptype = MATPRODUCT_PtAP;
207: }
208: break;
209: case MATPRODUCT_ABC:
210: /* TODO custom implementation ? */
211: break;
212: default: SETERRQ1(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[D->product->type]);
213: }
214: }
215: MatProductReplaceMats(Ain,Bin,Cin,D);
216: MatProductSetType(D,ptype);
217: MatProductSetFromOptions(D);
218: return(0);
219: }
221: PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
222: {
223: Mat_Transpose *Aa = (Mat_Transpose*)A->data;
226: *M = Aa->A;
227: return(0);
228: }
230: /*@
231: MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT
233: Logically collective on Mat
235: Input Parameter:
236: . A - the MATTRANSPOSE matrix
238: Output Parameter:
239: . M - the matrix object stored inside A
241: Level: intermediate
243: .seealso: MatCreateTranspose()
245: @*/
246: PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
247: {
254: PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));
255: return(0);
256: }
258: /*@
259: MatCreateTranspose - Creates a new matrix object that behaves like A'
261: Collective on Mat
263: Input Parameter:
264: . A - the (possibly rectangular) matrix
266: Output Parameter:
267: . N - the matrix that represents A'
269: Level: intermediate
271: Notes:
272: The transpose A' is NOT actually formed! Rather the new matrix
273: object performs the matrix-vector product by using the MatMultTranspose() on
274: the original matrix
276: .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()
278: @*/
279: PetscErrorCode MatCreateTranspose(Mat A,Mat *N)
280: {
282: PetscInt m,n;
283: Mat_Transpose *Na;
284: VecType vtype;
287: MatGetLocalSize(A,&m,&n);
288: MatCreate(PetscObjectComm((PetscObject)A),N);
289: MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
290: PetscLayoutSetUp((*N)->rmap);
291: PetscLayoutSetUp((*N)->cmap);
292: PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);
294: PetscNewLog(*N,&Na);
295: (*N)->data = (void*) Na;
296: PetscObjectReference((PetscObject)A);
297: Na->A = A;
299: (*N)->ops->destroy = MatDestroy_Transpose;
300: (*N)->ops->mult = MatMult_Transpose;
301: (*N)->ops->multadd = MatMultAdd_Transpose;
302: (*N)->ops->multtranspose = MatMultTranspose_Transpose;
303: (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose;
304: (*N)->ops->duplicate = MatDuplicate_Transpose;
305: (*N)->ops->getvecs = MatCreateVecs_Transpose;
306: (*N)->ops->axpy = MatAXPY_Transpose;
307: (*N)->ops->hasoperation = MatHasOperation_Transpose;
308: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
309: (*N)->assembled = PETSC_TRUE;
311: PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);
312: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
313: MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
314: MatGetVecType(A,&vtype);
315: MatSetVecType(*N,vtype);
316: MatSetUp(*N);
317: return(0);
318: }