diff --git a/docs/src/sparse/intro.md b/docs/src/sparse/intro.md index 5de81737f1..c5922ba2e0 100644 --- a/docs/src/sparse/intro.md +++ b/docs/src/sparse/intro.md @@ -75,6 +75,14 @@ mod_sym!(::SRow{ZZRingElem}, ::Integer) maximum(::typeof(abs), ::SRow{ZZRingElem}) ``` +### Conversion to/from julia and AbstractAlgebra types + +```@docs +Vector(r::SRow, n::Int) +sparse_row(A::MatElem) +dense_row(r::SRow, n::Int) +``` + ## Sparse matrices Let $R$ be a commutative ring. Sparse matrices with base ring $R$ are modelled by @@ -207,6 +215,7 @@ Other: sparse(::SMat) ZZMatrix(::SMat{ZZRingElem}) ZZMatrix(::SMat{T}) where {T <: Integer} -Array(::SMat{T}) where {T} +Matrix(::SMat) +Array(::SMat) ``` diff --git a/src/NumField/NfRel/NfRelNS.jl b/src/NumField/NfRel/NfRelNS.jl index e607638761..909eb99190 100644 --- a/src/NumField/NfRel/NfRelNS.jl +++ b/src/NumField/NfRel/NfRelNS.jl @@ -496,16 +496,6 @@ function minpoly_dense(a::NfRelNSElem) end end -function Base.Matrix(a::SMat) - A = zero_matrix(base_ring(a), nrows(a), ncols(a)) - for i = 1:nrows(a) - for (k, c) = a[i] - A[i, k] = c - end - end - return A -end - function minpoly_sparse(a::NfRelNSElem) K = parent(a) n = degree(K) diff --git a/src/Sparse.jl b/src/Sparse.jl index d7a786ca54..f0db9f798f 100644 --- a/src/Sparse.jl +++ b/src/Sparse.jl @@ -17,7 +17,7 @@ One (+one with trafo) is probably enough =# -import Base.push!, Base.max, Nemo.nbits, Base.Array, +import Base.push!, Base.max, Nemo.nbits, Base.Array, Base.Matrix, Base.hcat, Base.vcat, Base.max, Base.min diff --git a/src/Sparse/Matrix.jl b/src/Sparse/Matrix.jl index 65a7d09216..54912301da 100644 --- a/src/Sparse/Matrix.jl +++ b/src/Sparse/Matrix.jl @@ -691,37 +691,36 @@ end # ################################################################################ -function *(b::T, A::SMat{T}) where {T <: RingElem} - B = sparse_matrix(base_ring(A), 0, ncols(A)) +function *(b::T, A::SMat{T}) where {T} if iszero(b) - return B + return sparse_matrix(base_ring(A), nrows(A), ncols(A)) end - for a = A - push!(B, b*a) + B = sparse_matrix(base_ring(A), 0, ncols(A)) + for a in A + push!(B, b * a) end return B end -function *(b::Integer, A::SMat{T}) where T - return base_ring(A)(b)*A -end - -function *(b::ZZRingElem, A::SMat{T}) where {T <: RingElement} - return base_ring(A)(b)*A +function *(b, A::SMat) + return base_ring(A)(b) * A end -function *(b::ZZRingElem, A::SMat{ZZRingElem}) +function *(A::SMat{T}, b::T) where {T} if iszero(b) - return zero_matrix(SMat, FlintZZ, nrows(A), ncols(A)) + return sparse_matrix(base_ring(A), nrows(A), ncols(A)) end - B = sparse_matrix(base_ring(A)) - B.c = ncols(A) + B = sparse_matrix(base_ring(A), 0, ncols(A)) for a in A - push!(B, b * a) + push!(B, a * b) end return B end +function *(A::SMat, b) + return A * base_ring(A)(b) +end + ################################################################################ # # Submatrix @@ -1396,19 +1395,34 @@ function SparseArrays.sparse(A::SMat{T}) where T return SparseArrays.sparse(I, J, V) end +@doc raw""" + Matrix(A::SMat{T}) -> Matrix{T} + +The same matrix, but as a julia matrix. +""" +function Matrix(A::SMat) + M = elem_type(base_ring(A))[zero(base_ring(A)) for _ in 1:nrows(A), _ in 1:ncols(A)] + for i in 1:nrows(A) + for (k, c) in A[i] + M[i, k] = c + end + end + return M +end + @doc raw""" Array(A::SMat{T}) -> Matrix{T} The same matrix, but as a two-dimensional julia array. """ -function Array(A::SMat{T}) where T - R = zero_matrix(base_ring(A), A.r, A.c) - for i=1:nrows(A) - for j=1:length(A.rows[i].pos) - R[i, A.rows[i].pos[j]] = A.rows[i].values[j] +function Array(A::SMat) + M = elem_type(base_ring(A))[zero(base_ring(A)) for _ in 1:nrows(A), _ in 1:ncols(A)] + for i in 1:nrows(A) + for (k, c) in A[i] + M[i, k] = c end end - return R + return M end #TODO: write a kronnecker-row-product, this is THE diff --git a/src/Sparse/Row.jl b/src/Sparse/Row.jl index 859e73298a..79a10c479f 100644 --- a/src/Sparse/Row.jl +++ b/src/Sparse/Row.jl @@ -1,4 +1,6 @@ -export sparse_row, dot, scale_row!, add_scaled_row, permute_row +import Base.Vector + +export sparse_row, dot, scale_row!, add_scaled_row, permute_row, dense_row ################################################################################ # @@ -474,8 +476,8 @@ function *(b::T, A::SRow{T}) where T if iszero(b) return B end - for (p,v) = A - nv = v*b + for (p,v) in A + nv = b*v if !iszero(nv) # there are zero divisors - potentially push!(B.pos, p) push!(B.values, nv) @@ -484,13 +486,35 @@ function *(b::T, A::SRow{T}) where T return B end -function *(b::Integer, A::SRow{T}) where T +function *(b, A::SRow) if length(A.values) == 0 return sparse_row(base_ring(A)) end return base_ring(A)(b)*A end +function *(A::SRow{T}, b::T) where T + B = sparse_row(parent(b)) + if iszero(b) + return B + end + for (p,v) in A + nv = v*b + if !iszero(nv) # there are zero divisors - potentially + push!(B.pos, p) + push!(B.values, nv) + end + end + return B +end + +function *(A::SRow, b) + if length(A.values) == 0 + return sparse_row(base_ring(A)) + end + return A*base_ring(A)(b) +end + function div(A::SRow{T}, b::T) where T B = sparse_row(base_ring(A)) if iszero(b) @@ -676,3 +700,49 @@ Returns the smallest entry of $A$. function minimum(A::SRow) return minimum(A.values) end + +################################################################################ +# +# Conversion from/to julia and AbstractAlgebra types +# +################################################################################ + +@doc raw""" + Vector(a::SMat{T}, n::Int) -> Vector{T} + +The first `n` entries of `a`, as a julia vector. +""" +function Vector(r::SRow, n::Int) + A = elem_type(base_ring(r))[zero(base_ring(r)) for _ in 1:n] + for i in intersect(r.pos, 1:n) + A[i] = r[i] + end + return A +end + +@doc raw""" + sparse_row(A::MatElem) + +Convert `A` to a sparse row. +`nrows(A) == 1` must hold. +""" +function sparse_row(A::MatElem) + @assert nrows(A) == 1 + if ncols(A) == 0 + return sparse_row(base_ring(A)) + end + return Hecke.sparse_matrix(A)[1] +end + +@doc raw""" + dense_row(r::SRow, n::Int) + +Convert `r[1:n]` to a dense row, that is an AbstractAlgebra matrix. +""" +function dense_row(r::SRow, n::Int) + A = zero_matrix(base_ring(r), 1, n) + for i in intersect(r.pos, 1:n) + A[1,i] = r[i] + end + return A +end diff --git a/src/Sparse/Solve.jl b/src/Sparse/Solve.jl index 889499ea9d..c025033cc2 100644 --- a/src/Sparse/Solve.jl +++ b/src/Sparse/Solve.jl @@ -1,3 +1,5 @@ +import Nemo: can_solve, can_solve_with_solution + function can_solve_ut(A::SMat{T}, g::SRow{T}) where T <: Union{FieldElem, zzModRingElem} # Works also for non-square matrices #@hassert :HNF 1 ncols(A) == nrows(A) @@ -346,6 +348,14 @@ function solve(a::SMat{T}, b::SRow{T}) where T <: FieldElem return sol end +function can_solve(a::SMat{T}, b::SRow{T}) where T <: FieldElem + c = sparse_matrix(base_ring(b)) + push!(c, b) + + # b is a row, so this is always from the left + return can_solve(a, c, side = :left) +end + function can_solve_with_solution(a::SMat{T}, b::SRow{T}) where T <: FieldElem c = sparse_matrix(base_ring(b)) push!(c, b) @@ -383,6 +393,28 @@ function find_pivot(A::SMat{T}) where T <: RingElement return p end +function can_solve(A::SMat{T}, B::SMat{T}; side::Symbol = :right) where T <: FieldElement + @assert side == :right || side == :left "Unsupported argument :$side for side: Must be :left or :right." + K = base_ring(A) + if side == :right + # sparse matrices might have omitted zero rows, so checking compatibility of + # the dimensions does not really make sense (?) + #nrows(A) != nrows(B) && error("Incompatible matrices") + mu = hcat(A, B) + ncolsA = ncols(A) + ncolsB = ncols(B) + else # side == :left + #ncols(A) != ncols(B) && error("Incompatible matrices") + mu = hcat(transpose(A), transpose(B)) + ncolsA = nrows(A) # They are transposed + ncolsB = nrows(B) + end + + rk, mu = rref(mu, truncate = true) + p = find_pivot(mu) + return !any(let ncolsA = ncolsA; i -> i > ncolsA; end, p) +end + function can_solve_with_solution(A::SMat{T}, B::SMat{T}; side::Symbol = :right) where T <: FieldElement @assert side == :right || side == :left "Unsupported argument :$side for side: Must be :left or :right." K = base_ring(A) diff --git a/test/Sparse/Matrix.jl b/test/Sparse/Matrix.jl index cba8d16190..47dd6c0b6b 100644 --- a/test/Sparse/Matrix.jl +++ b/test/Sparse/Matrix.jl @@ -177,19 +177,32 @@ using Hecke.SparseArrays E = @inferred 0 * D @test E == zero_matrix(SMat, FlintZZ, 3) @test E == sparse_matrix(FlintZZ, 3, 3) + E = @inferred D * 0 + @test E == zero_matrix(SMat, FlintZZ, 3) + @test E == sparse_matrix(FlintZZ, 3, 3) E = @inferred BigInt(2) * D @test E == sparse_matrix(FlintZZ, [2 10 6; 0 0 0; 0 2 0]) + E = @inferred D * BigInt(2) + @test E == sparse_matrix(FlintZZ, [2 10 6; 0 0 0; 0 2 0]) E = @inferred ZZRingElem(2) * D @test E == sparse_matrix(FlintZZ, [2 10 6; 0 0 0; 0 2 0]) + E = @inferred D * ZZRingElem(2) + @test E == sparse_matrix(FlintZZ, [2 10 6; 0 0 0; 0 2 0]) R = residue_ring(FlintZZ, 6) D = sparse_matrix(R, [1 2 2; 0 0 1; 2 2 2]) E = @inferred ZZRingElem(3) * D @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0]) + E = @inferred D * ZZRingElem(3) + @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0]) E = @inferred Int(3) * D @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0]) + E = @inferred D * Int(3) + @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0]) E = @inferred R(3) * D @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0]) + E = @inferred D * R(3) + @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0]) # Submatrix @@ -286,6 +299,8 @@ using Hecke.SparseArrays # Conversion to julia types D = sparse_matrix(FlintZZ, [1 5 3; 0 -10 0; 0 1 0]) + @test Matrix(D) == ZZRingElem[1 5 3; 0 -10 0; 0 1 0] + @test Array(D) == ZZRingElem[1 5 3; 0 -10 0; 0 1 0] E = SparseArrays.sparse(D) @test Matrix(E) == ZZRingElem[1 5 3; 0 -10 0; 0 1 0] @test Array(E) == ZZRingElem[1 5 3; 0 -10 0; 0 1 0] diff --git a/test/Sparse/Row.jl b/test/Sparse/Row.jl index cf3ccbe6fc..3a0cc489d4 100644 --- a/test/Sparse/Row.jl +++ b/test/Sparse/Row.jl @@ -100,7 +100,9 @@ for T in [Int, BigInt, ZZRingElem] b = T(2) B = @inferred b * A - @test B == map_entries(x -> T(2) * x, A) + @test B == map_entries(x -> b * x, A) + B = @inferred A * b + @test B == map_entries(x -> x * b, A) b = T(2) B = @inferred div(A, b) @@ -154,5 +156,11 @@ C = sparse_row(FlintZZ, [1, 2, 4, 5], ZZRingElem[-10, 100, 1, 1]) @test minimum(C) == ZZRingElem(-10) - + # Conversion + A = sparse_row(FlintZZ, [1, 3, 4, 5], ZZRingElem[-5, 2, -10, 1]) + @test Vector(A, 3) == ZZRingElem[-5, 0, 2] + @test Vector(A, 6) == ZZRingElem[-5, 0, 2, -10, 1, 0] + @test dense_row(A, 3) == matrix(FlintZZ, 1, 3, [-5, 0, 2]) + @test dense_row(A, 6) == matrix(FlintZZ, 1, 6, [-5, 0, 2, -10, 1, 0]) + @test sparse_row(dense_row(A, 6)) == A end diff --git a/test/Sparse/Solve.jl b/test/Sparse/Solve.jl index d429451d07..aaa3fbcd38 100644 --- a/test/Sparse/Solve.jl +++ b/test/Sparse/Solve.jl @@ -30,8 +30,10 @@ N = matrix(FlintQQ, rand([0,0,0,0,0,0,0,0,0,0,1], r, 2)) Ns = sparse_matrix(N) fl, sol = can_solve_with_solution(Ms, Ns, side = :right) + @test fl == can_solve(Ms, Ns, side = :right) @test nnz(sol) == sum(length(sol.rows[i].values) for i = 1:nrows(sol); init = 0) fl2, sol2 = can_solve_with_solution(M, N, side = :right) + @test fl2 == can_solve(M, N, side = :right) @test fl == fl2 if fl @test M*matrix(sol) == N