LU-based solve for matrix sizes 4:14

Author: schmrlng

src/solve.jl

@@ -1,7 +1,4 @@
-@inline (\)(a::StaticMatrix, b::StaticVector) = solve(Size(a), Size(b), a, b)
-
-# TODO: Ineffecient but requires some infrastructure (e.g. LU or QR) to make efficient so we fall back on inv for now
-@inline solve(::Size, ::Size, a, b) = inv(a) * b
+@inline (\)(a::StaticMatrix, b::StaticVecOrMat) = solve(Size(a), Size(b), a, b)
 
 @inline function solve(::Size{(1,1)}, ::Size{(1,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
     @inbounds return similar_type(b, typeof(a[1] \ b[1]))(a[1] \ b[1])
@@ -28,3 +25,23 @@ end
             (a[1,2]*a[3,1] - a[1,1]*a[3,2])*b[2] +
             (a[1,1]*a[2,2] - a[1,2]*a[2,1])*b[3]) / d )
 end
+
+@generated function solve(::Size{Sa}, ::Size{Sb}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVecOrMat{Tb}) where {Sa, Sb, Ta, Tb}
+    if Sa[end] != Sb[1]
+        throw(DimensionMismatch("right hand side B needs first dimension of size $(Sa[end]), has size $Sb"))
+    end
+    LinearAlgebra.checksquare(a)
+    if prod(Sa) ≤ 14*14
+        quote
+            @_inline_meta
+            L, U, p = lu(a)
+            U \ (L \ $(length(Sb) > 1 ? :(b[p,:]) : :(b[p])))
+        end
+    else
+        quote
+            @_inline_meta
+            T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/one(Ta))
+            similar_type(b, T)(Matrix(a) \ b)
+        end
+    end
+end

test/solve.jl

@@ -1,7 +1,7 @@
 using StaticArrays, Compat.Test
 
 @testset "Solving linear system" begin
-    @testset "Problem size: $n x $n. Matrix type: $m. Element type: $elty" for n in (1,2,3,4),
+    @testset "Problem size: $n x $n. Matrix type: $m. Element type: $elty" for n in (1,2,3,4,5,14,15,50),
             (m, v) in ((SMatrix{n,n}, SVector{n}), (MMatrix{n,n}, MVector{n})),
                 elty in (Float64, Int)
 
@@ -16,3 +16,14 @@ using StaticArrays, Compat.Test
     @test_throws DimensionMismatch m1\v
     @test_throws DimensionMismatch m1\m2
 end
+
+@testset "Solving linear system (multiple RHS)" begin
+    @testset "Problem size: $n x $n. Matrix type: $m1. Element type: $elty" for n in (1,2,3,4,5,14,15,50),
+            (m1, m2) in ((SMatrix{n,n}, SMatrix{n,2}), (MMatrix{n,n}, MMatrix{n,2})),
+                elty in (Float64, Int)
+
+        A = elty.(rand(-99:2:99, n, n))
+        b = A * elty.(rand(2:5, n, 2))
+        @test m1(A)\m2(b) ≈ A\b
+    end
+end