Fix A\B where A is nonsquare

Fixes #606

Author: Chris Foster

src/solve.jl

@@ -1,17 +1,17 @@
-@inline (\)(a::StaticMatrix, b::StaticVecOrMat) = solve(Size(a), Size(b), 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}
+@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])
 end
 
-@inline function solve(::Size{(2,2)}, ::Size{(2,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
+@inline function _solve(::Size{(2,2)}, ::Size{(2,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
     d = det(a)
     T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/d)
     @inbounds return similar_type(b, T)((a[2,2]*b[1] - a[1,2]*b[2])/d,
                                         (a[1,1]*b[2] - a[2,1]*b[1])/d)
 end
 
-@inline function solve(::Size{(3,3)}, ::Size{(3,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
+@inline function _solve(::Size{(3,3)}, ::Size{(3,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
     d = det(a)
     T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/d)
     @inbounds return similar_type(b, T)(
@@ -28,12 +28,12 @@ end
 
 for Sa in [(2,2), (3,3)]  # not needed for Sa = (1, 1);
     @eval begin
-        @inline function solve(::Size{$Sa}, ::Size{Sb}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticMatrix{<:Any, <:Any, Tb}) where {Sb, Ta, Tb}
+        @inline function _solve(::Size{$Sa}, ::Size{Sb}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticMatrix{<:Any, <:Any, Tb}) where {Sb, Ta, Tb}
             d = det(a)
             T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/d)
             c = similar(b, T)
             for col = 1:Sb[2]
-                @inbounds c[:, col] = solve(Size($Sa), Size($Sa[1],), a, b[:, col])
+                @inbounds c[:, col] = _solve(Size($Sa), Size($Sa[1],), a, b[:, col])
             end
             return similar_type(b, T)(c)
         end
@@ -42,22 +42,28 @@ 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"))
+@generated function _solve(::Size{Sa}, ::Size{Sb}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVecOrMat{Tb}) where {Sa, Sb, Ta, Tb}
+    if Sa[1] != Sb[1]
+        return quote
+            throw(DimensionMismatch("Left and right hand side first dimensions do not match in backdivide (got sizes $Sa and $Sb)"))
+        end
     end
-    LinearAlgebra.checksquare(a)
-    if prod(Sa) ≤ 14*14
+    if prod(Sa) ≤ 14*14 && Sa[1] == Sa[2]
+        # TODO: Consider triangular special cases as in Base?
         quote
             @_inline_meta
             LUp = lu(a)
             LUp.U \ (LUp.L \ $(length(Sb) > 1 ? :(b[LUp.p,:]) : :(b[LUp.p])))
         end
+        # TODO: Could also use static QR here if `a` is nonsquare.
+        # Requires that we implement \(::StaticArrays.QR,::StaticVecOrMat)
     else
+        # Fall back to LinearAlgebra, but carry across the statically known size.
+        outsize = length(Sb) == 1 ? Size(Sa[2]) : Size(Sa[2],Sb[end])
         quote
             @_inline_meta
             T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/one(Ta))
-            similar_type(b, T)(Matrix(a) \ b)
+            similar_type(b, T, $outsize)(Matrix(a) \ b)
         end
     end
 end

test/solve.jl

@@ -41,6 +41,17 @@ end
 
     end
 
+    # Solve with non-square left hand sides (#606)
+    m1 = @SMatrix[0.2 0.3
+                  0.0 0.1
+                  0.5 0.1]
+    m2 = @SVector[1,2,3]
+    @test @inferred(m1\m2) ≈ Array(m1)\Array(m2)
+    m2 = @SMatrix[1 4
+                  2 5
+                  3 6]
+    @test @inferred(m1\m2) ≈ Array(m1)\Array(m2)
+
     @testset "Mixed static/dynamic" begin
         m2 = @SMatrix([0.2 0.3; 0.0 0.1])
         for m1 in (@SMatrix([1.0 0; 0 1.0]), @SMatrix([1.0 0; 1.0 1.0]),