Fix mul! for uninitialized matrices and add interface tests

.github/workflows/Tests.yml

@@ -22,12 +22,20 @@ jobs:
   tests:
     name: "Tests"
     strategy:
+      fail-fast: false
       matrix:
         version:
           - "1"
           - "lts"
           - "pre"
+        group:
+          - Core
+          - nopre
+        exclude:
+          - version: "pre"
+            group: nopre
     uses: "SciML/.github/.github/workflows/tests.yml@v1"
     with:
       julia-version: "${{ matrix.version }}"
-    secrets: "inherit"
+      group: "${{ matrix.group }}"
+    secrets: "inherit"

docs/Project.toml

@@ -1,6 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-SparseBandedMatrices = "b42f87ec-cfbe-4098-ba17-1affa37bb2d8"
+SparseBandedMatrices = "bd59d7e1-4699-4102-944e-d05209cb92aa"
 
 [compat]
 Documenter = "1"

docs/make.jl

@@ -1,8 +1,13 @@
 using SparseBandedMatrices
 using Documenter
 
-cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
-cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
+# Copy manifest files if they exist (they are generated by CI)
+# Create assets directory if needed
+if isfile("./docs/Manifest.toml") || isfile("./docs/Project.toml")
+    mkpath("./docs/src/assets")
+end
+isfile("./docs/Manifest.toml") && cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
+isfile("./docs/Project.toml") && cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
 
 makedocs(;
     modules = [SparseBandedMatrices],

src/SparseBandedMatrices.jl

@@ -169,7 +169,11 @@ function LinearAlgebra.mul!(C::Matrix{T}, A::SparseBandedMatrix{T}, B::Matrix{T}
     @assert size(A, 2) == size(B, 1)
     @assert size(A, 1) == size(C, 1)
     @assert size(B, 2) == size(C, 2)
-    C .*= b
+    if iszero(b)
+        fill!(C, zero(T))
+    else
+        C .*= b
+    end
 
     rows, cols = size(A)
     @inbounds for (ind, location) in enumerate(A.indices)
@@ -195,13 +199,17 @@ function LinearAlgebra.mul!(C::Matrix{T}, A::SparseBandedMatrix{T}, B::Matrix{T}
     return C
 end
 
-# C = Cb + aBA
+# C = C*b + a*B*A
 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)
 
-    C .*= b
+    if iszero(b)
+        fill!(C, zero(T))
+    else
+        C .*= b
+    end
 
     rows, cols = size(B)
     @inbounds for (ind, location) in enumerate(B.indices)
@@ -227,7 +235,12 @@ function LinearAlgebra.mul!(C::SparseBandedMatrix{T}, A::SparseBandedMatrix{T},
     @assert size(A, 1) == size(C, 1)
     @assert size(B, 2) == size(C, 2)
 
-    C .*= b
+    if iszero(b)
+        empty!(C.indices)
+        empty!(C.diags)
+    else
+        C .*= b
+    end
 
     rows_a, cols_a = size(A)
     rows_b, cols_b = size(B)
@@ -271,7 +284,11 @@ function LinearAlgebra.mul!(C::Matrix{T}, A::SparseBandedMatrix{T}, B::SparseBan
     @assert size(A, 1) == size(C, 1)
     @assert size(B, 2) == size(C, 2)
 
-    C .*= b
+    if iszero(b)
+        fill!(C, zero(T))
+    else
+        C .*= b
+    end
 
     rows_a, cols_a = size(A)
     rows_b, cols_b = size(B)

test/Project.toml

@@ -1,7 +1,16 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+SparseBandedMatrices = "bd59d7e1-4699-4102-944e-d05209cb92aa"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+Aqua = "0.8"
+LinearAlgebra = "1"
+Pkg = "1"
+Random = "1"
+SafeTestsets = "0.1"
+Test = "1"

test/interface.jl

@@ -0,0 +1,155 @@
+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
+
+    # 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
+
+    # 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 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 = ones(BigFloat, 3, 2)
+        C = zeros(BigFloat, 3, 2)
+
+        mul!(C, A, B, BigFloat(1.0), BigFloat(0.0))
+        expected = Matrix(A) * B
+        @test isapprox(C, expected)
+    end
+
+    # 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)
+
+        B = ones(BigFloat, 3, 2)
+        C = A * B
+        expected = Matrix(A) * B
+
+        @test eltype(C) == BigFloat
+        @test isapprox(C, expected)
+    end
+
+    # 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 = ones(BigFloat, 2, 3)
+        C = B * A
+        expected = B * Matrix(A)
+
+        @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)
+    end
+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
+    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
+
+    @testset "eltype" begin
+        @test eltype(A) == Float64
+    end
+
+    @testset "axes" begin
+        @test axes(A) == (Base.OneTo(5), Base.OneTo(5))
+    end
+
+    @testset "firstindex/lastindex" begin
+        @test firstindex(A) == 1
+        @test lastindex(A) == 25
+    end
+
+    @testset "iteration" begin
+        vals = collect(A)
+        @test length(vals) == 25
+        @test vals[1] == 1.0  # A[1,1]
+    end
+
+    @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
+
+    @testset "IndexStyle" begin
+        @test Base.IndexStyle(typeof(A)) == IndexCartesian()
+    end
+end

test/nopre/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseBandedMatrices = "bd59d7e1-4699-4102-944e-d05209cb92aa"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+JET = "0.9, 0.10, 0.11"
+LinearAlgebra = "1"
+Test = "1"

test/nopre/runtests.jl

@@ -0,0 +1,5 @@
+using Test
+
+@testset "nopre tests" begin
+    include("jet.jl")
+end