/-Operator for LU (#583)

* Added rdiv for LU type. Added a test too.

Author: judober

src/lu.jl

@@ -141,3 +141,5 @@ end
 
 \(F::LU, v::AbstractVector) = F.U \ (F.L \ v[F.p])
 \(F::LU, B::AbstractMatrix) = F.U \ (F.L \ B[F.p,:])
+
+/(B::AbstractMatrix, F::LU) = @inbounds ((B/F.U)/F.L)[:,invperm(F.p)]

src/util.jl

@@ -101,3 +101,10 @@ TrivialView(a::AbstractArray{T,N}) where {T,N} = TrivialView{typeof(a),T,N}(a)
 # """
 @inline drop_sdims(a::StaticArrayLike) = TrivialView(a)
 @inline drop_sdims(a) = a
+
+Base.@propagate_inbounds function invperm(p::StaticVector)
+    # in difference to base, this does not check if p is a permutation (every value unique)
+     ip = similar(p)
+     ip[p] = 1:length(p)
+     similar_type(p)(ip)
+end

test/lu.jl

@@ -38,3 +38,15 @@ using StaticArrays, Test, LinearAlgebra
     l_u = l*u
     @test l*u ≈ a[p,:]
 end
+
+@testset "LU division ($m×$n)" for m in [1:4..., 15], n in [1:4..., 15]
+    a = SMatrix{m,m,Float64}(rand(Float64,m,m))
+    a_lu = lu(a)
+    b_col = SMatrix{m,n,Float64}(rand(Float64,m,n))
+    b_line = SMatrix{n,m,Float64}(rand(Float64,n,m))
+
+    # test if / and \ work with lu:
+    @test a\b_col ≈ a_lu\b_col
+    @test b_line/a ≈ b_line/a_lu
+
+end

test/runtests.jl

@@ -45,7 +45,7 @@ Random.seed!(44); include("eigen.jl")
 include("expm.jl")
 include("sqrtm.jl")
 include("lyap.jl")
-include("lu.jl")
+Random.seed!(42); include("lu.jl")
 Random.seed!(42); include("qr.jl")
 Random.seed!(42); include("chol.jl") # hermitian_type(::Type{Any}) for block algorithm
 include("deque.jl")