Implement Butterfly Factorization method using RecursiveFactorization.jl

benchmarks/lu.jl

@@ -1,7 +1,8 @@
 using BenchmarkTools, Random, VectorizationBase
 using LinearAlgebra, LinearSolve, MKL_jll
+using RecursiveFactorization
+
 nc = min(Int(VectorizationBase.num_cores()), Threads.nthreads())
-BLAS.set_num_threads(nc)
 BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
 
 function luflop(m, n = m; innerflop = 2)
@@ -24,10 +25,10 @@ algs = [
     RFLUFactorization(),
     MKLLUFactorization(),
     FastLUFactorization(),
-    SimpleLUFactorization()
+    SimpleLUFactorization(),
+    ButterflyFactorization(Val(true))
 ]
 res = [Float64[] for i in 1:length(algs)]
-
 ns = 4:8:500
 for i in 1:length(ns)
     n = ns[i]
@@ -65,3 +66,4 @@ p
 
 savefig("lubench.png")
 savefig("lubench.pdf")
+

ext/LinearSolveRecursiveFactorizationExt.jl

@@ -1,7 +1,7 @@
 module LinearSolveRecursiveFactorizationExt
 
 using LinearSolve: LinearSolve, userecursivefactorization, LinearCache, @get_cacheval,
-                   RFLUFactorization, RF32MixedLUFactorization, default_alias_A,
+                   RFLUFactorization, ButterflyFactorization, RF32MixedLUFactorization, default_alias_A,
                    default_alias_b
 using LinearSolve.LinearAlgebra, LinearSolve.ArrayInterface, RecursiveFactorization
 using SciMLBase: SciMLBase, ReturnCode
@@ -19,7 +19,6 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::RFLUFactorization
         end
         fact = RecursiveFactorization.lu!(A, ipiv, Val(P), Val(T), check = false)
         cache.cacheval = (fact, ipiv)
-
         if !LinearAlgebra.issuccess(fact)
             return SciMLBase.build_linear_solution(
                 alg, cache.u, nothing, cache; retcode = ReturnCode.Failure)
@@ -105,4 +104,32 @@ function SciMLBase.solve!(
         alg, cache.u, nothing, cache; retcode = ReturnCode.Success)
 end
 
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::ButterflyFactorization;
+        kwargs...)
+    A = cache.A
+    A = convert(AbstractMatrix, A)
+    b = cache.b
+    M, N = size(A)
+    B, U, V = cache.cacheval[2], cache.cacheval[3], cache.cacheval[4]
+    if cache.isfresh
+        @assert M==N "A must be square"
+        U, V, F, out = RecursiveFactorization.🦋workspace(A, b, B, U, V, alg.thread)
+        cache.cacheval = (A, B, U, V, F)
+        cache.isfresh = false
+        if (M % 4 != 0)
+            b = [b; rand(4 - M % 4)]
+        end
+    end
+    A, B, U, V, F = cache.cacheval
+    sol = V * (F \ (U * b))
+    out .= @view sol[1:M]    
+   SciMLBase.build_linear_solution(alg, out, nothing, cache)
+end
+
+function LinearSolve.init_cacheval(alg::ButterflyFactorization, A, b, u, Pl, Pr, maxiters::Int,
+        abstol, reltol, verbose::Bool, assumptions::LinearSolve.OperatorAssumptions)
+    A, A, A', A, RecursiveFactorization.lu!(rand(1, 1), alg.thread)
 end
+
+end
+

src/LinearSolve.jl

@@ -422,7 +422,7 @@ for kralg in (Krylov.lsmr!, Krylov.craigmr!)
 end
 for alg in (:LUFactorization, :FastLUFactorization, :SVDFactorization,
     :GenericFactorization, :GenericLUFactorization, :SimpleLUFactorization,
-    :RFLUFactorization, :UMFPACKFactorization, :KLUFactorization, :SparspakFactorization,
+    :RFLUFactorization, :ButterflyFactorization, :UMFPACKFactorization, :KLUFactorization, :SparspakFactorization,
     :DiagonalFactorization, :CholeskyFactorization, :BunchKaufmanFactorization,
     :CHOLMODFactorization, :LDLtFactorization, :AppleAccelerateLUFactorization,
     :MKLLUFactorization, :MetalLUFactorization, :CUSOLVERRFFactorization)
@@ -464,7 +464,7 @@ cudss_loaded(A) = false
 is_cusparse(A) = false
 
 export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
-       GenericLUFactorization, SimpleLUFactorization, RFLUFactorization,
+       GenericLUFactorization, SimpleLUFactorization, RFLUFactorization, ButterflyFactorization,
        NormalCholeskyFactorization, NormalBunchKaufmanFactorization,
        UMFPACKFactorization, KLUFactorization, FastLUFactorization, FastQRFactorization,
        SparspakFactorization, DiagonalFactorization, CholeskyFactorization,

src/extension_algs.jl

@@ -254,6 +254,29 @@ function RFLUFactorization(; pivot = Val(true), thread = Val(true), throwerror =
     RFLUFactorization(pivot, thread; throwerror)
 end
 
+"""
+`ButterflyFactorization()`
+
+A fast pure Julia LU-factorization implementation
+using RecursiveFactorization.jl. This method utilizes a butterly 
+factorization approach rather than pivoting. 
+"""
+struct ButterflyFactorization{T} <: AbstractDenseFactorization
+    thread::Val{T}
+    function ButterflyFactorization(::Val{T}; throwerror = true) where {T}
+        if !userecursivefactorization(nothing)
+            throwerror &&
+                error("ButterflyFactorization requires that RecursiveFactorization.jl is loaded, i.e. `using RecursiveFactorization`")
+        end
+        new{T}()
+    end
+end
+
+function ButterflyFactorization(; thread = Val(true), throwerror = true)
+    ButterflyFactorization(thread; throwerror)
+end
+
+
 # There's no options like pivot here.
 # But I'm not sure it makes sense as a GenericFactorization
 # since it just uses `LAPACK.getrf!`.

test/butterfly.jl

@@ -0,0 +1,35 @@
+using LinearAlgebra, LinearSolve
+using Test
+using RecursiveFactorization
+
+@testset "Random Matricies" begin
+    for i in 490 : 510
+        A = rand(i, i)
+        b = rand(i)
+        prob = LinearProblem(A, b)
+        x = solve(prob, ButterflyFactorization())
+        @test norm(A * x .- b) <= 1e-4
+    end
+end
+
+function wilkinson(N)
+    A = zeros(N, N)
+    A[1:(N+1):N*N] .= 1
+    A[:, end] .= 1
+    for n in 1:(N - 1)
+        for r in (n + 1):N
+            @inbounds A[r, n] = -1
+        end
+    end
+    A
+end
+
+@testset "Wilkinson" begin
+    for i in 790 : 810
+        A = wilkinson(i)
+        b = rand(i)
+        prob = LinearProblem(A, b)
+        x = solve(prob, ButterflyFactorization())
+        @test norm(A * x .- b) <= 1e-10
+    end
+end