[CUSPARSE] Update the interface for triangular solves

[amontoison]

@dkarrasch
It seems that NVIDIA reintroduced in-place triangular solves in CUSPARSE with recent CUDA toolkits.
I also found a few typos that you added in #1946. Can you review the PR?

Can you also confirm that ldiv!(C, A, B), ldiv!(A, B) and A \ B and automatically defined for the four triangular wrappers if we define generic_trimatdiv!?

[amontoison]

@dkarrasch
Is possible to modify the function

generic_trimatdiv!(C, uploc, isunitc, tfun, A, B)

into

generic_trimatdiv!(C, uploc, isunitc, tfunA, A, tfunB, B)

?

If the right-hand side B is transposed, the current version of generic_trimatdiv! is not working on Julia 1.10.
Julia doesn't dispatch to it on some cases and we have scalar indexing.

[dkarrasch]

No, at least that it is not designed like that. The reason is that triangular solves are typically 2-arg only. But you could create your own.

generic_trimatdiv!(C, uploc, isunitc, tfun, A, B::SomeCuMatrix) = # as usual
generic_trimatdiv!(C, uploc, isunitc, tfun, A, B::Transpose{<:Any,<:SomeCuMatrix}) =
    _generic_trimatdiv!(C, uploc, isunitc, tfunA, A, tfunB, B)

(with the other slots type-annotated as appropriate), which then calls whatever CUSPARSE offers.

[amontoison]

I already did that but in the case that the triangular solves are in-place, C is also in a Transpose / Adjoint and ldiv! doesn't dispatch to generic_trimatdiv!.

We can use 2-arg methods again with CUDA v12.x so it could simplify many things.
I changed everything last year to use the 3-args ldiv! because the 2-args ldiv! was not anymore supported but it's fixed now.

[dkarrasch]

The dispatch in LinearAlgebra is as follows: if you call ldiv!(T, B), then this gets expanded to the "three-arg" generic_trimatdiv! with the character for potential transposition. This is where packages like this one dispatch one the underlying storage. All the BLAS methods then check if C === B, and in that case ignore B completely and solve in-place of C. So it is a 2-arg solve eventually. For (Open)BLAS, I don't think you can call ldiv! on a transpose rhs directly (I mean you can, but you will end up in a generic solve routine). You can only call T\ transpose(B), and then transpose(B) gets materialized and then mutated. The same happens here: sm2! and sv2! are 2-arg solves.

[amontoison]

Thanks @dkarrasch! I understand how the dispatch with generic_trimatdiv! works now.
I solved the issue by adding the following three three methods:

function LinearAlgebra.generic_trimatdiv!(C::DenseCuMatrix{T}, uploc, isunitc, tfun, A, B::AdjOrTrans{T,<:DenseCuMatrix{T}}) where {T<:BlasFloat}
  ...
end

function LinearAlgebra.generic_trimatdiv!(C::Transpose{T,<:DenseCuMatrix{T}}, uploc, isunitc, tfun, A, B::Transpose{T,<:DenseCuMatrix{T}}) where {T<:BlasFloat}
  ...
end

function LinearAlgebra.generic_trimatdiv!(C::Adjoint{T,<:DenseCuMatrix{T}}, uploc, isunitc, tfun, A, B::Adjoint{T,<:DenseCuMatrix{T}}) where {T<:BlasFloat}
  ...
end