generic_matvecmul complains about deprecated +

[acroy]

We noticed the problem today in SciML/ODE.jl#29. A boiled-down version is

B = Array(typeof([1.,2.]),2);
fill!(B,[1.,2.]);
Base.LinAlg.generic_matvecmul('N', [1. 1.], B)

which gives

WARNING: x::Number + A::Array is deprecated, use x .+ A instead.
 in + at deprecated.jl:26
 in generic_matvecmul at linalg/matmul.jl:310
 in generic_matvecmul at linalg/matmul.jl:274
1-element Array{Any,1}:
 [2.0,4.0]

The problem appears to be that generic_matvecmul (and also generic_matmatmul) uses zero(R) to initialize the variable holding the sum. R is some promoted type, which happens to be Any in the case above and zero(Any)==0. Adding the first product, which is a Vector, triggers the depreciation warning.

[JeffBezanson]

In this case the warning is asking us whether we really mean to take a product of a matrix of numbers and a vector of vectors. I don't know what to do here, but I think it's cool that there is a use case for this.

[acroy]

One possibility might be to use zero(eltype(A))*B[1] as a value for zero (A is the matrix and B the vector)? It would not be perfect, because the return type of generic_matvecmul is still Array{Any,1}, but at least there would be no warning.

[andreasnoack]

It's a bit tricky. Usually, I would infer the right return type by something like typeof(zero(eltype(A))+zero(eltype(B))) but that is not feasible when eltype(B) is a vector. It is possible to do something like what @acroy suggests, e.g. typeof(A[1]*B[1]), but then I cannot infer the return type when the arrays are empty. Maybe the size of arrays should be a parameter.

Another question is whether you really want to store your B as a vector of vectors. I don't know the ODE code, so maybe it is the best solution, but for the matvec product it appears as a slow solution.

[acroy]

Probably we shouldn't use matvec with vectors of vectors in ODE.jl, but that is a different question. I think it would be nice to have a generic matvecmul, which works for all types that support the necessary operations (eltype(A) * eltype(B), zero and + for the product type), even if it is slow.

Why is typeof(zero(eltype(A))*zero(eltype(B))) not feasible for eltype(B) being a vector? Would it be problematic to define zero{T}(x::Vector{T})=zeros(eltype(T), size(x)) (or something better)?

[andreasnoack]

I agree. If possible, if would be great to support fully generic matvecmul and matmatmul. The problem is that zero(Vector) doesn't make sense because the size is only specified for an instanceof a type, not the type itself. When B is empty, there is no instance of the type to use for zero, but only information about the type, i.e. eltype(B).

[acroy]

I see. Maybe one could have a ZeroArray type similar to the new UniformScaling?

[andreasnoack]

I see your point, but would prefer some other solution because a ZeroVector doesn't seem that useful to me. For higher dimensions it is also not obvious what the size should be.

I realise that it there is an issue for FixedArrays so maybe that is the way to go.

[andreasnoack]

With something like this

import Base: convert, getindex, setindex!, similar, size, zero
type FixedSizeArray{T,N,S} <: DenseArray{T,N}
    data::Array{T,N}
end
FixedSizeArray{T,N}(A::Array{T,N}) = FixedSizeArray{T,N,size(A)}(A)

convert(T::Type, A::FixedSizeArray) = convert(T, A.data)
convert{T,N,S}(::Type{FixedSizeArray{T,N,S}}, A::Array{T,N}) = FixedSizeArray(A)

similar(A::FixedSizeArray) = FixedSizeArray(similar(A.data))
similar{T}(A::FixedSizeArray, ::Type{T}) = FixedSizeArray(similar(A.data, T))
similar{T}(A::FixedSizeArray, ::Type{T}, n::Dims) = FixedSizeArray(similar(A.data, T, n))

getindex(A::FixedSizeArray,i::Integer) = getindex(A.data, i)
getindex(A::FixedSizeArray,i::Range1{Int}) = getindex(A.data, i)
getindex(A::FixedSizeArray,i::Integer,j::Integer) = getindex(A.data,i,j)

setindex!(A::FixedSizeArray, x, i::Integer) = setindex!(A.data, x, i)
setindex!(A::FixedSizeArray, x, i::Integer, j::Integer) = setindex!(A.data, x, i, j)

size{T,N,S}(A::FixedSizeArray{T,N,S}) = S

zero{T,N,S}(::Type{FixedSizeArray{T,N,S}}) = FixedSizeArray(zeros(T,S))

you can almost get what you want. A slightly modifies version of your example is

julia> B=Array(FixedSizeArray{Float64,1,(2,)},2);
julia> fill!(B,FixedSizeArray([1.,2.]));
julia> Base.LinAlg.generic_matvecmul(similar(B,1),'N', [1. 1.], B)
1-element Array{FixedSizeArray{Float64,1,(2,)},1}:
 [2.0,4.0]

[acroy]

Thanks! FixedArrays might actually be a (performant) solution for us - I will look into it.

For the general problem here, I would say that if the underlying issue is that zero(Vector) is not defined, there should be an according error message. Changing the return type to typeof(zero(eltype(A))*zero(eltype(B))) is probably the right thing, this is also consistent with requiring +, * and zero for matvecmul/matmatmul to work. (BTW: Why is actually zero(Any)==0?)

[andreasnoack]

I agree, but it would actually have to be something like typeof(zero(eltype(A))*zero(eltype(B))+zero(eltype(A))*zero(eltype(B))) because Bool*Bool=Bool but Bool+Bool=Int.

The definition of zero(Any) is in operators.jl where many methods are defined for Any and I guess it has just been convenient to have the fallbacks like that. I also agree with you on this and I think zero(Any) shouldn't be defined. Restricting the method to Number would be better, I think, even though you could argue for restricting it further.

[andreasnoack]

Deprecation of Array+Number was reverted.