Add LU factorization implementation

- Include LU factorization module in SimpleLinearAlgebra
- Update Makefile to use wildcard patterns for source and test files
- Add LU factorization test suite to test runner

Author: yangjunjie0320

Makefile

@@ -1,20 +1,20 @@
 jl:=julia
 
-init: src/SimpleLinearAlgebra.jl
+init: src/*jl
 	$(jl) -e 'd=pwd(); @assert isdir(d);\
 	using Pkg: activate, instantiate, develop, precompile; \
 	activate(d); instantiate(); activate(joinpath(d, "examples")); \
 	develop(path=d); instantiate(); precompile();'
 	@echo "Environment initialized at: $$PWD"
 
-update: src/SimpleLinearAlgebra.jl
+update: src/*jl
 	$(jl) -e 'd=pwd(); @assert isdir(d);\
 	using Pkg: activate, update, develop, precompile; \
 	activate(d); update(); activate(joinpath(d, "examples")); \
 	update(); precompile();'
 	@echo "Environment updated at: $$PWD"
 
-test: test/runtests.jl
+test: test/*jl
 	@echo "Running tests at: $$PWD"
 	$(jl) -e 'd=pwd(); @assert isdir(d);\
 	using Pkg: activate, test; \

src/SimpleLinearAlgebra.jl

@@ -11,4 +11,5 @@ module SimpleLinearAlgebra
 
     include("forward_substitution.jl")
     include("back_substitution.jl")
+    include("lufact.jl")
 end

src/lufact.jl

@@ -0,0 +1,73 @@
+struct LUFactorizationProblem <: LinearSystemProblemMixin
+    a::AbstractMatrix
+    tol::Real
+end
+
+struct LUFactorizationSolution <: LinearSystemSolutionMixin
+    l::AbstractMatrix
+    u::AbstractMatrix
+
+    function LUFactorizationSolution(l, u)
+        if !istril(l)
+            ArgumentError("Matrix must be lower triangular") |> throw
+        end
+
+        if !istriu(u)
+            ArgumentError("Matrix must be upper triangular") |> throw
+        end
+
+        return new(l, u)
+    end
+end
+
+LUFactorization = LUFactorizationProblem
+
+function elemination_step(a::AbstractMatrix{T}, k::Integer) where T
+    # given a matrix A, perform the k-th elimination step.
+    # M A = N, which M is a lower triangular matrix with ones on the diagonal,
+    # and N is the matrix A with the k-th column eliminated.
+    # The inverse of M is simply given by Minv = 2I - m
+
+    n = size(a, 1)
+    @assert size(a, 2) == n
+    @assert 1 <= k <= n-1
+
+    m = identity_matrix(T, n)
+    for i in k+1:n
+        m[i, k] = -a[i, k] / a[k, k]
+    end
+    return m
+end
+
+function kernel(prob::LUFactorizationProblem)
+    u = copy(prob.a)
+    n = size(u, 1)
+    @assert size(u, 2) == n
+
+    tol = prob.tol
+
+    l = zero(u)
+
+    l[1:n+1:end] .= 1
+
+    for k in 1:n-1
+
+        if abs(u[k, k]) < tol
+            ArgumentError("Gaussian elimination failed") |> throw
+        end
+
+        for i=k+1:n
+            l[i, k] = u[i, k] / u[k, k]
+        end
+
+        for j in k+1:n
+            for i in k+1:n
+                u[i, j] -= l[i, k] * u[k, j]
+            end
+        end
+    end
+
+    return LUFactorizationSolution(l, triu(u))
+end
+
+export LUFactorization, kernel

test/lufact.jl

@@ -0,0 +1,14 @@
+@testset "LU factorization" begin
+    n = 10
+    tol = 1e-10
+    a = rand(n, n)
+
+    p = LUFactorization(a, tol)
+    s = kernel(p)
+    l = s.l
+    u = s.u
+
+    @test istril(l)
+    @test istriu(u)
+    @test maximum(abs, a - l * u) < tol
+end

test/runtests.jl

@@ -8,3 +8,7 @@ end
 @testset "backward substitution" begin
     include("backward_substitution.jl")
 end
+
+@testset "LU factorization" begin
+    include("lufact.jl")
+end