Refactor LU factorization with multiple implementation versions

- Introduce two LU factorization implementations (LUFactVersion1 and LUFactVersion2)
- Update abstract type names from LinearSystem* to Problem/SolutionMixin
- Modify kernel methods for both LU factorization versions
- Update test suite to cover both implementation versions
- Set LUFactorizationProblem as an alias for LUFactVersion2

Author: yangjunjie0320

src/SimpleLinearAlgebra.jl

@@ -2,10 +2,10 @@ module SimpleLinearAlgebra
 
     using LinearAlgebra
 
-    abstract type LinearSystemProblemMixin end
-    abstract type LinearSystemSolutionMixin end
+    abstract type ProblemMixin end
+    abstract type SolutionMixin end
 
-    function kernel(prob::LinearSystemProblemMixin)
+    function kernel(prob::ProblemMixin)
         MethodError(kernel, (typeof(prob),)) |> throw
     end
 

src/back_substitution.jl

@@ -1,4 +1,4 @@
-struct BackSubstitutionProblem <: LinearSystemProblemMixin
+struct BackSubstitutionProblem <: ProblemMixin
     u::AbstractMatrix 
     b::AbstractVector
     tol::Real
@@ -11,7 +11,7 @@ struct BackSubstitutionProblem <: LinearSystemProblemMixin
     end
 end
 
-struct BackSubstitutionSolution <: LinearSystemSolutionMixin
+struct BackSubstitutionSolution <: SolutionMixin
     x::AbstractVector
 end
 

src/forward_substitution.jl

@@ -1,4 +1,4 @@
-struct ForwardSubstitutionProblem <: LinearSystemProblemMixin
+struct ForwardSubstitutionProblem <: ProblemMixin
     l::AbstractMatrix
     b::AbstractVector
     tol::Real
@@ -13,7 +13,7 @@ end
 
 ForwardSubstitution = ForwardSubstitutionProblem
 
-struct ForwardSubstitutionSolution <: LinearSystemSolutionMixin
+struct ForwardSubstitutionSolution <: SolutionMixin
     x::AbstractVector
 end
 

src/lufact.jl

@@ -1,9 +1,4 @@
-struct LUFactorizationProblem <: LinearSystemProblemMixin
-    a::AbstractMatrix
-    tol::Real
-end
-
-struct LUFactorizationSolution <: LinearSystemSolutionMixin
+struct LUFactorizationSolution <: SolutionMixin
     l::AbstractMatrix
     u::AbstractMatrix
 
@@ -20,54 +15,72 @@ struct LUFactorizationSolution <: LinearSystemSolutionMixin
     end
 end
 
-LUFactorization = LUFactorizationProblem
+struct LUFactVersion1 <: ProblemMixin
+    a::AbstractMatrix
+    tol::Real
+end
 
-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
+function kernel(p::LUFactVersion1)
+    tol = p.tol
+    u = copy(p.a)
+    n = size(u, 1)
 
-    n = size(a, 1)
-    @assert size(a, 2) == n
-    @assert 1 <= k <= n-1
+    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
 
-    m = identity_matrix(T, n)
-    for i in k+1:n
-        m[i, k] = -a[i, k] / a[k, k]
+        m_k = zero(u)
+        m_k[1:n+1:end] .= 1
+        for i in k+1:n
+            m_k[i, k] = -u[i, k] / u[k, k]
+        end
+
+        l = l * inv(m_k)
+        u .= m_k * u
     end
-    return m
+
+    return LUFactorizationSolution(tril(l), triu(u))
 end
 
-function kernel(prob::LUFactorizationProblem)
-    u = copy(prob.a)
-    n = size(u, 1)
-    @assert size(u, 2) == n
 
-    tol = prob.tol
 
-    l = zero(u)
+struct LUFactVersion2 <: ProblemMixin
+    a::AbstractMatrix
+    tol::Real
+end
 
-    l[1:n+1:end] .= 1
+function kernel(p::LUFactVersion2)
+    tol = p.tol
+    u = copy(p.a)
+    n = size(u, 1)
 
-    for k in 1:n-1
+    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
+        for i in 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
+            for i=k+1:n
                 u[i, j] -= l[i, k] * u[k, j]
             end
         end
     end
 
-    return LUFactorizationSolution(l, triu(u))
+    return LUFactorizationSolution(tril(l), triu(u))
 end
 
-export LUFactorization, kernel
+export LUFactVersion1, LUFactVersion2
+
+LUFactorizationProblem = LUFactVersion2
+export LUFactorizationProblem

test/lufact.jl

@@ -1,9 +1,39 @@
-@testset "LU factorization" begin
+@testset "LU factorization version 1" begin
+    n = 4
+    tol = 1e-10
+    a = rand(n, n)
+
+    p = LUFactVersion1(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
+
+@testset "LU factorization version 2" begin
     n = 10
     tol = 1e-10
     a = rand(n, n)
 
-    p = LUFactorization(a, tol)
+    p = LUFactVersion2(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
+
+@testset "LU factorization" begin
+    n = 20
+    tol = 1e-10
+    a = rand(n, n)
+
+    p = LUFactorizationProblem(a, tol)
     s = kernel(p)
     l = s.l
     u = s.u