SparseMatrixCSC * Vector not supported with StaticArray eltypes

[pablosanjose]

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.0-DEV.56 (2019-08-27)
 _/ |\__'_|_|_|\__'_|  |  Commit 9725fb4c3f (9 days old master)
|__/                   |

julia> using SparseArrays, StaticArrays

julia> a = rand(SMatrix{3,3,Float64,9}, 2, 2); v = rand(SVector{3,Float64}, 2);

julia> a * v
2-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [1.3476421609948477, 1.2721798108571298, 1.5096254759870544]
 [0.7505145227234096, 1.3942491617395882, 1.8458982908945654]

julia> sparse(a) * v
ERROR: MethodError: no method matching one(::Type{SArray{Tuple{3},Float64,1,3}})
Closest candidates are:
  one(::Type{Missing}) at missing.jl:103
  one(::BitArray{2}) at bitarray.jl:400
  one(::Missing) at missing.jl:100
  ...
Stacktrace:
 [1] *(::SparseMatrixCSC{SArray{Tuple{3,3},Float64,2,9},Int64}, ::Array{SArray{Tuple{3},Float64,1,3},1}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/linalg.jl:54
 [2] top-level scope at REPL[4]:1

The problem is more general than this specific case (mat * mat also fails with these eltypes) , and boils down to SparseArrays.mul! corresponding to the five-argument version mul!(C, A, B, α, β) (as opposed to the three-argument version for dense matrices). It is assumed in SparseArray.:* that T = promote_op(matprod, eltype(A), eltype(B)) has a α = one(T) and a β = zero(T) that are Numbers. In this case the one function is not even defined for a general SArrays, and it is actually never a Number. The zero(::SArray) is defined, but again is not a number but an SArray.

It seems to me that a good way to make this work in general is to simply replace one(T) by 1 and zero(T) by 0 (Ints, or perhaps true and false) in the following lines (and similar ones), and rely on type conversion in mul! to do the right thing

julia/stdlib/SparseArrays/src/linalg.jl

Lines 79 to 80 in 36c59c1

	*(adjA::Adjoint{<:Any,<:AbstractSparseMatrixCSC{TA,S}}, x::StridedVector{Tx}) where {TA,S,Tx} =
	(T = promote_op(matprod, TA, Tx); mul!(similar(x, T, size(adjA, 1)), adjA, x, one(T), zero(T)))

This change does not fail any tests locally for me and fixes this issue. I can prepare a PR if this approach is acceptable.

[pablosanjose]

I went ahead and started preparing a PR, as I need this for my work. I have a question though: how does one go about testing stuff in Base using non-stdlib packages such as StaticArrays? Can I just write "using StaticArrays" and test away?

[timholy]

Can I just write "using StaticArrays" and test away?

Unfortunately no. For OffsetArrays we developed a mini-package just for testing: https://github.com/JuliaLang/julia/blob/master/test/testhelpers/OffsetArrays.jl. The reason is that packages may evolve away from testing the thing that you originally wanted to test.

In this case I can replicate your error with

ad, vd = Matrix.(a), Vector.(v)

so I don't think this is a major limitation for you.

[pablosanjose]

Hi Tim,
there is a difference though: zero(::Type{Vector}) is not defined, while zero(::Type{SVector{N,T}}) is. Unfortunately this method (zero(eltype(C))) is used in the 5-arg mul! to initialize the target vector in the case β=0. Which means that the solution proposed above is not enough for your example (one would need an additional fix for eltypes without a zero method).

Regarding testing, I understand the problem. If I include a zero fix too, we can just test against Arrays as you suggest.

[timholy]

If I include a zero fix too

I don't think you can fix that. How long of a zero vector do you create? zero is supposed to create the additive identity. y + zero(y) works, but y + zero(typeof(y)) seemingly can't work in general unless the length is encoded in typeof(y).

[pablosanjose]

Yes, the idea would be to initialize the target vector C in a different way, not using zero(typeof(y)). One way is to do rmul!(C,0) instead of fill!(C, zero(eltype(C))). The problem is that we lose a tiny bit of performance that way...

EDIT: I'm talking about this spot

julia/stdlib/SparseArrays/src/linalg.jl

Lines 34 to 42 in 36c59c1

	function mul!(C::StridedVecOrMat, A::AbstractSparseMatrixCSC, B::Union{StridedVector,AdjOrTransStridedMatrix}, α::Number, β::Number)
	size(A, 2) == size(B, 1) || throw(DimensionMismatch())
	size(A, 1) == size(C, 1) || throw(DimensionMismatch())
	size(B, 2) == size(C, 2) || throw(DimensionMismatch())
	nzv = nonzeros(A)
	rv = rowvals(A)
	if β != 1
	β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
	end