Skip to content

Commit

Permalink
update eco10 example
Browse files Browse the repository at this point in the history
  • Loading branch information
@sumiya11
sumiya11 committed Jan 29, 2024
1 parent de6aa66 commit 0ec642c
Show file tree
Hide file tree
Showing 7 changed files with 360 additions and 164 deletions.
145 changes: 6 additions & 139 deletions src/f4/linalg/backend_threaded.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,7 @@ function linalg_deterministic_sparse_threaded!(
@log level = -3 "linalg_deterministic_sparse!"
@log level = -3 matrix_string_repr(matrix)
# Reduce CD with AB
if true
linalg_reduce_matrix_lower_part_threaded_cas!(matrix, basis, arithmetic)
elseif false
linalg_reduce_matrix_lower_part_threaded_task_lock_free!(matrix, basis, arithmetic)
else
linalg_reduce_matrix_lower_part_threaded_task_cas!(matrix, basis, arithmetic)
end
linalg_reduce_matrix_lower_part_threaded_cas!(matrix, basis, arithmetic)
# Interreduce CD
linalg_interreduce_matrix_pivots!(matrix, basis, arithmetic)
true
Expand All @@ -42,8 +36,6 @@ end
# A B
# 0 D'
# The new pivots in the D' block are known, but possibly not fully interreduced.
#
# NOTE: This function is incorrect.
function linalg_reduce_matrix_lower_part_threaded_cas!(
matrix::MacaulayMatrix{CoeffType},
basis::Basis{CoeffType},
Expand All @@ -62,137 +54,12 @@ function linalg_reduce_matrix_lower_part_threaded_cas!(
resize!(matrix.some_coeffs, nlow)
resize!(matrix.sentinels, ncols)

# If there is a pivot with the leading column i, then sentinels[i] = 1.
# Otherwise, sentinels[i] = 0.
# Once sentinels[i] becomes 1, this equality is unchanged.
# NOTE: for all valid indices, isassigned(sentinels, i) must hold.
sentinels = matrix.sentinels
@inbounds for i in 1:ncols
sentinels[i] = 0
end
@inbounds for i in 1:nup
sentinels[matrix.upper_rows[i][1]] = 1
end

# Allocate thread-local buffers
buffers_row = map(_ -> zeros(AccumType, ncols), 1:nthreads())

# NOTE: by default, @threads uses the :dynamic execution schedule, which
# does not guarantee that threadid() is constant within one iteration
@inbounds Base.Threads.@threads for i in 1:nlow
t_id = threadid()
t_local_row = buffers_row[t_id]
new_sparse_row_support, new_sparse_row_coeffs =
linalg_new_empty_sparse_row(CoeffType)

# Select a row to be reduced from the lower part of the matrix
# NOTE: no copy of coefficients is needed
sparse_row_support = matrix.lower_rows[i]
sparse_row_coeffs = basis.coeffs[row_index_to_coeffs[i]]
@invariant length(sparse_row_support) == length(sparse_row_coeffs)

# Load the coefficients into a dense array
linalg_load_sparse_row!(t_local_row, sparse_row_support, sparse_row_coeffs)

first_nnz_column = sparse_row_support[1]

# In a somewhat ideal scenario, the while loop below executes just once
success = false
@inbounds while !success
@invariant 1 <= first_nnz_column <= length(t_local_row)

# Note that this call may only mutate the first three arguments, and
# the arguments do not overlap in memory with each other
zeroed = linalg_reduce_dense_row_by_pivots_sparse!(
new_sparse_row_support,
new_sparse_row_coeffs,
t_local_row,
matrix,
basis,
pivots,
first_nnz_column,
ncols,
arithmetic,
tmp_pos=-1
)

# If the row is fully reduced
zeroed && break

# At this point, we have discovered a row that is not reduced to
# zero (yet). Thus, two outcomes are possible:
# 1. sentinels[...] is 0, meaning that this row is potentially a new
# pivot row. We try to set sentinels[...] to 1, update the pivots
# and the matrix, and break out of the loop.
# If assigning sentinels[...] to 1 fails, go to 2.
# 2. sentinels[...] is already 1, meaning that this row can be
# further reduced on the following iteration

# Sync point. Everything before this point becomes visible to
# all other threads once they reach this point.
# NOTE: Note the absense of a total ordering on atomic operations
old, success = Atomix.replace!(
Atomix.IndexableRef(sentinels, (Int(new_sparse_row_support[1]),)),
Int8(0),
Int8(1),
Atomix.acquire_release,
Atomix.acquire
)

@invariant length(new_sparse_row_support) == length(new_sparse_row_coeffs)

if success
# new pivot is formed, wake up from the inner loop
@invariant iszero(old)

linalg_normalize_row!(new_sparse_row_coeffs, arithmetic)

matrix.some_coeffs[i] = new_sparse_row_coeffs
matrix.lower_to_coeffs[new_sparse_row_support[1]] = i
pivots[new_sparse_row_support[1]] = new_sparse_row_support
else
# go for another iteration
first_nnz_column = new_sparse_row_support[1]
end
end

new_sparse_row_support, new_sparse_row_coeffs =
linalg_new_empty_sparse_row(CoeffType)
end

true
end

# Given a matrix of the following form, where A is in REF,
# A B
# C D
# reduces the lower part CD with respect to the pivots in the upper part AB.
# As a result, a matrix of the following form is produced:
# A B
# 0 D'
# The new pivots in the D' block are known, but possibly not fully interreduced.
function linalg_reduce_matrix_lower_part_threaded_cas_maybe_correct!(
matrix::MacaulayMatrix{CoeffType},
basis::Basis{CoeffType},
arithmetic::AbstractArithmetic{AccumType, CoeffType}
) where {CoeffType <: Coeff, AccumType <: Coeff}
if nthreads() == 1
@log level = 0 """
Using multi-threaded linear algebra with nthreads() == 1.
Something probably went wrong."""
end
_, ncols = size(matrix)
nup, nlow = matrix_nrows_filled(matrix)

# Prepare the matrix
pivots, row_index_to_coeffs = linalg_prepare_matrix_pivots!(matrix)
resize!(matrix.some_coeffs, nlow)
resize!(matrix.sentinels, ncols)

# If there is a pivot with the leading column i, then sentinels[i] = 1.
# If there is a pivot with the leading column i, then we set
# sentinels[i] = 1.
# If sentinels[i] = 1, and the pivot is synced, then we set
# sentinels[i] = 2.
# Otherwise, sentinels[i] = 0.
# Once sentinels[i] becomes 1, this equality is unchanged.
# NOTE: for all valid indices, isassigned(sentinels, i) must hold.
# Once sentinels[i] becomes 1 or 2, this equality is unchanged.
sentinels = matrix.sentinels
@inbounds for i in 1:ncols
sentinels[i] = 0
Expand Down
Loading

0 comments on commit 0ec642c

Please sign in to comment.