Add ldiv! for LU decomposition (#1532)

Author: SBuercklin

lib/cusolver/linalg.jl

@@ -338,6 +338,16 @@ LinearAlgebra.ipiv2perm(v::CuVector{T}, maxi::Integer) where T =
 
 end
 
+function LinearAlgebra.ldiv!(F::LU{T,<:StridedCuMatrix{T}}, B::CuVecOrMat{T}) where {T}
+    return getrs!('N', F.factors, F.ipiv, B)
+end
+
+# LinearAlgebra.jl defines a comparable method with these restrictions on T, so we need
+#   to define with the same type parameters to resolve method ambiguity for these cases
+function LinearAlgebra.ldiv!(F::LU{T,<:StridedCuMatrix{T}}, B::CuVecOrMat{T}) where {T <: Union{Float32, Float64, ComplexF32, ComplexF64}}
+    return getrs!('N', F.factors, F.ipiv, B)
+end
+
 ## cholesky
 
 if VERSION >= v"1.8-"

test/cusolver/dense.jl

@@ -579,6 +579,15 @@ k = 1
 
             @test_throws LinearAlgebra.SingularException lu(CUDA.zeros(elty,n,n))
         end
+        @testset "lu ldiv! elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
+            A = rand(elty, m, m)
+            B = rand(elty, m, m)
+            A_d = CuArray(A)
+            B_d = CuArray(B)
+            lu_cpu = lu(A)
+            lu_gpu = lu(A_d)
+            @test ldiv!(lu_cpu, B) ≈ collect(ldiv!(lu_gpu, B_d))
+        end
     end
 end