Do not use zero(...) in generic_lu! algorithm to support Unitful (#38659

) * make generic_lu! accept unitful matrices * Add Furlong test for lu
JuliaLang · Dec 5, 2020 · 65382c7 · 65382c7 · nanosoldier · Dec 6, 2020
1 parent 8629e21
commit 65382c7
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 3 deletions.
diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl
@@ -139,9 +139,9 @@ function generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot} = Val(true);
         for k = 1:minmn
             # find index max
             kp = k
-            if Pivot
-                amax = abs(zero(T))
-                for i = k:m
+            if Pivot && k < m
+                amax = abs(A[k, k])
+                for i = k+1:m
                     absi = abs(A[i,k])
                     if absi > amax
                         kp = i

diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl
@@ -324,6 +324,22 @@ include("trickyarithmetic.jl")
     @test B isa LinearAlgebra.LU{ElT,Matrix{ElT}}
 end
 
+# dimensional correctness:
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
+using .Main.Furlongs
+
+@testset "lu factorization with dimension type" begin
+    n = 4
+    A = Matrix(Furlong(1.0) * I, n, n)
+    F = lu(A).factors
+    @test Diagonal(F) == Diagonal(A)
+    # upper triangular part has a unit Furlong{1}
+    @test all(x -> typeof(x) == Furlong{1, Float64}, F[i,j] for j=1:n for i=1:j)
+    # lower triangular part is unitless Furlong{0}
+    @test all(x -> typeof(x) == Furlong{0, Float64}, F[i,j] for j=1:n for i=j+1:n)
+end
+
 @testset "Issue #30917. Determinant of integer matrix" begin
     @test det([1 1 0 0 1 0 0 0
                1 0 1 0 0 1 0 0