More functionality for sparse matrices and rows (#1221)

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)
 ```
 

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)

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
 

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

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

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)

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]

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

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