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 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 ierr = MatMultTranspose(Na->A,x,y);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 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 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 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 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 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 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 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 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 */ 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 221 PetscErrorCode MatGetDiagonal_Transpose(Mat A,Vec v) 222 { 223 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 224 PetscErrorCode ierr; 225 226 PetscFunctionBegin; 227 ierr = MatGetDiagonal(Aa->A,v);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 PetscErrorCode MatConvert_Transpose(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 232 { 233 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 234 PetscErrorCode ierr; 235 PetscBool flg; 236 237 PetscFunctionBegin; 238 ierr = MatHasOperation(Aa->A,MATOP_TRANSPOSE,&flg);CHKERRQ(ierr); 239 if (flg) { 240 Mat B; 241 242 ierr = MatTranspose(Aa->A,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 243 if (reuse != MAT_INPLACE_MATRIX) { 244 ierr = MatConvert(B,newtype,reuse,newmat);CHKERRQ(ierr); 245 ierr = MatDestroy(&B);CHKERRQ(ierr); 246 } else { 247 ierr = MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr); 248 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 249 } 250 } else { /* use basic converter as fallback */ 251 ierr = MatConvert_Basic(A,newtype,reuse,newmat);CHKERRQ(ierr); 252 } 253 PetscFunctionReturn(0); 254 } 255 256 PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M) 257 { 258 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 259 260 PetscFunctionBegin; 261 *M = Aa->A; 262 PetscFunctionReturn(0); 263 } 264 265 /*@ 266 MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT 267 268 Logically collective on Mat 269 270 Input Parameter: 271 . A - the MATTRANSPOSE matrix 272 273 Output Parameter: 274 . M - the matrix object stored inside A 275 276 Level: intermediate 277 278 .seealso: MatCreateTranspose() 279 280 @*/ 281 PetscErrorCode MatTransposeGetMat(Mat A,Mat *M) 282 { 283 PetscErrorCode ierr; 284 285 PetscFunctionBegin; 286 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 287 PetscValidType(A,1); 288 PetscValidPointer(M,2); 289 ierr = PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 /*@ 294 MatCreateTranspose - Creates a new matrix object that behaves like A' 295 296 Collective on Mat 297 298 Input Parameter: 299 . A - the (possibly rectangular) matrix 300 301 Output Parameter: 302 . N - the matrix that represents A' 303 304 Level: intermediate 305 306 Notes: 307 The transpose A' is NOT actually formed! Rather the new matrix 308 object performs the matrix-vector product by using the MatMultTranspose() on 309 the original matrix 310 311 .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate() 312 313 @*/ 314 PetscErrorCode MatCreateTranspose(Mat A,Mat *N) 315 { 316 PetscErrorCode ierr; 317 PetscInt m,n; 318 Mat_Transpose *Na; 319 VecType vtype; 320 321 PetscFunctionBegin; 322 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 323 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 324 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 325 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr); 326 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr); 327 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr); 328 329 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 330 (*N)->data = (void*) Na; 331 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 332 Na->A = A; 333 334 (*N)->ops->destroy = MatDestroy_Transpose; 335 (*N)->ops->mult = MatMult_Transpose; 336 (*N)->ops->multadd = MatMultAdd_Transpose; 337 (*N)->ops->multtranspose = MatMultTranspose_Transpose; 338 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose; 339 (*N)->ops->duplicate = MatDuplicate_Transpose; 340 (*N)->ops->getvecs = MatCreateVecs_Transpose; 341 (*N)->ops->axpy = MatAXPY_Transpose; 342 (*N)->ops->hasoperation = MatHasOperation_Transpose; 343 (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose; 344 (*N)->ops->getdiagonal = MatGetDiagonal_Transpose; 345 (*N)->ops->convert = MatConvert_Transpose; 346 (*N)->assembled = PETSC_TRUE; 347 348 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);CHKERRQ(ierr); 349 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);CHKERRQ(ierr); 350 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 351 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr); 352 ierr = MatSetVecType(*N,vtype);CHKERRQ(ierr); 353 #if defined(PETSC_HAVE_DEVICE) 354 ierr = MatBindToCPU(*N,A->boundtocpu);CHKERRQ(ierr); 355 #endif 356 ierr = MatSetUp(*N);CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359