diff --git a/src/Sparse/StructuredGauss.jl b/src/Sparse/StructuredGauss.jl index 7e5e210e3b..9bb3cbd1b5 100644 --- a/src/Sparse/StructuredGauss.jl +++ b/src/Sparse/StructuredGauss.jl @@ -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] @@ -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 @@ -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) @@ -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) @@ -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) @@ -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 @@ -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 @@ -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) @@ -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) @@ -376,7 +376,7 @@ 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 @@ -384,7 +384,7 @@ function eliminate_and_update_field!(best_single_row::Int64, SG::data_StructGaus @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 @@ -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 @@ -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] @@ -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] @@ -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 diff --git a/src/Sparse/ZZStructuredGauss.jl b/src/Sparse/ZZStructuredGauss.jl index 91dba17425..4724b20b02 100644 --- a/src/Sparse/ZZStructuredGauss.jl +++ b/src/Sparse/ZZStructuredGauss.jl @@ -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: @@ -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 \ No newline at end of file + 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 \ No newline at end of file