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 (*F)->factortype = ftype; 247 PetscCall(MatDestroy(&FA)); 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 251 static PetscErrorCode MatDestroy_Transpose(Mat N) 252 { 253 Mat A; 254 255 PetscFunctionBegin; 256 PetscCall(MatShellGetContext(N, &A)); 257 PetscCall(MatDestroy(&A)); 258 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL)); 259 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL)); 260 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL)); 261 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatFactorGetSolverType_C", NULL)); 262 PetscFunctionReturn(PETSC_SUCCESS); 263 } 264 265 static PetscErrorCode MatGetInfo_Transpose(Mat N, MatInfoType flag, MatInfo *info) 266 { 267 Mat A; 268 269 PetscFunctionBegin; 270 PetscCall(MatShellGetContext(N, &A)); 271 PetscCall(MatGetInfo(A, flag, info)); 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 static PetscErrorCode MatFactorGetSolverType_Transpose(Mat N, MatSolverType *type) 276 { 277 Mat A; 278 279 PetscFunctionBegin; 280 PetscCall(MatShellGetContext(N, &A)); 281 PetscCall(MatFactorGetSolverType(A, type)); 282 PetscFunctionReturn(PETSC_SUCCESS); 283 } 284 285 static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m) 286 { 287 Mat A, C; 288 289 PetscFunctionBegin; 290 PetscCall(MatShellGetContext(N, &A)); 291 PetscCall(MatDuplicate(A, op, &C)); 292 PetscCall(MatCreateTranspose(C, m)); 293 if (op == MAT_COPY_VALUES) { 294 PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN)); 295 PetscCall(MatPropagateSymmetryOptions(A, C)); 296 } 297 PetscCall(MatDestroy(&C)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300 301 static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has) 302 { 303 Mat A; 304 305 PetscFunctionBegin; 306 PetscCall(MatShellGetContext(mat, &A)); 307 *has = PETSC_FALSE; 308 if (op == MATOP_MULT || op == MATOP_MULT_ADD) { 309 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has)); 310 } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 311 PetscCall(MatHasOperation(A, MATOP_MULT, has)); 312 } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE; 313 PetscFunctionReturn(PETSC_SUCCESS); 314 } 315 316 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D) 317 { 318 Mat A, B, C, Ain, Bin, Cin; 319 PetscBool Aistrans, Bistrans, Cistrans; 320 PetscInt Atrans, Btrans, Ctrans; 321 MatProductType ptype; 322 323 PetscFunctionBegin; 324 MatCheckProduct(D, 1); 325 A = D->product->A; 326 B = D->product->B; 327 C = D->product->C; 328 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans)); 329 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans)); 330 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans)); 331 PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen"); 332 Atrans = 0; 333 Ain = A; 334 while (Aistrans) { 335 Atrans++; 336 PetscCall(MatTransposeGetMat(Ain, &Ain)); 337 PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans)); 338 } 339 Btrans = 0; 340 Bin = B; 341 while (Bistrans) { 342 Btrans++; 343 PetscCall(MatTransposeGetMat(Bin, &Bin)); 344 PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans)); 345 } 346 Ctrans = 0; 347 Cin = C; 348 while (Cistrans) { 349 Ctrans++; 350 PetscCall(MatTransposeGetMat(Cin, &Cin)); 351 PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans)); 352 } 353 Atrans = Atrans % 2; 354 Btrans = Btrans % 2; 355 Ctrans = Ctrans % 2; 356 ptype = D->product->type; /* same product type by default */ 357 if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0; 358 if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0; 359 if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0; 360 361 if (Atrans || Btrans || Ctrans) { 362 ptype = MATPRODUCT_UNSPECIFIED; 363 switch (D->product->type) { 364 case MATPRODUCT_AB: 365 if (Atrans && Btrans) { /* At * Bt we do not have support for this */ 366 /* TODO custom implementation ? */ 367 } else if (Atrans) { /* At * B */ 368 ptype = MATPRODUCT_AtB; 369 } else { /* A * Bt */ 370 ptype = MATPRODUCT_ABt; 371 } 372 break; 373 case MATPRODUCT_AtB: 374 if (Atrans && Btrans) { /* A * Bt */ 375 ptype = MATPRODUCT_ABt; 376 } else if (Atrans) { /* A * B */ 377 ptype = MATPRODUCT_AB; 378 } else { /* At * Bt we do not have support for this */ 379 /* TODO custom implementation ? */ 380 } 381 break; 382 case MATPRODUCT_ABt: 383 if (Atrans && Btrans) { /* At * B */ 384 ptype = MATPRODUCT_AtB; 385 } else if (Atrans) { /* At * Bt we do not have support for this */ 386 /* TODO custom implementation ? */ 387 } else { /* A * B */ 388 ptype = MATPRODUCT_AB; 389 } 390 break; 391 case MATPRODUCT_PtAP: 392 if (Atrans) { /* PtAtP */ 393 /* TODO custom implementation ? */ 394 } else { /* RARt */ 395 ptype = MATPRODUCT_RARt; 396 } 397 break; 398 case MATPRODUCT_RARt: 399 if (Atrans) { /* RAtRt */ 400 /* TODO custom implementation ? */ 401 } else { /* PtAP */ 402 ptype = MATPRODUCT_PtAP; 403 } 404 break; 405 case MATPRODUCT_ABC: 406 /* TODO custom implementation ? */ 407 break; 408 default: 409 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]); 410 } 411 } 412 PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D)); 413 PetscCall(MatProductSetType(D, ptype)); 414 PetscCall(MatProductSetFromOptions(D)); 415 PetscFunctionReturn(PETSC_SUCCESS); 416 } 417 418 static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v) 419 { 420 Mat A; 421 422 PetscFunctionBegin; 423 PetscCall(MatShellGetContext(N, &A)); 424 PetscCall(MatGetDiagonal(A, v)); 425 PetscFunctionReturn(PETSC_SUCCESS); 426 } 427 428 static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str) 429 { 430 Mat a, b; 431 432 PetscFunctionBegin; 433 PetscCall(MatShellGetContext(A, &a)); 434 PetscCall(MatShellGetContext(B, &b)); 435 PetscCall(MatCopy(a, b, str)); 436 PetscFunctionReturn(PETSC_SUCCESS); 437 } 438 439 static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat) 440 { 441 Mat A; 442 PetscScalar vscale = 1.0, vshift = 0.0; 443 PetscBool flg; 444 445 PetscFunctionBegin; 446 PetscCall(MatShellGetContext(N, &A)); 447 PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg)); 448 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 */ 449 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"); 450 PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatAXPY() has been called on the input Mat"); 451 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"); 452 PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalSet() has been called on the input Mat"); 453 vscale = ((Mat_Shell *)N->data)->vscale; 454 vshift = ((Mat_Shell *)N->data)->vshift; 455 } 456 if (flg) { 457 Mat B; 458 459 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B)); 460 if (reuse != MAT_INPLACE_MATRIX) { 461 PetscCall(MatConvert(B, newtype, reuse, newmat)); 462 PetscCall(MatDestroy(&B)); 463 } else { 464 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B)); 465 PetscCall(MatHeaderReplace(N, &B)); 466 } 467 } else { /* use basic converter as fallback */ 468 flg = (PetscBool)(N->ops->getrow != NULL); 469 PetscCall(MatConvert_Basic(N, newtype, reuse, newmat)); 470 } 471 if (flg) { 472 PetscCall(MatScale(*newmat, vscale)); 473 PetscCall(MatShift(*newmat, vshift)); 474 } 475 PetscFunctionReturn(PETSC_SUCCESS); 476 } 477 478 static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M) 479 { 480 PetscFunctionBegin; 481 PetscCall(MatShellGetContext(N, M)); 482 PetscFunctionReturn(PETSC_SUCCESS); 483 } 484 485 /*@ 486 MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL` 487 488 Logically Collective 489 490 Input Parameter: 491 . A - the `MATTRANSPOSEVIRTUAL` matrix 492 493 Output Parameter: 494 . M - the matrix object stored inside `A` 495 496 Level: intermediate 497 498 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()` 499 @*/ 500 PetscErrorCode MatTransposeGetMat(Mat A, Mat *M) 501 { 502 PetscFunctionBegin; 503 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 504 PetscValidType(A, 1); 505 PetscAssertPointer(M, 2); 506 PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M)); 507 PetscFunctionReturn(PETSC_SUCCESS); 508 } 509 510 /*MC 511 MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix 512 513 Level: advanced 514 515 Developer Notes: 516 This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code 517 518 Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage 519 520 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`, 521 `MATNORMALHERMITIAN`, `MATNORMAL` 522 M*/ 523 524 /*@ 525 MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A' 526 527 Collective 528 529 Input Parameter: 530 . A - the (possibly rectangular) matrix 531 532 Output Parameter: 533 . N - the matrix that represents A' 534 535 Level: intermediate 536 537 Note: 538 The transpose A' is NOT actually formed! Rather the new matrix 539 object performs the matrix-vector product by using the `MatMultTranspose()` on 540 the original matrix 541 542 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`, 543 `MATNORMALHERMITIAN` 544 @*/ 545 PetscErrorCode MatCreateTranspose(Mat A, Mat *N) 546 { 547 VecType vtype; 548 549 PetscFunctionBegin; 550 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 551 PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap))); 552 PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap))); 553 PetscCall(MatSetType(*N, MATSHELL)); 554 PetscCall(MatShellSetContext(*N, A)); 555 PetscCall(PetscObjectReference((PetscObject)A)); 556 557 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 558 PetscCall(MatGetVecType(A, &vtype)); 559 PetscCall(MatSetVecType(*N, vtype)); 560 #if defined(PETSC_HAVE_DEVICE) 561 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 562 #endif 563 PetscCall(MatSetUp(*N)); 564 565 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Transpose)); 566 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Transpose)); 567 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Transpose)); 568 PetscCall(MatShellSetOperation(*N, MATOP_LUFACTOR, (void (*)(void))MatLUFactor_Transpose)); 569 PetscCall(MatShellSetOperation(*N, MATOP_CHOLESKYFACTOR, (void (*)(void))MatCholeskyFactor_Transpose)); 570 PetscCall(MatShellSetOperation(*N, MATOP_GET_FACTOR, (void (*)(void))MatGetFactor_Transpose)); 571 PetscCall(MatShellSetOperation(*N, MATOP_GETINFO, (void (*)(void))MatGetInfo_Transpose)); 572 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Transpose)); 573 PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_Transpose)); 574 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Transpose)); 575 PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Transpose)); 576 PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_Transpose)); 577 578 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose)); 579 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose)); 580 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatFactorGetSolverType_C", MatFactorGetSolverType_Transpose)); 581 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable)); 582 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable)); 583 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable)); 584 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL)); 585 PetscFunctionReturn(PETSC_SUCCESS); 586 } 587