Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
fix
Browse files Browse the repository at this point in the history
kalmarek committed Jul 8, 2024
1 parent ecdff74 commit 4f1f281
Showing 2 changed files with 104 additions and 11 deletions.
24 changes: 13 additions & 11 deletions .github/workflows/examples.yml
Original file line number Diff line number Diff line change
@@ -9,9 +9,7 @@ jobs:
strategy:
matrix:
version:
- '1.6'
- '1'
- 'nightly'
os:
- ubuntu-latest
arch:
@@ -35,7 +33,7 @@ jobs:
${{ runner.os }}-test-${{ env.cache-name }}-
${{ runner.os }}-test-
${{ runner.os }}-
- name: Example project instantiation
- name: Examples project instantiation
shell: julia --project=examples {0}
run: |
using Pkg
@@ -46,21 +44,25 @@ jobs:
- name: Example C₂ linear action
shell: julia --project=examples {0}
run: |
include("examples/preamble.jl");
include("examples/ex_C2.linear.jl")
wd = walkdir(pwd())
@info first(wd)
@info isdir("examples")
@info isfile("examples/preamble.jl")
include(joinpath(pwd(), "examples/preamble.jl"))
include(joinpath(pwd(), "examples/ex_C2.linear.jl"))
- name: Example a polynomial with Sym(4)-symmetry
shell: julia --project=examples {0}
run: |
include("examples/preamble.jl");
include("examples/ex_S4.jl")
include("./examples/preamble.jl");
include("./examples/ex_S4.jl")
- name: Example Motzkin polynomial
shell: julia --project=examples {0}
run: |
include("examples/preamble.jl");
include("examples/ex_motzkin.jl")
include("./examples/preamble.jl");
include("./examples/ex_motzkin.jl")
- name: Example Robinson Form
shell: julia --project=examples {0}
run: |
include("examples/preamble.jl");
include("examples/ex_robinson_form.jl")
include("./examples/preamble.jl");
include("./examples/ex_robinson_form.jl")
91 changes: 91 additions & 0 deletions bareiss.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
# Bareiss algorithm
function bareiss_update!(zero!, M::Matrix, k, swapto, pivot, prev_pivot)
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

function bareiss_update!(zero!, M::AbstractMatrix, k, swapto, pivot, prev_pivot)
V = @view 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)
V = @view M[k+1:end, :]
V .= exactdiv.(V * pivot - M[k+1:end, swapto[2]] * M[k, :]', prev_pivot)
zero!(M, k+1:size(M, 1), swapto[2])
end

bareiss_zero!(M, i, j) = M[i,j] .= zero(eltype(M))

function find_pivot_col(M, i)
p = findfirst(!iszero, @view M[i,i:end])
p === nothing && return nothing
idx = CartesianIndex(i, p + i - 1)
(idx, M[idx])
end

function find_pivot_any(M, i)
p = findfirst(!iszero, @view M[i:end,i:end])
p === nothing && return nothing
idx = p + CartesianIndex(i - 1, i - 1)
(idx, M[idx])
end

const bareiss_colswap = (Base.swapcols!, Base.swaprows!, bareiss_update!, bareiss_zero!)
const bareiss_virtcolswap = ((M,i,j)->nothing, Base.swaprows!, bareiss_update_virtual_colswap!, bareiss_zero!)

"""
bareiss!(M)
Perform Bareiss's fraction-free row-reduction algorithm on the matrix `M`.
Optionally, a specific pivoting method may be specified.
"""
function bareiss!(M::AbstractMatrix,
(swapcols!, swaprows!, update!, zero!) = bareiss_colswap;
find_pivot=find_pivot_any)
prev = one(eltype(M))
n = size(M, 1)
for k in 1:n
r = find_pivot(M, k)
r === nothing && return k - 1
(swapto, pivot) = r
if CartesianIndex(k, k) != swapto
swapcols!(M, k, swapto[2])
swaprows!(M, k, swapto[1])
end
update!(zero!, M, k, swapto, pivot, prev)
prev = pivot
end
return n
end

"""
det_bareiss!(M)
Calculates the determinant of a matrix using the
[Bareiss Algorithm](https://en.wikipedia.org/wiki/Bareiss_algorithm) using
inplace operations.
# Examples
```jldoctest
julia> M = [1 0; 2 2]
2×2 Matrix{Int64}:
1 0
2 2
julia> LinearAlgebra.det_bareiss!(M)
2
```
"""
function det_bareiss!(M)
n = checksquare(M)
parity = true
swaprows!(M, i, j) = (i != j && (parity = !parity); Base.swaprows!(M, i, j))
swapcols!(M, i, j) = (i != j && (parity = !parity); Base.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)
rank != n && return zero(eltype(M))
return parity ? M[n,n] : -M[n, n]
end

0 comments on commit 4f1f281

Please sign in to comment.