1
2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
3
4 typedef struct {
5 Mat A;
6 } Mat_Transpose;
7
MatMult_Transpose(Mat N,Vec x,Vec y)8 PetscErrorCode MatMult_Transpose(Mat N,Vec x,Vec y)
9 {
10 Mat_Transpose *Na = (Mat_Transpose*)N->data;
11 PetscErrorCode ierr;
12
13 PetscFunctionBegin;
14 ierr = MatMultTranspose(Na->A,x,y);CHKERRQ(ierr);
15 PetscFunctionReturn(0);
16 }
17
MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)18 PetscErrorCode MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
19 {
20 Mat_Transpose *Na = (Mat_Transpose*)N->data;
21 PetscErrorCode ierr;
22
23 PetscFunctionBegin;
24 ierr = MatMultTransposeAdd(Na->A,v1,v2,v3);CHKERRQ(ierr);
25 PetscFunctionReturn(0);
26 }
27
MatMultTranspose_Transpose(Mat N,Vec x,Vec y)28 PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y)
29 {
30 Mat_Transpose *Na = (Mat_Transpose*)N->data;
31 PetscErrorCode ierr;
32
33 PetscFunctionBegin;
34 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr);
35 PetscFunctionReturn(0);
36 }
37
MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)38 PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
39 {
40 Mat_Transpose *Na = (Mat_Transpose*)N->data;
41 PetscErrorCode ierr;
42
43 PetscFunctionBegin;
44 ierr = MatMultAdd(Na->A,v1,v2,v3);CHKERRQ(ierr);
45 PetscFunctionReturn(0);
46 }
47
MatDestroy_Transpose(Mat N)48 PetscErrorCode MatDestroy_Transpose(Mat N)
49 {
50 Mat_Transpose *Na = (Mat_Transpose*)N->data;
51 PetscErrorCode ierr;
52
53 PetscFunctionBegin;
54 ierr = MatDestroy(&Na->A);CHKERRQ(ierr);
55 ierr = PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);CHKERRQ(ierr);
56 ierr = PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);CHKERRQ(ierr);
57 ierr = PetscFree(N->data);CHKERRQ(ierr);
58 PetscFunctionReturn(0);
59 }
60
MatDuplicate_Transpose(Mat N,MatDuplicateOption op,Mat * m)61 PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
62 {
63 Mat_Transpose *Na = (Mat_Transpose*)N->data;
64 PetscErrorCode ierr;
65
66 PetscFunctionBegin;
67 if (op == MAT_COPY_VALUES) {
68 ierr = MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr);
69 } else if (op == MAT_DO_NOT_COPY_VALUES) {
70 ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr);
71 ierr = MatTranspose(*m,MAT_INPLACE_MATRIX,m);CHKERRQ(ierr);
72 } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
73 PetscFunctionReturn(0);
74 }
75
MatCreateVecs_Transpose(Mat A,Vec * r,Vec * l)76 PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
77 {
78 Mat_Transpose *Aa = (Mat_Transpose*)A->data;
79 PetscErrorCode ierr;
80
81 PetscFunctionBegin;
82 ierr = MatCreateVecs(Aa->A,l,r);CHKERRQ(ierr);
83 PetscFunctionReturn(0);
84 }
85
MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)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;
92 PetscErrorCode ierr;
93
94 PetscFunctionBegin;
95 ierr = MatAXPY(M,a,N,str);CHKERRQ(ierr);
96 PetscFunctionReturn(0);
97 }
98
MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool * has)99 PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
100 {
101 Mat_Transpose *X = (Mat_Transpose*)mat->data;
102 PetscErrorCode ierr;
103 PetscFunctionBegin;
104
105 *has = PETSC_FALSE;
106 if (op == MATOP_MULT) {
107 ierr = MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);CHKERRQ(ierr);
108 } else if (op == MATOP_MULT_TRANSPOSE) {
109 ierr = MatHasOperation(X->A,MATOP_MULT,has);CHKERRQ(ierr);
110 } else if (op == MATOP_MULT_ADD) {
111 ierr = MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);CHKERRQ(ierr);
112 } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
113 ierr = MatHasOperation(X->A,MATOP_MULT_ADD,has);CHKERRQ(ierr);
114 } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
115 PetscFunctionReturn(0);
116 }
117
118 /* used by hermitian transpose */
MatProductSetFromOptions_Transpose(Mat D)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;
125 PetscErrorCode ierr;
126
127 PetscFunctionBegin;
128 MatCheckProduct(D,1);
129 A = D->product->A;
130 B = D->product->B;
131 C = D->product->C;
132 ierr = PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);CHKERRQ(ierr);
133 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);CHKERRQ(ierr);
134 ierr = PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);CHKERRQ(ierr);
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 ierr = MatTransposeGetMat(Ain,&Ain);CHKERRQ(ierr);
141 ierr = PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);CHKERRQ(ierr);
142 }
143 Btrans = 0;
144 Bin = B;
145 while (Bistrans) {
146 Btrans++;
147 ierr = MatTransposeGetMat(Bin,&Bin);CHKERRQ(ierr);
148 ierr = PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);CHKERRQ(ierr);
149 }
150 Ctrans = 0;
151 Cin = C;
152 while (Cistrans) {
153 Ctrans++;
154 ierr = MatTransposeGetMat(Cin,&Cin);CHKERRQ(ierr);
155 ierr = PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&Cistrans);CHKERRQ(ierr);
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;
164
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 ierr = MatProductReplaceMats(Ain,Bin,Cin,D);CHKERRQ(ierr);
216 ierr = MatProductSetType(D,ptype);CHKERRQ(ierr);
217 ierr = MatProductSetFromOptions(D);CHKERRQ(ierr);
218 PetscFunctionReturn(0);
219 }
220
MatTransposeGetMat_Transpose(Mat A,Mat * M)221 PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
222 {
223 Mat_Transpose *Aa = (Mat_Transpose*)A->data;
224
225 PetscFunctionBegin;
226 *M = Aa->A;
227 PetscFunctionReturn(0);
228 }
229
230 /*@
231 MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT
232
233 Logically collective on Mat
234
235 Input Parameter:
236 . A - the MATTRANSPOSE matrix
237
238 Output Parameter:
239 . M - the matrix object stored inside A
240
241 Level: intermediate
242
243 .seealso: MatCreateTranspose()
244
245 @*/
MatTransposeGetMat(Mat A,Mat * M)246 PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
247 {
248 PetscErrorCode ierr;
249
250 PetscFunctionBegin;
251 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
252 PetscValidType(A,1);
253 PetscValidPointer(M,2);
254 ierr = PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));CHKERRQ(ierr);
255 PetscFunctionReturn(0);
256 }
257
258 /*@
259 MatCreateTranspose - Creates a new matrix object that behaves like A'
260
261 Collective on Mat
262
263 Input Parameter:
264 . A - the (possibly rectangular) matrix
265
266 Output Parameter:
267 . N - the matrix that represents A'
268
269 Level: intermediate
270
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
275
276 .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()
277
278 @*/
MatCreateTranspose(Mat A,Mat * N)279 PetscErrorCode MatCreateTranspose(Mat A,Mat *N)
280 {
281 PetscErrorCode ierr;
282 PetscInt m,n;
283 Mat_Transpose *Na;
284 VecType vtype;
285
286 PetscFunctionBegin;
287 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
288 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
289 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
290 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr);
291 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr);
292 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr);
293
294 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr);
295 (*N)->data = (void*) Na;
296 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
297 Na->A = A;
298
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;
310
311 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);CHKERRQ(ierr);
312 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);CHKERRQ(ierr);
313 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr);
314 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr);
315 ierr = MatSetVecType(*N,vtype);CHKERRQ(ierr);
316 ierr = MatSetUp(*N);CHKERRQ(ierr);
317 PetscFunctionReturn(0);
318 }
319