Merge pull request #3 from YingboMa/myb/generic

Make LU more generic

Author: YingboMa

.gitignore

@@ -1,3 +1,4 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+Manifest.toml

Project.toml

@@ -0,0 +1,17 @@
+name = "RecursiveFactorization"
+uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
+authors = ["Yingbo Ma <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[compat]
+julia = "1"
+
+[extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test", "Random"]

src/lu.jl

@@ -1,20 +1,23 @@
-using LinearAlgebra: BlasInt, LU, UnitLowerTriangular, ldiv!, BLAS, checknonsingular
+using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, ldiv!, BLAS, checknonsingular
 
-lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true);
-    check::Bool = true, blocksize::Integer = 16) = lu!(copy(A), pivot;
-                                                        check = check, blocksize = blocksize)
+function lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true);
+            check::Bool = true, blocksize::Integer = 16)
+    lu!(copy(A), pivot; check = check, blocksize = blocksize)
+end
 
-lu!(A, pivot::Union{Val{false}, Val{true}} = Val(true);
-    check::Bool = true, blocksize::Integer = 16) = lu!(copy(A), Vector{BlasInt}(undef, min(size(A)...)), pivot;
-                                                        check = check, blocksize = blocksize)
+function lu!(A, pivot::Union{Val{false}, Val{true}} = Val(true);
+             check::Bool = true, blocksize::Integer = 16)
+    lu!(A, Vector{BlasInt}(undef, min(size(A)...)), pivot;
+        check = check, blocksize = blocksize)
+end
 
 function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
              pivot::Union{Val{false}, Val{true}} = Val(true);
              check::Bool=true, blocksize::Integer=16) where T
     info = Ref(zero(BlasInt))
     m, n = size(A)
     mnmin = min(m, n)
-    if isconcretetype(T)
+    if T <: BlasFloat && A isa StridedArray
         reckernel!(A, pivot, m, mnmin, ipiv, info, blocksize)
         if m < n # fat matrix
             # [AL AR]
@@ -115,15 +118,15 @@ end
     Modified from https://github.com/JuliaLang/julia/blob/b56a9f07948255dfbe804eef25bdbada06ec2a57/stdlib/LinearAlgebra/src/lu.jl
     License is MIT: https://julialang.org/license
 =#
-function _generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot}, ipiv, info) where {Pivot,T}
+function _generic_lufact!(A, ::Val{Pivot}, ipiv, info) where Pivot
     m, n = size(A)
     minmn = length(ipiv)
     @inbounds begin
         for k = 1:minmn
             # find index max
             kp = k
             if Pivot
-              amax = abs(zero(T))
+              amax = abs(zero(eltype(A)))
               for i = k:m
                   absi = abs(A[i,k])
                   if absi > amax

test/runtests.jl

@@ -6,8 +6,8 @@ using  Random
 
 Random.seed!(12)
 
-baselu = LinearAlgebra.lu
-mylu = RecursiveFactorization.lu
+const baselu = LinearAlgebra.lu
+const mylu = RecursiveFactorization.lu
 
 function testlu(A, MF, BF)
     @test MF.info == BF.info