Merge pull request #1554 from SciML/myb/optimize

Optimization

Author: YingboMa

src/structural_transformation/bareiss.jl

@@ -114,20 +114,40 @@ else
     end
 end
 
+function bareiss_update!(zero!, M::StridedMatrix, k, swapto, pivot, prev_pivot::Base.BitInteger)
+    flag = zero(prev_pivot)
+    prev_pivot = Base.MultiplicativeInverses.SignedMultiplicativeInverse(prev_pivot)
+    @inbounds for i in k+1:size(M, 2)
+        Mki = M[k,i]
+        @simd ivdep for j in k+1:size(M, 1)
+            M[j,i], r = divrem(M[j,i]*pivot - M[j,k]*Mki, prev_pivot)
+            flag = flag | r
+        end
+    end
+    iszero(flag) || error("Overflow occurred")
+    zero!(M, k+1:size(M, 1), k)
+end
+
 function bareiss_update!(zero!, M::StridedMatrix, k, swapto, pivot, prev_pivot)
-    for i in k+1:size(M, 2), j in k+1:size(M, 1)
+    @inbounds for i in k+1:size(M, 2), j in k+1:size(M, 1)
         M[j,i] = exactdiv(M[j,i]*pivot - M[j,k]*M[k,i], prev_pivot)
     end
     zero!(M, k+1:size(M, 1), k)
 end
 
 @views function bareiss_update!(zero!, M::AbstractMatrix, k, swapto, pivot, prev_pivot)
+    if prev_pivot isa Base.BitInteger
+        prev_pivot = Base.MultiplicativeInverses.SignedMultiplicativeInverse(prev_pivot)
+    end
     V = M[k+1:end, k+1:end]
     V .= exactdiv.(V .* pivot .- M[k+1:end, k] * M[k, k+1:end]', prev_pivot)
     zero!(M, k+1:size(M, 1), k)
 end
 
 function bareiss_update_virtual_colswap!(zero!, M::AbstractMatrix, k, swapto, pivot, prev_pivot)
+    if prev_pivot isa Base.BitInteger
+        prev_pivot = Base.MultiplicativeInverses.SignedMultiplicativeInverse(prev_pivot)
+    end
     V = @view M[k+1:end, :]
     V .= @views exactdiv.(V .* pivot .- M[k+1:end, swapto[2]] * M[k, :]', prev_pivot)
     zero!(M, k+1:size(M, 1), swapto[2])

src/systems/alias_elimination.jl

@@ -324,13 +324,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     return ag, mm
 end
 
-iszeroterm(v_types, v) = v_types[v] == 0
-isirreducible(v_types, v) = v_types[v] == KEEP
-isalias(v_types, v) = v_types[v] > 0 && !isirreducible(v_types, v)
-alias(v_types, v) = v_types[v]
-negalias(v_types, v) = -v_types[v]
-
-function exactdiv(a::Integer, b::Integer)
+function exactdiv(a::Integer, b)
     d, r = divrem(a, b)
     @assert r == 0
     return d