add tests for bareiss, fix precompile (#1455)

Co-authored-by: Christopher Rackauckas <[email protected]>

Author: oscardssmith

src/ModelingToolkit.jl

@@ -107,6 +107,8 @@ Get the set of parameters variables for the given system.
 """
 function parameters end
 
+# this has to be included early to deal with depency issues
+include("structural_transformation/bareiss.jl")
 include("bipartite_graph.jl")
 using .BipartiteGraphs
 
@@ -116,9 +118,6 @@ include("parameters.jl")
 include("utils.jl")
 include("domains.jl")
 
-# Code that should eventually go elsewhere, but is here for fow
-include("structural_transformation/bareiss.jl")
-
 include("systems/abstractsystem.jl")
 include("systems/connectors.jl")
 

src/structural_transformation/bareiss.jl

@@ -2,111 +2,116 @@
 
 using LinearAlgebra
 using SparseArrays
-using SparseArrays: AbstractSparseMatrixCSC
+using SparseArrays: AbstractSparseMatrixCSC, getcolptr
 
 macro swap(a, b)
     esc(:(($a, $b) = ($b, $a)))
 end
 
-function swaprows!(a::AbstractMatrix, i, j)
-    i == j && return
-    rows = axes(a,1)
-    @boundscheck i in rows || throw(BoundsError(a, (:,i)))
-    @boundscheck j in rows || throw(BoundsError(a, (:,j)))
-    for k in axes(a,2)
-        @inbounds a[i,k],a[j,k] = a[j,k],a[i,k]
+# https://github.com/JuliaLang/julia/pull/42678
+@static if VERSION > v"1.8.0-DEV.762"
+    import Base: swaprows!
+else
+    function swaprows!(a::AbstractMatrix, i, j)
+        i == j && return
+        rows = axes(a,1)
+        @boundscheck i in rows || throw(BoundsError(a, (:,i)))
+        @boundscheck j in rows || throw(BoundsError(a, (:,j)))
+        for k in axes(a,2)
+            @inbounds a[i,k],a[j,k] = a[j,k],a[i,k]
+        end
     end
-end
-function Base.circshift!(a::AbstractVector, shift::Integer)
-    n = length(a)
-    n == 0 && return
-    shift = mod(shift, n)
-    shift == 0 && return
-    reverse!(a, 1, shift)
-    reverse!(a, shift+1, length(a))
-    reverse!(a)
-    return a
-end
-function Base.swapcols!(A::AbstractSparseMatrixCSC, i, j)
-    i == j && return
-
-    # For simplicitly, let i denote the smaller of the two columns
-    j < i && @swap(i, j)
-
-    colptr = getcolptr(A)
-    irow = colptr[i]:(colptr[i+1]-1)
-    jrow = colptr[j]:(colptr[j+1]-1)
-
-    function rangeexchange!(arr, irow, jrow)
-        if length(irow) == length(jrow)
-            for (a, b) in zip(irow, jrow)
-                @inbounds @swap(arr[i], arr[j])
+    function Base.circshift!(a::AbstractVector, shift::Integer)
+        n = length(a)
+        n == 0 && return
+        shift = mod(shift, n)
+        shift == 0 && return
+        reverse!(a, 1, shift)
+        reverse!(a, shift+1, length(a))
+        reverse!(a)
+        return a
+    end
+    function Base.swapcols!(A::AbstractSparseMatrixCSC, i, j)
+        i == j && return
+
+        # For simplicitly, let i denote the smaller of the two columns
+        j < i && @swap(i, j)
+
+        colptr = getcolptr(A)
+        irow = colptr[i]:(colptr[i+1]-1)
+        jrow = colptr[j]:(colptr[j+1]-1)
+
+        function rangeexchange!(arr, irow, jrow)
+            if length(irow) == length(jrow)
+                for (a, b) in zip(irow, jrow)
+                    @inbounds @swap(arr[i], arr[j])
+                end
+                return
             end
-            return
+            # This is similar to the triple-reverse tricks for
+            # circshift!, except that we have three ranges here,
+            # so it ends up being 4 reverse calls (but still
+            # 2 overall reversals for the memory range). Like
+            # circshift!, there's also a cycle chasing algorithm
+            # with optimal memory complexity, but the performance
+            # tradeoffs against this implementation are non-trivial,
+            # so let's just do this simple thing for now.
+            # See https://github.com/JuliaLang/julia/pull/42676 for
+            # discussion of circshift!-like algorithms.
+            reverse!(@view arr[irow])
+            reverse!(@view arr[jrow])
+            reverse!(@view arr[(last(irow)+1):(first(jrow)-1)])
+            reverse!(@view arr[first(irow):last(jrow)])
         end
-        # This is similar to the triple-reverse tricks for
-        # circshift!, except that we have three ranges here,
-        # so it ends up being 4 reverse calls (but still
-        # 2 overall reversals for the memory range). Like
-        # circshift!, there's also a cycle chasing algorithm
-        # with optimal memory complexity, but the performance
-        # tradeoffs against this implementation are non-trivial,
-        # so let's just do this simple thing for now.
-        # See https://github.com/JuliaLang/julia/pull/42676 for
-        # discussion of circshift!-like algorithms.
-        reverse!(@view arr[irow])
-        reverse!(@view arr[jrow])
-        reverse!(@view arr[(last(irow)+1):(first(jrow)-1)])
-        reverse!(@view arr[first(irow):last(jrow)])
-    end
-    rangeexchange!(rowvals(A), irow, jrow)
-    rangeexchange!(nonzeros(A), irow, jrow)
+        rangeexchange!(rowvals(A), irow, jrow)
+        rangeexchange!(nonzeros(A), irow, jrow)
 
-    if length(irow) != length(jrow)
-        @inbounds colptr[i+1:j] .+= length(jrow) - length(irow)
+        if length(irow) != length(jrow)
+            @inbounds colptr[i+1:j] .+= length(jrow) - length(irow)
+        end
+        return nothing
     end
-    return nothing
-end
-function swaprows!(A::AbstractSparseMatrixCSC, i, j)
-    # For simplicitly, let i denote the smaller of the two rows
-    j < i && @swap(i, j)
-
-    rows = rowvals(A)
-    vals = nonzeros(A)
-    for col = 1:size(A, 2)
-        rr = nzrange(A, col)
-        iidx = searchsortedfirst(@view(rows[rr]), i)
-        has_i = iidx <= length(rr) && rows[rr[iidx]] == i
-
-        jrange = has_i ? (iidx:last(rr)) : rr
-        jidx = searchsortedlast(@view(rows[jrange]), j)
-        has_j = jidx != 0 && rows[jrange[jidx]] == j
-
-        if !has_j && !has_i
-            # Has neither row - nothing to do
-            continue
-        elseif has_i && has_j
-            # This column had both i and j rows - swap them
-            @swap(vals[rr[iidx]], vals[jrange[jidx]])
-        elseif has_i
-            # Update the rowval and then rotate both nonzeros
-            # and the remaining rowvals into the correct place
-            rows[rr[iidx]] = j
-            jidx == 0 && continue
-            rotate_range = rr[iidx]:jrange[jidx]
-            circshift!(@view(vals[rotate_range]), -1)
-            circshift!(@view(rows[rotate_range]), -1)
-        else
-            # Same as i, but in the opposite direction
-            @assert has_j
-            rows[jrange[jidx]] = i
-            iidx > length(rr) && continue
-            rotate_range = rr[iidx]:jrange[jidx]
-            circshift!(@view(vals[rotate_range]), 1)
-            circshift!(@view(rows[rotate_range]), 1)
+    function swaprows!(A::AbstractSparseMatrixCSC, i, j)
+        # For simplicitly, let i denote the smaller of the two rows
+        j < i && @swap(i, j)
+
+        rows = rowvals(A)
+        vals = nonzeros(A)
+        for col = 1:size(A, 2)
+            rr = nzrange(A, col)
+            iidx = searchsortedfirst(@view(rows[rr]), i)
+            has_i = iidx <= length(rr) && rows[rr[iidx]] == i
+
+            jrange = has_i ? (iidx:last(rr)) : rr
+            jidx = searchsortedlast(@view(rows[jrange]), j)
+            has_j = jidx != 0 && rows[jrange[jidx]] == j
+
+            if !has_j && !has_i
+                # Has neither row - nothing to do
+                continue
+            elseif has_i && has_j
+                # This column had both i and j rows - swap them
+                @swap(vals[rr[iidx]], vals[jrange[jidx]])
+            elseif has_i
+                # Update the rowval and then rotate both nonzeros
+                # and the remaining rowvals into the correct place
+                rows[rr[iidx]] = j
+                jidx == 0 && continue
+                rotate_range = rr[iidx]:jrange[jidx]
+                circshift!(@view(vals[rotate_range]), -1)
+                circshift!(@view(rows[rotate_range]), -1)
+            else
+                # Same as i, but in the opposite direction
+                @assert has_j
+                rows[jrange[jidx]] = i
+                iidx > length(rr) && continue
+                rotate_range = rr[iidx]:jrange[jidx]
+                circshift!(@view(vals[rotate_range]), 1)
+                circshift!(@view(rows[rotate_range]), 1)
+            end
         end
+        return nothing
     end
-    return nothing
 end
 
 function bareiss_update!(zero!, M::StridedMatrix, k, swapto, pivot, prev_pivot)

test/structural_transformation/bareiss.jl

@@ -0,0 +1,28 @@
+using SparseArrays
+using ModelingToolkit
+import ModelingToolkit: bareiss!, find_pivot_col, bareiss_update!, swaprows!
+import Base: swapcols!
+
+function det_bareiss!(M)
+    parity = 1
+    _swaprows!(M, i, j) = (i != j && (parity = -parity); swaprows!(M, i, j))
+    _swapcols!(M, i, j) = (i != j && (parity = -parity); swapcols!(M, i, j))
+    # We only look at the last entry, so we don't care that the sub-diagonals are
+    # garbage.
+    zero!(M, i, j) = nothing
+    rank = bareiss!(M, (_swapcols!, _swaprows!, bareiss_update!, zero!);
+        find_pivot=find_pivot_col)
+    return parity * M[end,end]
+end
+
+@testset "bareiss tests" begin
+    # copy gives a dense matrix
+    @testset "bareiss tests: $T" for T in (copy, sparse)
+        # matrix determinent pairs
+        for (M, d) in ((BigInt[9 1 8 0; 0 0 8 7; 7 6 8 3; 2 9 7 7], -1),
+                   (BigInt[1 big(2)^65+1; 3 4], 4-3*(big(2)^65+1)))
+            # test that the determinent was correctly computed
+            @test det_bareiss!(T(M)) == d
+        end
+    end
+end

test/structural_transformation/runtests.jl

@@ -3,3 +3,4 @@ using SafeTestsets
 @safetestset "Utilities" begin include("utils.jl") end
 @safetestset "Index Reduction & SCC" begin include("index_reduction.jl") end
 @safetestset "Tearing" begin include("tearing.jl") end
+@safetestset "Bareiss" begin include("bareiss.jl") end