1 #include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/ 2 3 static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y) 4 { 5 Mat A; 6 7 PetscFunctionBegin; 8 PetscCall(MatShellGetContext(N, &A)); 9 PetscCall(MatMultTranspose(A, x, y)); 10 PetscFunctionReturn(PETSC_SUCCESS); 11 } 12 13 static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y) 14 { 15 Mat A; 16 17 PetscFunctionBegin; 18 PetscCall(MatShellGetContext(N, &A)); 19 PetscCall(MatMult(A, x, y)); 20 PetscFunctionReturn(PETSC_SUCCESS); 21 } 22 23 static PetscErrorCode MatSolve_Transpose_LU(Mat N, Vec b, Vec x) 24 { 25 Mat A; 26 27 PetscFunctionBegin; 28 PetscCall(MatShellGetContext(N, &A)); 29 PetscCall(MatSolveTranspose(A, b, x)); 30 PetscFunctionReturn(PETSC_SUCCESS); 31 } 32 33 static PetscErrorCode MatSolveAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x) 34 { 35 Mat A; 36 37 PetscFunctionBegin; 38 PetscCall(MatShellGetContext(N, &A)); 39 PetscCall(MatSolveTransposeAdd(A, b, y, x)); 40 PetscFunctionReturn(PETSC_SUCCESS); 41 } 42 43 static PetscErrorCode MatSolveTranspose_Transpose_LU(Mat N, Vec b, Vec x) 44 { 45 Mat A; 46 47 PetscFunctionBegin; 48 PetscCall(MatShellGetContext(N, &A)); 49 PetscCall(MatSolve(A, b, x)); 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode MatSolveTransposeAdd_Transpose_LU(Mat N, Vec b, Vec y, Vec x) 54 { 55 Mat A; 56 57 PetscFunctionBegin; 58 PetscCall(MatShellGetContext(N, &A)); 59 PetscCall(MatSolveAdd(A, b, y, x)); 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 static PetscErrorCode MatMatSolve_Transpose_LU(Mat N, Mat B, Mat X) 64 { 65 Mat A; 66 67 PetscFunctionBegin; 68 PetscCall(MatShellGetContext(N, &A)); 69 PetscCall(MatMatSolveTranspose(A, B, X)); 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 73 static PetscErrorCode MatMatSolveTranspose_Transpose_LU(Mat N, Mat B, Mat X) 74 { 75 Mat A; 76 77 PetscFunctionBegin; 78 PetscCall(MatShellGetContext(N, &A)); 79 PetscCall(MatMatSolve(A, B, X)); 80 PetscFunctionReturn(PETSC_SUCCESS); 81 } 82 83 static PetscErrorCode MatLUFactor_Transpose(Mat N, IS row, IS col, const MatFactorInfo *minfo) 84 { 85 Mat A; 86 87 PetscFunctionBegin; 88 PetscCall(MatShellGetContext(N, &A)); 89 PetscCall(MatLUFactor(A, col, row, minfo)); 90 PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (void (*)(void))MatSolve_Transpose_LU)); 91 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_LU)); 92 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_LU)); 93 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_LU)); 94 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_LU)); 95 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_LU)); 96 PetscFunctionReturn(PETSC_SUCCESS); 97 } 98 99 static PetscErrorCode MatSolve_Transpose_Cholesky(Mat N, Vec b, Vec x) 100 { 101 Mat A; 102 103 PetscFunctionBegin; 104 PetscCall(MatShellGetContext(N, &A)); 105 PetscCall(MatSolveTranspose(A, b, x)); 106 PetscFunctionReturn(PETSC_SUCCESS); 107 } 108 109 static PetscErrorCode MatSolveAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x) 110 { 111 Mat A; 112 113 PetscFunctionBegin; 114 PetscCall(MatShellGetContext(N, &A)); 115 PetscCall(MatSolveTransposeAdd(A, b, y, x)); 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 119 static PetscErrorCode MatSolveTranspose_Transpose_Cholesky(Mat N, Vec b, Vec x) 120 { 121 Mat A; 122 123 PetscFunctionBegin; 124 PetscCall(MatShellGetContext(N, &A)); 125 PetscCall(MatSolve(A, b, x)); 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 static PetscErrorCode MatSolveTransposeAdd_Transpose_Cholesky(Mat N, Vec b, Vec y, Vec x) 130 { 131 Mat A; 132 133 PetscFunctionBegin; 134 PetscCall(MatShellGetContext(N, &A)); 135 PetscCall(MatSolveAdd(A, b, y, x)); 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 static PetscErrorCode MatMatSolve_Transpose_Cholesky(Mat N, Mat B, Mat X) 140 { 141 Mat A; 142 143 PetscFunctionBegin; 144 PetscCall(MatShellGetContext(N, &A)); 145 PetscCall(MatMatSolveTranspose(A, B, X)); 146 PetscFunctionReturn(PETSC_SUCCESS); 147 } 148 149 static PetscErrorCode MatMatSolveTranspose_Transpose_Cholesky(Mat N, Mat B, Mat X) 150 { 151 Mat A; 152 153 PetscFunctionBegin; 154 PetscCall(MatShellGetContext(N, &A)); 155 PetscCall(MatMatSolve(A, B, X)); 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158 159 static PetscErrorCode MatCholeskyFactor_Transpose(Mat N, IS perm, const MatFactorInfo *minfo) 160 { 161 Mat A; 162 163 PetscFunctionBegin; 164 PetscCall(MatShellGetContext(N, &A)); 165 PetscCall(MatCholeskyFactor(A, perm, minfo)); 166 PetscCall(MatShellSetOperation(N, MATOP_SOLVE, (void (*)(void))MatSolve_Transpose_Cholesky)); 167 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_Cholesky)); 168 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_Cholesky)); 169 PetscCall(MatShellSetOperation(N, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_Cholesky)); 170 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_Cholesky)); 171 PetscCall(MatShellSetOperation(N, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_Cholesky)); 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174 175 static PetscErrorCode MatLUFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info) 176 { 177 Mat A, FA; 178 179 PetscFunctionBegin; 180 PetscCall(MatShellGetContext(N, &A)); 181 PetscCall(MatShellGetContext(F, &FA)); 182 PetscCall(MatLUFactorNumeric(FA, A, info)); 183 PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (void (*)(void))MatSolve_Transpose_LU)); 184 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_LU)); 185 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_LU)); 186 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_LU)); 187 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_LU)); 188 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_LU)); 189 PetscFunctionReturn(PETSC_SUCCESS); 190 } 191 192 static PetscErrorCode MatLUFactorSymbolic_Transpose(Mat F, Mat N, IS row, IS col, const MatFactorInfo *info) 193 { 194 Mat A, FA; 195 196 PetscFunctionBegin; 197 PetscCall(MatShellGetContext(N, &A)); 198 PetscCall(MatShellGetContext(F, &FA)); 199 PetscCall(MatLUFactorSymbolic(FA, A, row, col, info)); 200 PetscCall(MatShellSetOperation(F, MATOP_LUFACTOR_NUMERIC, (void (*)(void))MatLUFactorNumeric_Transpose)); 201 PetscFunctionReturn(PETSC_SUCCESS); 202 } 203 204 static PetscErrorCode MatCholeskyFactorNumeric_Transpose(Mat F, Mat N, const MatFactorInfo *info) 205 { 206 Mat A, FA; 207 208 PetscFunctionBegin; 209 PetscCall(MatShellGetContext(N, &A)); 210 PetscCall(MatShellGetContext(F, &FA)); 211 PetscCall(MatCholeskyFactorNumeric(FA, A, info)); 212 PetscCall(MatShellSetOperation(F, MATOP_SOLVE, (void (*)(void))MatSolve_Transpose_Cholesky)); 213 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_ADD, (void (*)(void))MatSolveAdd_Transpose_Cholesky)); 214 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE, (void (*)(void))MatSolveTranspose_Transpose_Cholesky)); 215 PetscCall(MatShellSetOperation(F, MATOP_SOLVE_TRANSPOSE_ADD, (void (*)(void))MatSolveTransposeAdd_Transpose_Cholesky)); 216 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE, (void (*)(void))MatMatSolve_Transpose_Cholesky)); 217 PetscCall(MatShellSetOperation(F, MATOP_MAT_SOLVE_TRANSPOSE, (void (*)(void))MatMatSolveTranspose_Transpose_Cholesky)); 218 PetscFunctionReturn(PETSC_SUCCESS); 219 } 220 221 static PetscErrorCode MatCholeskyFactorSymbolic_Transpose(Mat F, Mat N, IS perm, const MatFactorInfo *info) 222 { 223 Mat A, FA; 224 225 PetscFunctionBegin; 226 PetscCall(MatShellGetContext(N, &A)); 227 PetscCall(MatShellGetContext(F, &FA)); 228 PetscCall(MatCholeskyFactorSymbolic(FA, A, perm, info)); 229 PetscCall(MatShellSetOperation(F, MATOP_CHOLESKY_FACTOR_NUMERIC, (void (*)(void))MatCholeskyFactorNumeric_Transpose)); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 static PetscErrorCode MatGetFactor_Transpose(Mat N, MatSolverType type, MatFactorType ftype, Mat *F) 234 { 235 Mat A, FA; 236 237 PetscFunctionBegin; 238 PetscCall(MatShellGetContext(N, &A)); 239 PetscCall(MatGetFactor(A, type, ftype, &FA)); 240 PetscCall(MatCreateTranspose(FA, F)); 241 if (ftype == MAT_FACTOR_LU) PetscCall(MatShellSetOperation(*F, MATOP_LUFACTOR_SYMBOLIC, (void (*)(void))MatLUFactorSymbolic_Transpose)); 242 else if (ftype == MAT_FACTOR_CHOLESKY) { 243 PetscCall(MatShellSetOperation(*F, MATOP_CHOLESKY_FACTOR_SYMBOLIC, (void (*)(void))MatCholeskyFactorSymbolic_Transpose)); 244 PetscCall(MatPropagateSymmetryOptions(A, FA)); 245 } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Support for factor type %s not implemented in MATTRANSPOSEVIRTUAL", MatFactorTypes[ftype]); 246 PetscCall(MatDestroy(&FA)); 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 static PetscErrorCode MatDestroy_Transpose(Mat N) 251 { 252 Mat A; 253 254 PetscFunctionBegin; 255 PetscCall(MatShellGetContext(N, &A)); 256 PetscCall(MatDestroy(&A)); 257 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL)); 258 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL)); 259 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL)); 260 PetscFunctionReturn(PETSC_SUCCESS); 261 } 262 263 static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m) 264 { 265 Mat A, C; 266 267 PetscFunctionBegin; 268 PetscCall(MatShellGetContext(N, &A)); 269 PetscCall(MatDuplicate(A, op, &C)); 270 PetscCall(MatCreateTranspose(C, m)); 271 if (op == MAT_COPY_VALUES) { 272 PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN)); 273 PetscCall(MatPropagateSymmetryOptions(A, C)); 274 } 275 PetscCall(MatDestroy(&C)); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has) 280 { 281 Mat A; 282 283 PetscFunctionBegin; 284 PetscCall(MatShellGetContext(mat, &A)); 285 *has = PETSC_FALSE; 286 if (op == MATOP_MULT || op == MATOP_MULT_ADD) { 287 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has)); 288 } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 289 PetscCall(MatHasOperation(A, MATOP_MULT, has)); 290 } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE; 291 PetscFunctionReturn(PETSC_SUCCESS); 292 } 293 294 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D) 295 { 296 Mat A, B, C, Ain, Bin, Cin; 297 PetscBool Aistrans, Bistrans, Cistrans; 298 PetscInt Atrans, Btrans, Ctrans; 299 MatProductType ptype; 300 301 PetscFunctionBegin; 302 MatCheckProduct(D, 1); 303 A = D->product->A; 304 B = D->product->B; 305 C = D->product->C; 306 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans)); 307 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans)); 308 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans)); 309 PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen"); 310 Atrans = 0; 311 Ain = A; 312 while (Aistrans) { 313 Atrans++; 314 PetscCall(MatTransposeGetMat(Ain, &Ain)); 315 PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans)); 316 } 317 Btrans = 0; 318 Bin = B; 319 while (Bistrans) { 320 Btrans++; 321 PetscCall(MatTransposeGetMat(Bin, &Bin)); 322 PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans)); 323 } 324 Ctrans = 0; 325 Cin = C; 326 while (Cistrans) { 327 Ctrans++; 328 PetscCall(MatTransposeGetMat(Cin, &Cin)); 329 PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans)); 330 } 331 Atrans = Atrans % 2; 332 Btrans = Btrans % 2; 333 Ctrans = Ctrans % 2; 334 ptype = D->product->type; /* same product type by default */ 335 if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0; 336 if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0; 337 if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0; 338 339 if (Atrans || Btrans || Ctrans) { 340 ptype = MATPRODUCT_UNSPECIFIED; 341 switch (D->product->type) { 342 case MATPRODUCT_AB: 343 if (Atrans && Btrans) { /* At * Bt we do not have support for this */ 344 /* TODO custom implementation ? */ 345 } else if (Atrans) { /* At * B */ 346 ptype = MATPRODUCT_AtB; 347 } else { /* A * Bt */ 348 ptype = MATPRODUCT_ABt; 349 } 350 break; 351 case MATPRODUCT_AtB: 352 if (Atrans && Btrans) { /* A * Bt */ 353 ptype = MATPRODUCT_ABt; 354 } else if (Atrans) { /* A * B */ 355 ptype = MATPRODUCT_AB; 356 } else { /* At * Bt we do not have support for this */ 357 /* TODO custom implementation ? */ 358 } 359 break; 360 case MATPRODUCT_ABt: 361 if (Atrans && Btrans) { /* At * B */ 362 ptype = MATPRODUCT_AtB; 363 } else if (Atrans) { /* At * Bt we do not have support for this */ 364 /* TODO custom implementation ? */ 365 } else { /* A * B */ 366 ptype = MATPRODUCT_AB; 367 } 368 break; 369 case MATPRODUCT_PtAP: 370 if (Atrans) { /* PtAtP */ 371 /* TODO custom implementation ? */ 372 } else { /* RARt */ 373 ptype = MATPRODUCT_RARt; 374 } 375 break; 376 case MATPRODUCT_RARt: 377 if (Atrans) { /* RAtRt */ 378 /* TODO custom implementation ? */ 379 } else { /* PtAP */ 380 ptype = MATPRODUCT_PtAP; 381 } 382 break; 383 case MATPRODUCT_ABC: 384 /* TODO custom implementation ? */ 385 break; 386 default: 387 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]); 388 } 389 } 390 PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D)); 391 PetscCall(MatProductSetType(D, ptype)); 392 PetscCall(MatProductSetFromOptions(D)); 393 PetscFunctionReturn(PETSC_SUCCESS); 394 } 395 396 static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v) 397 { 398 Mat A; 399 400 PetscFunctionBegin; 401 PetscCall(MatShellGetContext(N, &A)); 402 PetscCall(MatGetDiagonal(A, v)); 403 PetscFunctionReturn(PETSC_SUCCESS); 404 } 405 406 static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str) 407 { 408 Mat a, b; 409 410 PetscFunctionBegin; 411 PetscCall(MatShellGetContext(A, &a)); 412 PetscCall(MatShellGetContext(B, &b)); 413 PetscCall(MatCopy(a, b, str)); 414 PetscFunctionReturn(PETSC_SUCCESS); 415 } 416 417 static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat) 418 { 419 Mat A; 420 PetscScalar vscale = 1.0, vshift = 0.0; 421 PetscBool flg; 422 423 PetscFunctionBegin; 424 PetscCall(MatShellGetContext(N, &A)); 425 PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg)); 426 if (flg || N->ops->getrow) { /* if this condition is false, MatConvert_Shell() will be called in MatConvert_Basic(), so the following checks are not needed */ 427 PetscCheck(!((Mat_Shell *)N->data)->zrows && !((Mat_Shell *)N->data)->zcols, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); 428 PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatAXPY() has been called on the input Mat"); 429 PetscCheck(!((Mat_Shell *)N->data)->left && !((Mat_Shell *)N->data)->right, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalScale() has been called on the input Mat"); 430 PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalSet() has been called on the input Mat"); 431 vscale = ((Mat_Shell *)N->data)->vscale; 432 vshift = ((Mat_Shell *)N->data)->vshift; 433 } 434 if (flg) { 435 Mat B; 436 437 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B)); 438 if (reuse != MAT_INPLACE_MATRIX) { 439 PetscCall(MatConvert(B, newtype, reuse, newmat)); 440 PetscCall(MatDestroy(&B)); 441 } else { 442 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B)); 443 PetscCall(MatHeaderReplace(N, &B)); 444 } 445 } else { /* use basic converter as fallback */ 446 flg = (PetscBool)(N->ops->getrow != NULL); 447 PetscCall(MatConvert_Basic(N, newtype, reuse, newmat)); 448 } 449 if (flg) { 450 PetscCall(MatScale(*newmat, vscale)); 451 PetscCall(MatShift(*newmat, vshift)); 452 } 453 PetscFunctionReturn(PETSC_SUCCESS); 454 } 455 456 static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M) 457 { 458 PetscFunctionBegin; 459 PetscCall(MatShellGetContext(N, M)); 460 PetscFunctionReturn(PETSC_SUCCESS); 461 } 462 463 /*@ 464 MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL` 465 466 Logically Collective 467 468 Input Parameter: 469 . A - the `MATTRANSPOSEVIRTUAL` matrix 470 471 Output Parameter: 472 . M - the matrix object stored inside `A` 473 474 Level: intermediate 475 476 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()` 477 @*/ 478 PetscErrorCode MatTransposeGetMat(Mat A, Mat *M) 479 { 480 PetscFunctionBegin; 481 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 482 PetscValidType(A, 1); 483 PetscAssertPointer(M, 2); 484 PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M)); 485 PetscFunctionReturn(PETSC_SUCCESS); 486 } 487 488 /*MC 489 MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix 490 491 Level: advanced 492 493 Developer Notes: 494 This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code 495 496 Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage 497 498 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`, 499 `MATNORMALHERMITIAN`, `MATNORMAL` 500 M*/ 501 502 /*@ 503 MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A' 504 505 Collective 506 507 Input Parameter: 508 . A - the (possibly rectangular) matrix 509 510 Output Parameter: 511 . N - the matrix that represents A' 512 513 Level: intermediate 514 515 Note: 516 The transpose A' is NOT actually formed! Rather the new matrix 517 object performs the matrix-vector product by using the `MatMultTranspose()` on 518 the original matrix 519 520 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`, 521 `MATNORMALHERMITIAN` 522 @*/ 523 PetscErrorCode MatCreateTranspose(Mat A, Mat *N) 524 { 525 VecType vtype; 526 527 PetscFunctionBegin; 528 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 529 PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap))); 530 PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap))); 531 PetscCall(MatSetType(*N, MATSHELL)); 532 PetscCall(MatShellSetContext(*N, A)); 533 PetscCall(PetscObjectReference((PetscObject)A)); 534 535 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 536 PetscCall(MatGetVecType(A, &vtype)); 537 PetscCall(MatSetVecType(*N, vtype)); 538 #if defined(PETSC_HAVE_DEVICE) 539 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 540 #endif 541 PetscCall(MatSetUp(*N)); 542 543 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Transpose)); 544 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Transpose)); 545 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Transpose)); 546 PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (void (*)(void))MatLUFactor_Transpose)); 547 PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (void (*)(void))MatCholeskyFactor_Transpose)); 548 PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (void (*)(void))MatGetFactor_Transpose)); 549 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Transpose)); 550 PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_Transpose)); 551 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Transpose)); 552 PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Transpose)); 553 PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_Transpose)); 554 555 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose)); 556 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose)); 557 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable)); 558 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable)); 559 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable)); 560 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL)); 561 PetscFunctionReturn(PETSC_SUCCESS); 562 } 563