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