[Mirror] MueLu: Support BlockedCrsMatrix in Utilities::Transpose#115
[Mirror] MueLu: Support BlockedCrsMatrix in Utilities::Transpose#115
Conversation
…rsMatrix Signed-off-by: malphil <malphil@sandia.gov>
Signed-off-by: malphil <malphil@sandia.gov>
Fix unit test to allow comparison via BlockedCrsMatrix::Merge Signed-off-by: malphil <malphil@sandia.gov>
|
CDash for AT1 results [Only accessible from Sandia networks] |
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request extends the Utilities::Transpose function to support BlockedCrsMatrix by recursively transposing individual blocks. It also introduces a new unit test, TransposeBlockedCrsMatrix, to verify the correctness of this implementation. The review feedback identifies several critical improvements: using full type names instead of short names in the template definition file to ensure compilation, adding null checks for individual blocks to prevent runtime crashes, and correcting the row entry count used when initializing the new blocked matrix. All review comments provide actionable feedback and should be addressed.
| RCP<ParameterList> transposeParams = params.is_null() ? rcp(new ParameterList) : rcp(new ParameterList(*params)); | ||
| transposeParams->set("sort", false); | ||
| A = transposer.createTranspose(transposeParams); | ||
| auto blockOp = rcp_dynamic_cast<BlockedCrsMatrix>(rcpFromRef(Op)); |
There was a problem hiding this comment.
The type BlockedCrsMatrix is a short name that is typically not available in template definition files (_def.hpp) unless explicitly defined or included via MueLu_UseShortNames.hpp. It is recommended to use the full type name Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> to ensure compilation across different environments.
auto blockOp = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Teuchos::rcpFromRef(Op));| auto rangeMaps = blockOp->getBlockedDomainMap(); | ||
| auto domainMaps = blockOp->getBlockedRangeMap(); | ||
|
|
||
| auto blockOpT = make_rcp<BlockedCrsMatrix>(rangeMaps, domainMaps, numEntPerRow); |
There was a problem hiding this comment.
Similar to the cast above, use the full type name for BlockedCrsMatrix. Additionally, using Teuchos::rcp(new ...) is more consistent with the existing code in this file than make_rcp.
auto blockOpT = Teuchos::rcp(new Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(rangeMaps, domainMaps, numEntPerRow));| auto A_ij = blockOp->getMatrix(row, col); | ||
| auto A_ij_T = Utilities::Transpose(*A_ij, false, label, params); | ||
| blockOpT->setMatrix(col, row, A_ij_T); |
There was a problem hiding this comment.
In a BlockedCrsMatrix, individual blocks can be null (representing zero blocks). Dereferencing A_ij without checking for null will cause a crash. The code should verify that the block exists before attempting to transpose it.
auto A_ij = blockOp->getMatrix(row, col);
if (!A_ij.is_null()) {
auto A_ij_T = Utilities::Transpose(*A_ij, false, label, params);
blockOpT->setMatrix(col, row, A_ij_T);
}| A = transposer.createTranspose(transposeParams); | ||
| auto blockOp = rcp_dynamic_cast<BlockedCrsMatrix>(rcpFromRef(Op)); | ||
| if (blockOp != Teuchos::null) { | ||
| auto numEntPerRow = blockOp->getLocalMaxNumRowEntries(); |
There was a problem hiding this comment.
blockOp->getLocalMaxNumRowEntries() returns the maximum number of non-zero scalar entries in any row. However, the BlockedCrsMatrix constructor expects nMaxEntriesPerRow to be the maximum number of blocks per row. Using blockOp->Cols() (the number of block columns) provides a more appropriate upper bound for the block count.
auto numEntPerRow = blockOp->Cols();2 participants
Automated mirror of upstream PR trilinos#15202 @trilinos/muelu