Do not use zero(...) in generic_lu! algorithm to support Unitful (Jul…

…iaLang#38659)

* make generic_lu! accept unitful matrices

* Add Furlong test for lu

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

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