fixes

Author: Shreyas-Ekanathan

Project.toml

@@ -10,7 +10,6 @@ Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SparseBandedMatrices = "bd59d7e1-4699-4102-944e-d05209cb92aa"
 StrideArraysCore = "7792a7ef-975c-4747-a70f-980b88e8d1da"
-TriangularSolve = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf"
 VectorizedRNG = "33b4df10-0173-11e9-2a0c-851a7edac40e"
 
 [compat]
@@ -20,7 +19,6 @@ Polyester = "0.3.2,0.4.1, 0.5, 0.6, 0.7"
 PrecompileTools = "1"
 SparseBandedMatrices = "1"
 StrideArraysCore = "0.5.5"
-TriangularSolve = "0.2"
 VectorizedRNG = "0.2.25"
 julia = "1.10"
 

src/butterflylu.jl

@@ -20,7 +20,6 @@ end
 struct 🦋workspace{T}
     A::Matrix{T}
     b::Vector{T}
-    B::Matrix{T}
     ws::Vector{T}
     U::Matrix{T}
     V::Matrix{T}
@@ -33,18 +32,17 @@ struct 🦋workspace{T}
             xn = 4 - M % 4
             b = [b; rand(xn)]
         end
-        B = similar(A)
         U, V = (similar(A), similar(A))
         ws = 🦋generate_random!(A)
-        new{eltype(A)}(A, b, B, ws, U, V, out)
+        materializeUV(U, V, ws)
+        new{eltype(A)}(A, b, ws, U, V, out)
     end
 end
 
 function 🦋lu!(workspace::🦋workspace, M, thread)
