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