Fix type conversion for weird orders (#226)

* Fix type conversion for weird orders

* Fix

* Add tests

Author: gdalle

Project.toml

@@ -1,7 +1,7 @@
 name = "SparseMatrixColorings"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
 authors = ["Guillaume Dalle", "Alexis Montoison"]
-version = "0.4.15"
+version = "0.4.16"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

src/order.jl

@@ -137,16 +137,17 @@ struct DegreeBuckets{T}
     positions::Vector{T}
 end
 
-function DegreeBuckets(::Type{T}, degrees::Vector{<:Integer}, dmax::Integer) where {T}
+function DegreeBuckets(::Type{T}, degrees::Vector{T}, dmax::Integer) where {T}
     # number of vertices per degree class
     deg_count = zeros(T, dmax + 1)
     for d in degrees
         deg_count[d + 1] += 1
     end
     # bucket limits
-    bucket_high = cumsum(deg_count)
-    bucket_low = vcat(zero(T), @view(bucket_high[1:(end - 1)]))
-    bucket_low .+= 1
+    bucket_high = accumulate(+, deg_count)
+    bucket_low = similar(bucket_high)
+    bucket_low[1] = 1
+    bucket_low[2:end] .= @view(bucket_high[1:(end - 1)]) .+ 1
     # assign each vertex to the correct position inside its degree class
     bucket_storage = similar(degrees, T)
     positions = similar(degrees, T)
@@ -253,7 +254,7 @@ function vertices(
     if degree_increasing(; degtype, direction)
         degrees = zeros(T, nb_vertices(g))
     else
-        degrees = [degree(g, v) for v in vertices(g)]
+        degrees = T[degree(g, v) for v in vertices(g)]
     end
     db = DegreeBuckets(T, degrees, maximum_degree(g))
     π = T[]
@@ -378,3 +379,14 @@ The `elimination_algorithm` must be an instance of `CliqueTrees.EliminationAlgor
 struct PerfectEliminationOrder{E} <: AbstractOrder
     elimination_algorithm::E
 end
+
+function all_orders()
+    return [
+        NaturalOrder(),
+        RandomOrder(),
+        LargestFirst(),
+        SmallestLast(),
+        IncidenceDegree(),
+        DynamicLargestFirst(),
+    ]
+end

test/allocations.jl

@@ -132,3 +132,33 @@ end
         test_noallocs_structured_decompression(1000; structure, partition, decompression)
     end
 end
+
+@testset "Single precision" begin
+    n, p = 10_000, 0.001
+    A64 = sparse(Symmetric(sprand(rng, Float32, 100, 100, 0.1)))
+    A32 = convert(SparseMatrixCSC{Float32,Int32}, A64)
+    @testset "$structure - $partition - $decompression" for (
+        structure, partition, decompression
+    ) in [
+        (:nonsymmetric, :column, :direct),
+        (:nonsymmetric, :row, :direct),
+        (:symmetric, :column, :direct),
+        (:symmetric, :column, :substitution),
+        (:nonsymmetric, :bidirectional, :direct),
+        (:nonsymmetric, :bidirectional, :substitution),
+    ]
+        @testset for order in
+                     (NaturalOrder(), RandomOrder(), LargestFirst(), DynamicLargestFirst())
+            problem = ColoringProblem(; structure, partition)
+            algo = GreedyColoringAlgorithm(order; decompression)
+            b64 = @b fast_coloring(A64, problem, algo)
+            b32 = @b fast_coloring(A32, problem, algo)
+            # check that we allocate no more than 50% + epsilon with Int32
+            if decompression == :direct
+                @test b32.bytes < 0.6 * b64.bytes
+            else
+                @test_broken b32.bytes < 0.6 * b64.bytes
+            end
+        end
+    end
+end;

test/type_stability.jl

@@ -3,7 +3,7 @@ using JET
 using LinearAlgebra
 using SparseArrays
 using SparseMatrixColorings
-using SparseMatrixColorings: matrix_versions, respectful_similar
+using SparseMatrixColorings: matrix_versions, respectful_similar, all_orders
 using StableRNGs
 using Test
 
@@ -178,25 +178,27 @@ end;
         (:nonsymmetric, :bidirectional, :direct),
         (:nonsymmetric, :bidirectional, :substitution),
     ]
-        result = coloring(
-            A,
-            ColoringProblem(; structure, partition),
-            GreedyColoringAlgorithm(; decompression);
-        )
-        if partition in (:column, :bidirectional)
-            @test eltype(column_colors(result)) == Int32
-            @test eltype(column_groups(result)[1]) == Int32
-        end
-        if partition in (:row, :bidirectional)
-            @test eltype(row_colors(result)) == Int32
-            @test eltype(row_groups(result)[1]) == Int32
-        end
-        if partition == :bidirectional
-            Br, Bc = compress(A, result)
-            @test decompress(Br, Bc, result) isa SparseMatrixCSC{Float32,Int32}
-        else
-            B = compress(A, result)
-            @test decompress(B, result) isa SparseMatrixCSC{Float32,Int32}
+        @testset for order in all_orders()
+            result = coloring(
+                A,
+                ColoringProblem(; structure, partition),
+                GreedyColoringAlgorithm(order; decompression);
+            )
+            if partition in (:column, :bidirectional)
+                @test eltype(column_colors(result)) == Int32
+                @test eltype(column_groups(result)[1]) == Int32
+            end
+            if partition in (:row, :bidirectional)
+                @test eltype(row_colors(result)) == Int32
+                @test eltype(row_groups(result)[1]) == Int32
+            end
+            if partition == :bidirectional
+                Br, Bc = compress(A, result)
+                @test decompress(Br, Bc, result) isa SparseMatrixCSC{Float32,Int32}
+            else
+                B = compress(A, result)
+                @test decompress(B, result) isa SparseMatrixCSC{Float32,Int32}
+            end
         end
     end
 end