-    (;A, b, B, ws, U, V, out) = workspace
-    🦋mul!(copyto!(B, A), ws)
-    materializeUV(U, V, ws)
-    F = RecursiveFactorization.lu!(B, Val(false), thread)
+    (;A, b, ws, U, V, out) = workspace
+    🦋mul!(A, ws)
+    F = RecursiveFactorization.lu!(A, Val(false), thread)
     sol = V * (F \ (U' * b))  
     out .= @view sol[1:M]  
     out
@@ -55,14 +53,14 @@ const butterfly_workspace = 🦋workspace;
 function 🦋mul_level!(A, u, v)
     M, N = size(A)
     @assert M == length(u) && N == length(v)
-    Mh = M >>> 1
-    Nh = N >>> 1
-    @turbo for n in 1 : Nh
-        for m in 1 : Mh
+    M_half = M >>> 1
+    N_half = N >>> 1
+    @turbo for n in 1 : N_half
+        for m in 1 : M_half
             A11 = A[m, n]
-            A21 = A[m + Mh, n]
-            A12 = A[m, n + Nh]
-            A22 = A[m + Mh, n + Nh]
+            A21 = A[m + M_half, n]
+            A12 = A[m, n + N_half]
+            A22 = A[m + M_half, n + N_half]
 
             T1 = A11 + A12
             T2 = A21 + A22
@@ -74,32 +72,32 @@ function 🦋mul_level!(A, u, v)
             C22 = T3 - T4
 
             u1 = u[m]
-            u2 = u[m + Mh]
+            u2 = u[m + M_half]
             v1 = v[n]
-            v2 = v[n + Nh]
+            v2 = v[n + N_half]
 
             A[m, n] = u1 * C11 * v1
-            A[m + Mh, n] = u2 * C21 * v1
-            A[m, n + Nh] = u1 * C12 * v2
-            A[m + Mh, n + Nh] = u2 * C22 * v2
+            A[m + M_half, n] = u2 * C21 * v1
+            A[m, n + N_half] = u1 * C12 * v2
+            A[m + M_half, n + N_half] = u2 * C22 * v2
         end
     end 
 end
 
 function 🦋mul!(A, uv)
     M, N = size(A)
     @assert M == N
-    Mh = M >>> 1
+    M_half = M >>> 1
 
-    U₁ = @view(uv[1:Mh]) 
-    V₁ = @view(uv[(Mh + 1):(M)]) 
-    U₂ = @view(uv[(1 + M):(M + Mh)]) 
-    V₂ = @view(uv[(1 + M + Mh):(2 * M)]) 
+    U₁ = @view(uv[1:M_half]) 
+    V₁ = @view(uv[(M_half + 1):(M)]) 
+    U₂ = @view(uv[(1 + M):(M + M_half)]) 
+    V₂ = @view(uv[(1 + M + M_half):(2 * M)]) 
 
-    🦋mul_level!(@view(A[1:Mh, 1:Mh]), U₁, V₁) 
-    🦋mul_level!(@view(A[Mh + 1:M, 1:Mh]), U₂, V₁) 
-    🦋mul_level!(@view(A[1:Mh, Mh + 1:M]), U₁, V₂) 
-    🦋mul_level!(@view(A[Mh + 1:M, Mh + 1:M]), U₂, V₂) 
+    🦋mul_level!(@view(A[1:M_half, 1:M_half]), U₁, V₁) 
+    🦋mul_level!(@view(A[M_half + 1:M, 1:M_half]), U₂, V₁) 
+    🦋mul_level!(@view(A[1:M_half, M_half + 1:M]), U₁, V₂) 
+    🦋mul_level!(@view(A[M_half + 1:M, M_half + 1:M]), U₂, V₂) 
 
     U = @view(uv[(1 + 2 * M):(3 * M)]) 
     V = @view(uv[(1 + 3 * M):(4 * M)]) 
@@ -128,7 +126,7 @@ function 🦋!(C::SparseBandedMatrix, A::Diagonal, B::Diagonal)
     C
 end
 
-function 🦋2!(C, A::Diagonal, B::Diagonal)
+function 🦋!(C, A::Diagonal, B::Diagonal)
     @assert size(A) == size(B)
     A1 = size(A, 1)
 
@@ -144,27 +142,27 @@ end
 
 function materializeUV(U, V, uv)
     M = size(U, 1)
-    Mh = M >>> 1   
+    M_half = M >>> 1   
 
-    U₁u, U₁l = diagnegbottom(@view(uv[1:Mh])) #Mh
-    U₂u, U₂l = diagnegbottom(@view(uv[(1 + 2 * Mh):(M + Mh)])) #Mh
-    V₁u, V₁l = diagnegbottom(@view(uv[(Mh + 1):(2 * Mh)])) #Mh
-    V₂u, V₂l = diagnegbottom(@view(uv[(1 + 3 * Mh):(2 * Mh + M)])) #Mh
+    U₁u, U₁l = diagnegbottom(@view(uv[1:M_half])) #M_half
+    U₂u, U₂l = diagnegbottom(@view(uv[(1 + 2 * M_half):(M + M_half)])) #M_half
+    V₁u, V₁l = diagnegbottom(@view(uv[(M_half + 1):(2 * M_half)])) #M_half
+    V₂u, V₂l = diagnegbottom(@view(uv[(1 + 3 * M_half):(2 * M_half + M)])) #M_half
     Uu, Ul = diagnegbottom(@view(uv[(1 + 2 * M):(3 * M)])) #M
     Vu, Vl = diagnegbottom(@view(uv[(1 + 3 * M):(4 * M)])) #M
 
     Bu2 = SparseBandedMatrix{typeof(uv[1])}(undef, M, M)
 
-    🦋2!(view(Bu2, 1 : Mh, 1 : Mh), U₁u, U₁l)
-    🦋2!(view(Bu2, Mh + 1: M, Mh + 1: M), U₂u, U₂l)
+    🦋!(view(Bu2, 1 : M_half, 1 : M_half), U₁u, U₁l)
+    🦋!(view(Bu2, M_half + 1: M, M_half + 1: M), U₂u, U₂l)
 
     Bu1 = SparseBandedMatrix{typeof(uv[1])}(undef, M, M)
     🦋!(Bu1, Uu, Ul)
 
     Bv2 = SparseBandedMatrix{typeof(uv[1])}(undef, M, M)
 
-    🦋2!(view(Bv2, 1 : Mh, 1 : Mh), V₁u, V₁l)
-    🦋2!(view(Bv2, Mh + 1: M, Mh + 1: M), V₂u, V₂l)
+    🦋!(view(Bv2, 1 : M_half, 1 : M_half), V₁u, V₁l)
+    🦋!(view(Bv2, M_half + 1: M, M_half + 1: M), V₂u, V₂l)
 
     Bv1 = SparseBandedMatrix{typeof(uv[1])}(undef, M, M)
     🦋!(Bv1, Vu, Vl)

test/runtests.jl

@@ -81,8 +81,9 @@ end
     for i in 790 : 810
         A = wilkinson(i)
         b = rand(i)
-        ws = RecursiveFactorization.🦋workspace(A, b)    
+        ws = RecursiveFactorization.🦋workspace(copy(A), copy(b))    
         out = RecursiveFactorization.🦋lu!(ws, i, Val(true))
         @test norm(A * out .- b) <= 1e-10
     end
 end
+