Skip to content

Commit

Permalink
some improvement with pointers
Browse files Browse the repository at this point in the history
  • Loading branch information
@BD-00
BD-00 committed Nov 28, 2024
1 parent 445739d commit e25a850
Show file tree
Hide file tree
Showing 2 changed files with 182 additions and 16 deletions.
32 changes: 17 additions & 15 deletions src/Sparse/StructuredGauss.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ mutable struct data_StructGauss{T}
end
end

function _col_list(A::SMat{T}) where T
function _col_list(A::SMat{T})::Vector{Vector{Int64}} where T
n = nrows(A)
m = ncols(A)
col_list = [Array{Int64}([]) for i = 1:m]
Expand All @@ -76,7 +76,7 @@ Return a tuple (\nu, N) consisting of the nullity \nu of A and a basis N
AN is the zero matrix.
"""

function structured_gauss(A::SMat{T}) where T <: ZZRingElem
function structured_gauss(A::SMat{T}) where T <: RingElem
SG = part_echolonize!(A)
return compute_kernel(SG)
end
Expand All @@ -103,7 +103,7 @@ end
#Build an upper triangular matrix for as many columns as possible compromising
#the loss of sparsity during this process.

function part_echolonize!(A::SMat{T}) where T <: ZZRingElem
function part_echolonize!(A::SMat{T})::data_StructGauss where T <: RingElem
A = delete_zero_rows!(A)
n = nrows(A)
m = ncols(A)
Expand All @@ -129,7 +129,7 @@ function part_echolonize!(A::SMat{T}) where T <: ZZRingElem
return SG
end

function part_echolonize_field!(A)
function part_echolonize_field!(A::SMat{T})::data_StructGauss where T <: FieldElem
A = delete_zero_rows!(A)
n = nrows(A)
m = ncols(A)
Expand All @@ -156,7 +156,7 @@ function part_echolonize_field!(A)
return SG
end

function single_rows_to_top!(SG::data_StructGauss)
function single_rows_to_top!(SG::data_StructGauss)::data_StructGauss
for i = 1:nrows(SG.A)
len = length(SG.A[i])
@assert !iszero(len)
Expand Down Expand Up @@ -230,7 +230,7 @@ function build_part_ref_field!(SG::data_StructGauss)
end
end

function find_best_single_row(SG::data_StructGauss)
function find_best_single_row(SG::data_StructGauss)::Int64
best_single_row = -1
best_col = NaN
best_val = NaN
Expand Down Expand Up @@ -268,7 +268,7 @@ function find_best_single_row(SG::data_StructGauss)
return best_single_row
end

function find_best_single_row_field(SG::data_StructGauss)
function find_best_single_row_field(SG::data_StructGauss)::Int64
best_single_row = -1
best_len = -1
for i = SG.base:SG.single_row_limit-1
Expand Down Expand Up @@ -338,7 +338,7 @@ function turn_heavy(SG::data_StructGauss)
SG.nlight -= SG.heavy_ext
end

function handle_new_light_weight!(i::Int64, SG::data_StructGauss)
function handle_new_light_weight!(i::Int64, SG::data_StructGauss)::data_StructGauss
w = SG.light_weight[i]
if w == 0
swap_i_into_base(i, SG)
Expand All @@ -354,7 +354,7 @@ function handle_new_light_weight!(i::Int64, SG::data_StructGauss)
return SG
end

function eliminate_and_update!(best_single_row::Int64, SG::data_StructGauss)
function eliminate_and_update!(best_single_row::Int64, SG::data_StructGauss)::data_StructGauss
@assert best_single_row != 0
best_row = deepcopy(SG.A[best_single_row])
best_col = find_light_entry(best_row.pos, SG.is_light_col)
Expand All @@ -376,15 +376,15 @@ function eliminate_and_update!(best_single_row::Int64, SG::data_StructGauss)
return SG
end

function eliminate_and_update_field!(best_single_row::Int64, SG::data_StructGauss)
function eliminate_and_update_field!(best_single_row::Int64, SG::data_StructGauss)::data_StructGauss
@assert best_single_row != 0
best_col = find_light_entry(SG.A[best_single_row].pos, SG.is_light_col)
@assert length(SG.col_list[best_col]) > 1
best_val = SG.A[best_single_row, best_col]
@assert !iszero(best_val)
scale_row!(SG.A, best_single_row, inv(best_val))
best_val = SG.A[best_single_row, best_col]
@assert isone(best_val)A.nnz
@assert isone(best_val)
best_row = deepcopy(SG.A[best_single_row])
best_col_idces = SG.col_list[best_col]
row_idx = 0
Expand All @@ -401,7 +401,7 @@ function eliminate_and_update_field!(best_single_row::Int64, SG::data_StructGaus
return SG
end

function find_row_to_add_on(row_idx::Int64, best_row::SRow{T}, best_col_idces::Vector{Int64}, SG::data_StructGauss) where T
function find_row_to_add_on(row_idx::Int64, best_row::SRow{T}, best_col_idces::Vector{Int64}, SG::data_StructGauss)::Int64 where T <: FieldElem
for L_row in best_col_idces[end:-1:1]
row_idx = SG.col_list_permi[L_row]
SG.A[row_idx] == best_row && continue
Expand All @@ -411,7 +411,7 @@ function find_row_to_add_on(row_idx::Int64, best_row::SRow{T}, best_col_idces::V
return row_idx
end

function add_to_eliminate!(L_row::Int64, row_idx::Int64, best_row::SRow{T}, best_col::Int64, best_val::ZZRingElem, SG::data_StructGauss) where T <: ZZRingElem
function add_to_eliminate!(L_row::Int64, row_idx::Int64, best_row::SRow{T}, best_col::Int64, best_val::RingElem, SG::data_StructGauss)::data_StructGauss where T <: RingElem
@assert L_row in SG.col_list[best_col]
@assert !(0 in SG.A[row_idx].values)
val = SG.A[row_idx, best_col]
Expand All @@ -437,7 +437,7 @@ function add_to_eliminate!(L_row::Int64, row_idx::Int64, best_row::SRow{T}, best
return SG
end

function add_to_eliminate_field!(L_row::Int64, row_idx::Int64, best_row::SRow{T}, best_col::Int64, best_val, SG::data_StructGauss) where T
function add_to_eliminate_field!(L_row::Int64, row_idx::Int64, best_row::SRow{T}, best_col::Int64, best_val::FieldElem, SG::data_StructGauss)::data_StructGauss where T <: FieldElem
@assert L_row in SG.col_list[best_col]
@assert !(0 in SG.A[row_idx].values)
val = SG.A[row_idx, best_col]
Expand All @@ -451,12 +451,14 @@ function add_to_eliminate_field!(L_row::Int64, row_idx::Int64, best_row::SRow{T}
end
end
@assert !(0 in SG.A[row_idx].values)
SG.A.nnz -= length(SG.A[row_idx])
Hecke.add_scaled_row!(best_row, SG.A[row_idx], -val)
SG.A.nnz += length(SG.A[row_idx])
@assert iszero(SG.A[row_idx, best_col])
return SG
end

function update_after_addition!(L_row::Int64, row_idx::Int64, best_col::Int64, SG::data_StructGauss)
function update_after_addition!(L_row::Int64, row_idx::Int64, best_col::Int64, SG::data_StructGauss)::data_StructGauss
SG.light_weight[row_idx] = 0
for c in SG.A[row_idx].pos
@assert c != best_col
Expand Down
166 changes: 165 additions & 1 deletion src/Sparse/ZZStructuredGauss.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#TODO: sizehint for heavy stuff can be given later, so don't initialize immediately
#TODO: new struct for second part
#TODO: add sizehint for Y?
#TODO: remove whitespaces at beginning of lines

#=
PROBLEMS:
Expand Down Expand Up @@ -381,4 +382,167 @@ mutable struct data_ZZStructGauss{T}
end
SG.A.nnz-=length(SG.A[SG.base])
empty!(SG.A[SG.base].pos), empty!(SG.A[SG.base].values)
end
end

function delete_zero_rows!(A::SMat{T}) where T <: ZZRingElem
for i = A.r:-1:1
if isempty(A[i].pos)
deleteat!(A.rows, i)
A.r-=1
end
end
return A
end

################################################################################
#
# Kernel Computation
#
################################################################################

#Compute the kernel corresponding to the non echolonized part from above and
#insert backwards using the triangular part to get the full kernel.

mutable struct data_ZZKernel
heavy_mapi::Vector{Int64}
heavy_map::Vector{Int64}

function data_ZZKernel(SG::data_ZZStructGauss, nheavy::Int64)
return new(
sizehint!(Int[], nheavy),
fill(0, ncols(SG.A))
)
end
end

function compute_kernel(SG, with_light = true)
Hecke.update_light_cols!(SG)
@assert SG.nlight >= 0
KER = Hecke.collect_dense_cols!(SG)
D = dense_matrix(SG, KER)
_nullity, _dense_kernel = nullspace(D)
l, K = Hecke.init_kernel(_nullity, _dense_kernel, SG, KER, with_light)
return compose_kernel(l, K, SG)
end

function update_light_cols!(SG)
for i = 1:ncols(SG.A)
if SG.is_light_col[i] && !isempty(SG.col_list[i])
SG.is_light_col[i] = false
SG.nlight -= 1
end
end
return SG
end

function collect_dense_cols!(SG)
m = ncols(SG.A)
nheavy = m - SG.nlight - SG.nsingle
KER = data_ZZKernel(SG, nheavy)
j = 1
for i = 1:m
if !SG.is_light_col[i] && !SG.is_single_col[i]
KER.heavy_map[i] = j
push!(KER.heavy_mapi, i)
j+=1
end
end
@assert length(KER.heavy_mapi) == nheavy
return KER
end

function dense_matrix(SG, KER)
D = ZZMatrix(nrows(SG.A) - SG.nsingle, length(KER.heavy_mapi))
for i in 1:length(SG.Y)
for j in 1:length(KER.heavy_mapi)
setindex!(D, SG.Y[i], i, j, KER.heavy_mapi[j])
end
end
return D
end

function dense_kernel(SG, KER)
ST = sparse_matrix(base_ring(SG.A), 0, nrows(SG.Y))
YT = transpose(delete_zero_rows!(SG.Y))
for j in KER.heavy_mapi
push!(ST, YT[j])
end
S = transpose(ST)
d, _dense_kernel = nullspace(matrix(S))
return size(_dense_kernel)[2], _dense_kernel
end

function init_kernel(_nullity, _dense_kernel, SG, KER, with_light=false)
R = base_ring(SG.A)
m = ncols(SG.A)
if with_light
l = _nullity+SG.nlight
else
l = _nullity
end
K = ZZMatrix(m, l)
#initialisation
ilight = 1
for i = 1:m
if SG.is_light_col[i]
if with_light
@assert ilight <= SG.nlight
K[i, _nullity+ilight] = one(R)
ilight +=1
end
else
j = KER.heavy_map[i]
if j>0
for c = 1:_nullity
ccall((:fmpz_set, Nemo.libflint), Cvoid, (Ptr{ZZRingElem}, Ptr{ZZRingElem}), Nemo.mat_entry_ptr(K, i, j), Nemo.mat_entry_ptr(_dense_kernel, j, c))
#K[i,c] = _dense_kernel[j, c]
end
end
end
end
return l, K
end

function compose_kernel(l, K, SG)
R = base_ring(SG.A)
n = nrows(SG.A)
for i = n:-1:1
c = SG.single_col[i]
if c>0
x = R(0)
y = zeros(R,l)
for idx in 1:length(SG.A[i])
cc = SG.A[i].pos[idx]
xx = SG.A[i].values[idx]
if cc == c
x = xx
else
for k = 1:l
kern_c = K[cc,k]
!iszero(kern_c) && (y[k]-=xx*kern_c)
end
end
end
for k = 1:l
x_inv = try
inv(x)
catch
R(0)
end
if iszero(x_inv)
K[:,k] *= x
K[c, k] = y[k]
else
K[c, k] = y[k]*x_inv
end
end
end
end
return l, K
end

#Warning: skips zero entries in sparse matrix
function Base.setindex!(A::ZZMatrix, B::SRow{ZZRingElem}, ar::Int64, ac::Int64, bj::Int64)
isnothing(findfirst(isequal(bj), B.pos)) && return
ccall((:fmpz_set, Nemo.libflint), Cvoid, (Ptr{ZZRingElem}, Ptr{Int}), Nemo.mat_entry_ptr(A, ar, ac), Hecke.get_ptr(B.values, findfirst(isequal(bj), B.pos)))
end

0 comments on commit e25a850

Please sign in to comment.