UniformScaling with SparseMatrixCSC{T,Int32} loses integer type

linearalgebra

#1

Is this a bug?

A = sparse(Int32[], Int32[], Float64[], 2, 2)
A + I # Loses Int32

A is of type SparseMatrixCSC{Float64,Int32} and A + I is of type SparseMatrixCSC{Float64,Int64}, losing the integer type.

Thanks


#2
julia> @which A+I
+(A::SparseMatrixCSC, J::UniformScaling) in SparseArrays at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/SparseArrays/src/sparsematrix.jl:3500

This is the code:

(+)(A::SparseMatrixCSC, J::UniformScaling) = A + sparse(J, size(A)...)

Now

julia> typeof(size(A))
Tuple{Int64,Int64}

julia> @which size(A)
size(S::SparseMatrixCSC) in SparseArrays at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/SparseArrays/src/sparsematrix.jl:36

size(S::SparseMatrixCSC) = (S.m, S.n)
struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Int                  # Number of rows
    n::Int                  # Number of columns

Maybe this should be:

struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Ti                  # Number of rows
    n::Ti                  # Number of columns

#3

I don’t think that’s enough because

sparse(I, Int32(3), Int32(3))

still has integer type Int64. Maybe it can be solved changing the definition of +, since

SparseMatrixCSC{Float64,Int32}(I, 3, 3)

works.

Should I open an issue? I think I can fix it.


#4

I don’t see how it is related to +.

julia> @which sparse(I, Int32(3), Int32(3))
sparse(s::UniformScaling, m::Integer, n::Integer) in SparseArrays at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/SparseArrays/src/sparsematrix.jl:1558

This is the code:

sparse(s::UniformScaling, dims::Dims{2}) = SparseMatrixCSC(s, dims)
sparse(s::UniformScaling, m::Integer, n::Integer) = sparse(s, Dims((m, n)))
help?> Dims
search: Dims ndims DimensionMismatch permutedims permutedims! PermutedDimsArray dropdims selectdim diagm spdiagm divrem dirname RoundingMode ENDIAN_BOM

  Dims{N}

  An NTuple of N Ints used to represent the dimensions of an AbstractArray.

julia> Dims
Tuple{Vararg{Int64,N}} where N

Maybe replacing

sparse(s::UniformScaling, dims::Dims{2}) = SparseMatrixCSC(s, dims)
sparse(s::UniformScaling, m::Integer, n::Integer) = sparse(s, Dims((m, n)))

by:

sparse(s::UniformScaling, m::Integer, n::Integer) = SparseMatrixCSC(s, promote(m,n))

and replacing in all constructors (part of concerned code unmodified below) of SparseMatrixCSC Integer by template Ti<:Integer:

SparseMatrixCSC{Tv,Ti}(s::UniformScaling, m::Integer, n::Integer) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(s, Dims((m, n)))
SparseMatrixCSC{Tv}(s::UniformScaling, m::Integer, n::Integer) where {Tv} = SparseMatrixCSC{Tv}(s, Dims((m, n)))
SparseMatrixCSC(s::UniformScaling, m::Integer, n::Integer) = SparseMatrixCSC(s, Dims((m, n)))
SparseMatrixCSC{Tv}(s::UniformScaling, dims::Dims{2}) where {Tv} = SparseMatrixCSC{Tv,Int}(s, dims)
SparseMatrixCSC(s::UniformScaling, dims::Dims{2}) = SparseMatrixCSC{eltype(s)}(s, dims)
function SparseMatrixCSC{Tv,Ti}(s::UniformScaling, dims::Dims{2}) where {Tv,Ti}
    @boundscheck first(dims) < 0 && throw(ArgumentError("first dimension invalid ($(first(dims)) < 0)"))
    @boundscheck last(dims) < 0 && throw(ArgumentError("second dimension invalid ($(last(dims)) < 0)"))
    iszero(s.λ) && return spzeros(Tv, Ti, dims...)
    m, n, k = dims..., min(dims...)
    nzval = fill!(Vector{Tv}(undef, k), Tv(s.λ))
    rowval = copyto!(Vector{Ti}(undef, k), 1:k)
    colptr = copyto!(Vector{Ti}(undef, n + 1), 1:(k + 1))
    for i in (k + 2):(n + 1) colptr[i] = (k + 1) end
    SparseMatrixCSC{Tv,Ti}(dims..., colptr, rowval, nzval)
end

But the fact that Dims is defined as an Int64 tuple may have other implications…


#5

I’m saying that instead of redefining SparseMatrixCSC{Tv,Ti} so that its sizes have type Ti, we simply change

(+)(A::SparseMatrixCSC, J::UniformScaling) = A + sparse(J, size(A)...) 

to something like

(+)(A::SparseMatrixCSC{Tv,Ti}, J::UniformScaling) where {Tv, Ti <: Integer} =
    A + SparseMatrixCSC{Tv,Ti}(J, size(A)...)

#6

But this would only work for +. What about -?
The rest of the code should be made consistent. The first thing bothering me is that Dims is defined with Int64 but SparseMatrix is allowed to have custom indices type. It should be made consistent.

julia> Dims
Tuple{Vararg{Int64,N}} where N

#7

Yes, - as well. And I + A also. All three are together in the code.


#8

While it might solve your particular problem, it may not be a definitive solution. There are other methods that use Int64 (eg via Dims) and the number of rows and cols being stored as Int64 makes it unusable for BigInt indices (which are not practical for dense stored arrays but may be useful for sparse arrays or lazy evaluation arrays).