diff --git a/src/SparseBandedMatrices.jl b/src/SparseBandedMatrices.jl index 4c9cc56..5b71f8c 100644 --- a/src/SparseBandedMatrices.jl +++ b/src/SparseBandedMatrices.jl @@ -53,32 +53,32 @@ diagonal indices range from 1 to N+M-1. The storage is optimized for fast matrix multiplication operations used in butterfly factorizations. """ struct SparseBandedMatrix{T} <: AbstractMatrix{T} - size::Tuple{Int, Int} - indices::Vector{Int} - diags::Vector{Vector{T}} - function SparseBandedMatrix{T}(::UndefInitializer, N, M) where {T} + size :: Tuple{Int, Int} + indices :: Vector{Int} + diags :: Vector{Vector{T}} + function SparseBandedMatrix{T}(::UndefInitializer, N, M) where T size = (N, M) indices = Int[] diags = Vector{T}[] - return new(size, indices, diags) + new(size, indices, diags) end - function SparseBandedMatrix{T}(ind_vals, diag_vals, N, M) where {T} + function SparseBandedMatrix{T}(ind_vals, diag_vals, N, M) where T size = (N, M) perm = sortperm(ind_vals) indices = ind_vals[perm] - for i in 1:(length(indices) - 1) + for i in 1 : length(indices) - 1 @assert indices[i] != indices[i + 1] end diags = diag_vals[perm] - return new(size, indices, diags) + new(size, indices, diags) end end -function Base.size(M::SparseBandedMatrix) - return M.size +function Base.size(M :: SparseBandedMatrix) + M.size end -function Base.getindex(M::SparseBandedMatrix{T}, i::Int, j::Int, I::Int...) where {T} +function Base.getindex(M :: SparseBandedMatrix{T}, i :: Int, j :: Int, I :: Int...) where T @boundscheck checkbounds(M, i, j, I...) rows, cols = size(M) wanted_ind = rows - i + j @@ -90,11 +90,11 @@ function Base.getindex(M::SparseBandedMatrix{T}, i::Int, j::Int, I::Int...) wher return M.diags[ind][i] end end - return zero(T) + zero(T) end -function Base.setindex!(M::SparseBandedMatrix{T}, val, i::Int, j::Int, I::Int...) where {T} - @boundscheck checkbounds(M, i, j, I...) +function Base.setindex!(M :: SparseBandedMatrix{T}, val, i :: Int, j :: Int, I :: Int...) where T + @boundscheck checkbounds(M, i, j, I...) rows = size(M, 1) wanted_ind = rows - i + j ind = searchsortedfirst(M.indices, wanted_ind) @@ -107,8 +107,8 @@ function Base.setindex!(M::SparseBandedMatrix{T}, val, i::Int, j::Int, I::Int... else M.diags[ind][i] = val isa T ? val : convert(T, val)::T end - return val -end + val + end """ setdiagonal!(M::SparseBandedMatrix{T}, diagvals, lower::Bool) where T @@ -141,7 +141,7 @@ setdiagonal!(A, [1.0, 2.0], false) # Throws - `ErrorException`: If `length(diagvals) > rows` (diagonal is too large for the matrix) """ -function setdiagonal!(M::SparseBandedMatrix{T}, diagvals, lower::Bool) where {T} +function setdiagonal!(M :: SparseBandedMatrix{T}, diagvals, lower :: Bool) where T rows, cols = size(M) if length(diagvals) > rows error("size of diagonal is too big for the matrix") @@ -161,98 +161,85 @@ function setdiagonal!(M::SparseBandedMatrix{T}, diagvals, lower::Bool) where {T} M.diags[ind][i] = diagvals[i] isa T ? diagvals[i] : convert(T, diagvals[i])::T end end - return diagvals + diagvals end # C = Cb + aAB -function LinearAlgebra.mul!(C::Matrix{T}, A::SparseBandedMatrix{T}, B::Matrix{T}, a::Number, b::Number) where {T} +function LinearAlgebra.mul!(C :: Matrix{T}, A:: SparseBandedMatrix{T}, B :: Matrix{T}, a :: Number, b :: Number) where T @assert size(A, 2) == size(B, 1) @assert size(A, 1) == size(C, 1) @assert size(B, 2) == size(C, 2) - if iszero(b) - fill!(C, zero(T)) - else - C .*= b - end + C.*=b rows, cols = size(A) @inbounds for (ind, location) in enumerate(A.indices) @threads for i in 1:length(A.diags[ind]) # value: diag[i] - # index in array: + # index in array: # if ind < rows(A), then index = (rows - loc + i, i) # else index = (i, loc - cols + i) val = A.diags[ind][i] * a - if location < rows - index_i = rows - location + i - index_j = i + if location < rows + index_i = rows - location + i + index_j = i else - index_i = i - index_j = location - cols + i + index_i = i + index_j = location - cols + i end #A[index_i, index_j] * B[index_j, j] = C[index_i, j] - for j in 1:size(B, 2) - C[index_i, j] = fma(val, B[index_j, j], C[index_i, j]) + for j in 1 : size(B, 2) + C[index_i, j] = muladd(val, B[index_j, j], C[index_i, j]) end end end - return C + C end -# C = C*b + a*B*A -function LinearAlgebra.mul!(C::Matrix{T}, A::Matrix{T}, B::SparseBandedMatrix{T}, a::Number, b::Number) where {T} +# C = Cb + aBA +function LinearAlgebra.mul!(C :: Matrix{T}, A:: Matrix{T}, B :: SparseBandedMatrix{T}, a :: Number, b :: Number) where T @assert size(A, 2) == size(B, 1) @assert size(A, 1) == size(C, 1) @assert size(B, 2) == size(C, 2) - if iszero(b) - fill!(C, zero(T)) - else - C .*= b - end + C.*=b rows, cols = size(B) @inbounds for (ind, location) in enumerate(B.indices) @threads for i in eachindex(B.diags[ind]) - val = B.diags[ind][i] * a - if location < rows - index_i = rows - location + i - index_j = i - else - index_i = i - index_j = location - cols + i - end - @simd for j in 1:size(A, 1) - C[j, index_j] = fma(val, A[j, index_i], C[j, index_j]) + val = B.diags[ind][i] * a + if location < rows + index_i = rows - location + i + index_j = i + else + index_i = i + index_j = location - cols + i + end + @simd for j in 1 : size(A, 1) + C[j, index_j] = muladd(val, A[j, index_i], C[j, index_j]) end end end - return C + C end -function LinearAlgebra.mul!(C::SparseBandedMatrix{T}, A::SparseBandedMatrix{T}, B::SparseBandedMatrix{T}, a::Number, b::Number) where {T} +function LinearAlgebra.mul!(C :: SparseBandedMatrix{T}, A:: SparseBandedMatrix{T}, B :: SparseBandedMatrix{T}, a :: Number, b :: Number) where T @assert size(A, 2) == size(B, 1) @assert size(A, 1) == size(C, 1) @assert size(B, 2) == size(C, 2) - if iszero(b) - empty!(C.indices) - empty!(C.diags) - else - C .*= b - end + C.*=b rows_a, cols_a = size(A) rows_b, cols_b = size(B) @inbounds for (ind_a, location_a) in enumerate(A.indices) @threads for i in eachindex(A.diags[ind_a]) val_a = A.diags[ind_a][i] * a - if location_a < rows_a - index_ia = rows_a - location_a + i - index_ja = i + if location_a < rows_a + index_ia = rows_a - location_a + i + index_ja = i else - index_ia = i - index_ja = location_a - cols_a + i + index_ia = i + index_ja = location_a - cols_a + i end min_loc = rows_b - index_ja + 1 max_loc = 2 * rows_b - index_ja @@ -261,46 +248,42 @@ function LinearAlgebra.mul!(C::SparseBandedMatrix{T}, A::SparseBandedMatrix{T}, # if ind < rows(A), then index = (rows - loc + i, i) #rows - loc + j = index_ja, j = index_ja - rows + loc # else index = (i, loc - cols + i) - # if location < rows(B), then + # if location < rows(B), then if location_b <= rows_b && location_b >= min_loc j = index_ja - rows_b + location_b index_jb = j val_b = B.diags[ind_b][j] - C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) + C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) elseif location_b > rows_b && location_b <= max_loc j = index_ja - index_jb = location_b - cols_b + j + index_jb = location_b - cols_b + j val_b = B.diags[ind_b][j] - C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) - end + C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) + end end end end - return C + C end -function LinearAlgebra.mul!(C::Matrix{T}, A::SparseBandedMatrix{T}, B::SparseBandedMatrix{T}, a::Number, b::Number) where {T} +function LinearAlgebra.mul!(C :: Matrix{T}, A:: SparseBandedMatrix{T}, B :: SparseBandedMatrix{T}, a :: Number, b :: Number) where T @assert size(A, 2) == size(B, 1) @assert size(A, 1) == size(C, 1) @assert size(B, 2) == size(C, 2) - if iszero(b) - fill!(C, zero(T)) - else - C .*= b - end + C.*=b rows_a, cols_a = size(A) rows_b, cols_b = size(B) @inbounds for (ind_a, location_a) in enumerate(A.indices) @threads for i in eachindex(A.diags[ind_a]) val_a = A.diags[ind_a][i] * a - if location_a < rows_a - index_ia = rows_a - location_a + i - index_ja = i + if location_a < rows_a + index_ia = rows_a - location_a + i + index_ja = i else - index_ia = i - index_ja = location_a - cols_a + i + index_ia = i + index_ja = location_a - cols_a + i end min_loc = rows_b - index_ja + 1 max_loc = 2 * rows_b - index_ja @@ -309,22 +292,22 @@ function LinearAlgebra.mul!(C::Matrix{T}, A::SparseBandedMatrix{T}, B::SparseBan # if ind < rows(A), then index = (rows - loc + i, i) #rows - loc + j = index_ja, j = index_ja - rows + loc # else index = (i, loc - cols + i) - # if location < rows(B), then + # if location < rows(B), then if location_b <= rows_b && location_b >= min_loc j = index_ja - rows_b + location_b index_jb = j val_b = B.diags[ind_b][j] - C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) + C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) elseif location_b > rows_b && location_b <= max_loc j = index_ja - index_jb = location_b - cols_b + j + index_jb = location_b - cols_b + j val_b = B.diags[ind_b][j] - C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) - end + C[index_ia, index_jb] = muladd(val_a, val_b, C[index_ia, index_jb]) + end end end end - return C + C end export SparseBandedMatrix, size, getindex, setindex!, setdiagonal!, mul! @@ -363,4 +346,4 @@ export SparseBandedMatrix, size, getindex, setindex!, setdiagonal!, mul! end end -end +end \ No newline at end of file diff --git a/test/interface.jl b/test/interface.jl index ce259e0..4056b12 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -1,155 +1,132 @@ -using Test -using SparseBandedMatrices -using LinearAlgebra - -@testset "BigFloat Support" begin - # Test 1: Basic constructor - @testset "Constructor" begin - A = SparseBandedMatrix{BigFloat}(undef, 5, 5) - @test size(A) == (5, 5) - @test eltype(A) == BigFloat - end - - # Test 2: setindex!/getindex - @testset "Indexing" begin - A = SparseBandedMatrix{BigFloat}(undef, 5, 5) - A[1, 1] = BigFloat(2.0) - A[2, 3] = BigFloat(3.14) - @test A[1, 1] == BigFloat(2.0) - @test typeof(A[1, 1]) == BigFloat - @test A[3, 3] == BigFloat(0.0) # Unset element returns zero - @test typeof(A[3, 3]) == BigFloat - end +using SparseBandedMatrices, LinearAlgebra, Test + +@testset "BigFloat support" begin + dim = 5 + diag_locs = [1, 3, 5, 7, 9] + diag_vals = Vector{Vector{BigFloat}}(undef, length(diag_locs)) + for (j, loc) in enumerate(diag_locs) + len = min(loc, 2 * dim - loc) + diag_vals[j] = BigFloat.(rand(len)) + end + + A = SparseBandedMatrix{BigFloat}(diag_locs, diag_vals, dim, dim) + @test eltype(A) == BigFloat + @test size(A) == (dim, dim) + + # Test getindex returns BigFloat + @test A[1, 1] isa BigFloat + @test A[2, 3] isa BigFloat # zero element + + # Test setindex! with BigFloat + A[1, 1] = big"42.0" + @test A[1, 1] == big"42.0" + + # Test mul! with BigFloat matrices + B = Matrix{BigFloat}(ones(BigFloat, dim, 3)) + C = Matrix{BigFloat}(zeros(BigFloat, dim, 3)) + mul!(C, A, B, big"1.0", big"0.0") + @test eltype(C) == BigFloat + + # Verify result + A_dense = zeros(BigFloat, dim, dim) + for i in 1:dim, j in 1:dim + A_dense[i, j] = A[i, j] + end + expected = A_dense * B + @test isapprox(C, expected) + + # Test mul! with Matrix * SparseBandedMatrix + B2 = Matrix{BigFloat}(ones(BigFloat, 3, dim)) + C2 = Matrix{BigFloat}(zeros(BigFloat, 3, dim)) + mul!(C2, B2, A, big"1.0", big"0.0") + expected2 = B2 * A_dense + @test isapprox(C2, expected2) +end - # Test 3: setdiagonal! - @testset "setdiagonal!" begin - A = SparseBandedMatrix{BigFloat}(undef, 5, 5) - diagvals = BigFloat[1.0, 2.0, 3.0] - setdiagonal!(A, diagvals, true) - setdiagonal!(A, BigFloat[4.0, 5.0], false) - # Verify the diagonals were set (checking actual values) - @test A[3, 1] == BigFloat(1.0) # first element of lower diagonal - @test A[4, 2] == BigFloat(2.0) # second element of lower diagonal - @test A[5, 3] == BigFloat(3.0) # third element of lower diagonal - @test A[1, 4] == BigFloat(4.0) # first element of upper diagonal - @test A[2, 5] == BigFloat(5.0) # second element of upper diagonal - end +@testset "ComplexF64 support" begin + A = SparseBandedMatrix{ComplexF64}(undef, 5, 5) + A[1, 1] = 2.0 + 1.0im + A[2, 2] = 3.0 - 2.0im + A[3, 3] = 4.0 + 0.5im - # Test 4: Constructor with diagonal values - @testset "Constructor with diagonals" begin - ind_vals = [2, 7] - diag_vals = [BigFloat[1.0, 2.0], BigFloat[3.0, 4.0, 5.0]] - B = SparseBandedMatrix{BigFloat}(ind_vals, diag_vals, 5, 5) - @test eltype(B) == BigFloat - end + @test eltype(A) == ComplexF64 + @test A[1, 1] == 2.0 + 1.0im + @test A[4, 4] == 0.0 + 0.0im # zero element - # Test 5: mul! with Matrix - @testset "mul! Matrix operations" begin - A = SparseBandedMatrix{BigFloat}(undef, 3, 3) - A[1, 1] = BigFloat(1.0) - A[2, 2] = BigFloat(2.0) - A[3, 3] = BigFloat(3.0) + B = Matrix{ComplexF64}(ones(ComplexF64, 5, 3)) + C = Matrix{ComplexF64}(zeros(ComplexF64, 5, 3)) - B = ones(BigFloat, 3, 2) - C = zeros(BigFloat, 3, 2) + mul!(C, A, B, 1.0 + 0im, 0.0 + 0im) + @test eltype(C) == ComplexF64 - mul!(C, A, B, BigFloat(1.0), BigFloat(0.0)) - expected = Matrix(A) * B - @test isapprox(C, expected) + # Verify result + A_dense = zeros(ComplexF64, 5, 5) + for i in 1:5, j in 1:5 + A_dense[i, j] = A[i, j] end + expected = A_dense * B + @test isapprox(C, expected) - # Test 6: * operator (this was the key bug fix) - @testset "* operator" begin - A = SparseBandedMatrix{BigFloat}(undef, 3, 3) - A[1, 1] = BigFloat(1.0) - A[2, 2] = BigFloat(2.0) - A[3, 3] = BigFloat(3.0) + # Test mul! with Matrix * SparseBandedMatrix + B2 = Matrix{ComplexF64}(ones(ComplexF64, 3, 5)) + C2 = Matrix{ComplexF64}(zeros(ComplexF64, 3, 5)) + mul!(C2, B2, A, 1.0 + 0im, 0.0 + 0im) + expected2 = B2 * A_dense + @test isapprox(C2, expected2) +end - B = ones(BigFloat, 3, 2) - C = A * B - expected = Matrix(A) * B +@testset "Float32 support" begin + A = SparseBandedMatrix{Float32}(undef, 5, 5) + A[1, 1] = 2.0f0 + A[2, 2] = 3.0f0 + A[3, 3] = 4.0f0 - @test eltype(C) == BigFloat - @test isapprox(C, expected) - end + @test eltype(A) == Float32 + @test A[1, 1] == 2.0f0 - # Test 7: Matrix * SparseBandedMatrix - @testset "Matrix * SparseBandedMatrix" begin - A = SparseBandedMatrix{BigFloat}(undef, 3, 3) - A[1, 1] = BigFloat(1.0) - A[2, 2] = BigFloat(2.0) - A[3, 3] = BigFloat(3.0) + B = Matrix{Float32}(ones(Float32, 5, 3)) + C = Matrix{Float32}(zeros(Float32, 5, 3)) - B = ones(BigFloat, 2, 3) - C = B * A - expected = B * Matrix(A) + mul!(C, A, B, 1.0f0, 0.0f0) + @test eltype(C) == Float32 - @test eltype(C) == BigFloat - @test isapprox(C, expected) - end - - # Test 8: SparseBandedMatrix * SparseBandedMatrix - @testset "SparseBandedMatrix * SparseBandedMatrix" begin - A = SparseBandedMatrix{BigFloat}(undef, 3, 3) - A[1, 1] = BigFloat(1.0) - A[2, 2] = BigFloat(2.0) - - B = SparseBandedMatrix{BigFloat}(undef, 3, 3) - B[1, 1] = BigFloat(2.0) - B[3, 3] = BigFloat(3.0) - - C = A * B - expected = Matrix(A) * Matrix(B) - @test isapprox(C, expected) + # Verify result + A_dense = zeros(Float32, 5, 5) + for i in 1:5, j in 1:5 + A_dense[i, j] = A[i, j] end + expected = A_dense * B + @test isapprox(C, expected) end -# Note: ComplexF64 is not fully supported because fma() is not defined for complex numbers. -# This is a known limitation of the current implementation which uses fma for performance. - -@testset "AbstractArray Interface" begin +@testset "AbstractArray interface" begin A = SparseBandedMatrix{Float64}(undef, 5, 5) - A[1, 1] = 1.0 - A[2, 2] = 2.0 - A[3, 3] = 3.0 - A[1, 3] = 0.5 - - @testset "size" begin - @test size(A) == (5, 5) - @test size(A, 1) == 5 - @test size(A, 2) == 5 - end - - @testset "length" begin - @test length(A) == 25 - end + A[1, 1] = 2.0 + A[2, 2] = 3.0 + A[3, 3] = 4.0 - @testset "eltype" begin - @test eltype(A) == Float64 - end + # Test size + @test size(A) == (5, 5) + @test size(A, 1) == 5 + @test size(A, 2) == 5 - @testset "axes" begin - @test axes(A) == (Base.OneTo(5), Base.OneTo(5)) - end + # Test length + @test length(A) == 25 - @testset "firstindex/lastindex" begin - @test firstindex(A) == 1 - @test lastindex(A) == 25 - end + # Test eltype + @test eltype(A) == Float64 - @testset "iteration" begin - vals = collect(A) - @test length(vals) == 25 - @test vals[1] == 1.0 # A[1,1] - end + # Test firstindex/lastindex + @test firstindex(A) == 1 + @test lastindex(A) == 25 - @testset "Matrix conversion" begin - B = Matrix(A) - @test typeof(B) == Matrix{Float64} - @test all(A[i, j] == B[i, j] for i in 1:5, j in 1:5) - end + # Test IndexStyle + @test Base.IndexStyle(typeof(A)) == Base.IndexCartesian() - @testset "IndexStyle" begin - @test Base.IndexStyle(typeof(A)) == IndexCartesian() + # Test iteration + count = 0 + for _ in A + count += 1 end + @test count == 25 end diff --git a/test/runtests.jl b/test/runtests.jl index 709ae97..cb93490 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,95 +1,85 @@ -using SafeTestsets, Test, Pkg - -const GROUP = get(ENV, "GROUP", "All") +using SafeTestsets, Test @testset "SparseBandedMatrices" begin - if GROUP == "All" || GROUP == "Core" - @safetestset "Quality Assurance" include("qa.jl") - @safetestset "Interface Compatibility" include("interface.jl") - - @safetestset "Constructors" begin - using SparseBandedMatrices - - A = SparseBandedMatrix{Float64}(undef, 5, 5) - A[1, 1] = 2 - @test A[1, 1] == 2.0 - A[4, 1] = 0 - @test A[4, 1] == 0.0 - A[1, 3] = 5 - @test A[1, 3] == 5.0 - - @test size(A) == (5, 5) - end - - @safetestset "Multiplication" begin - using SparseBandedMatrices, Random - dim = 5000 - x = rand(10:75) - diag_vals = Vector{Vector{Float64}}(undef, x) - diag_locs = randperm(dim * 2 - 1)[1:x] - for j in 1:x - diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) - end - - x_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) - x_dense = copy(x_butterfly) + @safetestset "Quality Assurance" include("qa.jl") + @safetestset "JET Static Analysis" include("jet.jl") + @safetestset "Interface Compatibility" include("interface.jl") + + @safetestset "Constructors" begin + using SparseBandedMatrices + + A = SparseBandedMatrix{Float64}(undef, 5, 5) + A[1, 1] = 2 + @test A[1, 1] == 2.0 + A[4, 1] = 0 + @test A[4, 1] == 0.0 + A[1, 3] = 5 + @test A[1, 3] == 5.0 + + @test size(A) == (5,5) + end - y = rand(dim, dim) - z = zeros(dim, dim) + @safetestset "Multiplication" begin + using SparseBandedMatrices, Random + dim = 5000 + x = rand(10:75) + diag_vals = Vector{Vector{Float64}}(undef, x) + diag_locs = randperm(dim * 2 - 1)[1:x] + for j in 1:x + diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) + end - @test isapprox(x_dense * y, x_butterfly * y) - @test isapprox(y * x_dense, y * x_butterfly) + x_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) + x_dense = copy(x_butterfly) - y = rand(10:75) - diag_vals = Vector{Vector{Float64}}(undef, y) - diag_locs = randperm(dim * 2 - 1)[1:y] - for j in 1:y - diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) - end + y = rand(dim, dim) + z = zeros(dim, dim) - y_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) - y_dense = copy(y_butterfly) + @test isapprox(x_dense * y, x_butterfly * y) + @test isapprox(y * x_dense, y * x_butterfly) - @test isapprox(x_butterfly * y_butterfly, x_dense * y_dense) + y = rand(10:75) + diag_vals = Vector{Vector{Float64}}(undef, y) + diag_locs = randperm(dim * 2 - 1)[1:y] + for j in 1:y + diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) end - @safetestset "Division" begin - using SparseBandedMatrices, Random - dim = 5000 - x = rand(10:75) - diag_vals = Vector{Vector{Float64}}(undef, x) - diag_locs = randperm(dim * 2 - 1)[1:x] - for j in 1:x - diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) - end + y_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) + y_dense = copy(y_butterfly) - x_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) - x_dense = copy(x_butterfly) + @test isapprox(x_butterfly * y_butterfly, x_dense * y_dense) + end - y = rand(dim, dim) - z = zeros(dim, dim) + @safetestset "Division" begin + using SparseBandedMatrices, Random + dim = 5000 + x = rand(10:75) + diag_vals = Vector{Vector{Float64}}(undef, x) + diag_locs = randperm(dim * 2 - 1)[1:x] + for j in 1:x + diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) + end - @test isapprox(x_dense / y, x_butterfly / y) - @test isapprox(y / x_dense, y / x_butterfly) + x_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) + x_dense = copy(x_butterfly) - y = rand(10:75) - diag_vals = Vector{Vector{Float64}}(undef, y) - diag_locs = randperm(dim * 2 - 1)[1:y] - for j in 1:y - diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) - end + y = rand(dim, dim) + z = zeros(dim, dim) - y_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) - y_dense = copy(y_butterfly) + @test isapprox(x_dense / y, x_butterfly / y) + @test isapprox(y / x_dense, y / x_butterfly) - @test isapprox(x_butterfly / y_butterfly, x_dense / y_dense) + y = rand(10:75) + diag_vals = Vector{Vector{Float64}}(undef, y) + diag_locs = randperm(dim * 2 - 1)[1:y] + for j in 1:y + diag_vals[j] = rand(min(diag_locs[j], 2 * dim - diag_locs[j])) end - end - if GROUP == "All" || GROUP == "nopre" - Pkg.activate(joinpath(@__DIR__, "nopre")) - Pkg.develop(PackageSpec(path = joinpath(@__DIR__, ".."))) - Pkg.instantiate() - @safetestset "JET Static Analysis" include("nopre/jet.jl") + y_butterfly = SparseBandedMatrix{Float64}(diag_locs, diag_vals, dim, dim) + y_dense = copy(y_butterfly) + + @test isapprox(x_butterfly / y_butterfly, x_dense / y_dense) end -end +end \ No newline at end of file