Ifpack2: Plug in HTS for thread-parallel sparse triangular solve

[mhoemmen]

@trilinos/ifpack2 @trilinos/tpetra

[mhoemmen]

This depended on #514. Ifpack2::LocalSparseTriangularSolver is the right place to put this, NOT Tpetra::CrsMatrix::localSolve. See #581.

[egboman]

I suppose the triangular solve codes might live in KokkosKernels or ShyLu, but used by Ifpack2?

[mhoemmen]

@egboman Main interface to local sparse triangular solves is in KokkosKernels. I'm wrapping MKL's sparse TRSM right now. Ifpack2::LocalSparseTriangularSolver just wraps the KokkosKernels interface.

[ambrad]

@mhoemmen, in RILUK::apply(), you wrote (Jan 2014): "C must have the same Map as D. We have to use a temp multivector, since our implementation of triangular solves does not allow its input and output to alias one another." However, I believe that is no longer true. CrsMatrix::localSolve's doc says: "...unlike multiply(), it is permissible for &X == &Y." Also, the tests pass without the temporary. Anything I'm missing? Thanks.

[mhoemmen]

@ambrad wrote:

... CrsMatrix::localSolve's doc says: "...unlike multiply(), it is permissible for &X == &Y."

Would that be true for a general parallel triangular solve, such as what HTS does?

I haven't had time to audit all the comments for correctness, alas ;-)

[ambrad]

@mhoemmen, my view is that if a particular trisolver can't handle it, that should be handled in the impl details of LocalSparseTriangularSolver, not in RILUK. That's what I was about to start doing, but then I thought to check whether it's still true for CrsMatrix::localSolve. I concluded that, at present, neither option (CrsMatrix::localSolve or HTS) requires &X != &Y.

[ambrad]

BTW, work for this issue is in this branch of my fork: https://github.com/ambrad/Trilinos/commits/hts1. (I'll squash it to one commit when I'm done.)

[mhoemmen]

@ambrad I don't particularly like the promise that aliasing works. &X == &Y does not suffice to check all possible aliasing cases (e.g., Tpetra::MultiVector that views a subset of columns of another Tpetra::MultiVector). Also, even if aliasing works for the local kernel, we might like to overlap communication and computation at the MPI level. That could affect whether aliasing is correct.

[ambrad]

@mhoemmen, that's fine, but do you agree that LocalSparseTriangularSolver is the right place to manage aliasing? I.e., RILUK should not have to know whether a temp MV needs to be made; it should be able to call the Trisolver and let it figure it out.

[mhoemmen]

@ambrad wrote:

... do you agree that LocalSparseTriangularSolver is the right place to manage aliasing?

Yes, I agree. I don't think we should ask KokkosKernels or third-party implementations to do that. The Ifpack2 interface should handle it. (It already does, for some other solvers.)

[ambrad]

@trilinos/ifpack2, @trilinos/shylu, I'm proposing to push the commit linked above. It is the first of I anticipate >=3 commits focusing on Ifpack2::LocalSparseTriangularSolver. @mhoemmen and I might still discuss some interface issues re: aliasing in LSTS, but I propose that I push this commit, anyway, since it gets a fair bit done, and include further decisions on aliasing in a future commit. Thoughts? Objections to my pushing this commit?

[mhoemmen]

@ambrad I think it's OK to push for now. The implementation eventually should tie the sparse triangular solve implementation to MatrixType::execution _space (e.g., cuSPARSE for Cuda, HTS for OpenMP, "Internal" for Serial). Thanks!

[srajama1]

@ambrad : Can we close this ? The issue is quite broad, but I assume the primary target is integrating HTS . Integrating cuSPARSE or an internal one is planned later this year and can be addressed in a separate issue.

[ambrad]

I agree, @srajama1. I will close this now, with GPU deliverables to go into a new issue.

[ambrad]

Sorry, we should ask @mhoemmen first. @mhoemmen, close this issue if you agree with putting GPU trisolver work into a new issue.

[mhoemmen]

@ambrad I agree; thanks! The new issue (for GPU sparse triangular solve) is #1790.