Well there’s your problem - b can hold elements of quite literally any type and there is no one definition of oneunit for all possible types. I don’t know how you generate b, but you’ll probably want to make sure it contains elements of a concrete type, in your case possibly Float64.
b is constructed as a matrix product of SymTridiagonal{Float64,Array{Float64,1}} and Array{Real,1} and the result was an array of Any which… should I think that’s strange?
Anyway, I just wrapped the definition of b with real, i.e.
Maybe there’s a missing method? If you can share some more code to reproduce this deterministically, you might be able to file an issue to get this fixed (or you’re doing something else weird, who knows :p).
I see your point, and the answer is not necessarily yes, but considering things like optimization through type stability, and… other stuff I’m sure exists in efficient and healthy coding in Julia… I suppose the answer is no.
I had the (incorrect) mindset that if I initialize with Real then it would well handle anything I feed the vector elements.
Oh, in my original context, there were specific values to assign the vector in a loop (complicated so I couldn’t use comprehensions (or could I…)) and in my MWE above, I used rand to grab any random values, what the values were didn’t matter.
You can always put the hard work of generating those values into a function and call that in the comprehension - in fact, it’ll probably work better than hard-coding Real since the compiler can infer exactly what types are actually needed (as long as the generating function is somewhat type stable, of course)
Nevertheless, I think the compiler should be able to infer Real here
Well… not if it unnecessarily restricts the type, like in your case If you change your mind or the return type of the function called in your comprehension changes, your code will recompile to fit those types and still be fast. If you keep it abstractly typed, you’ll incur performance overhead from having to actually check the runtime type of the elements.
No, that’s not what I mean In fact, “restricts” was the wrong word to use, it’s more like writing Array{Real,1} unnecessarily opens your the type of the array? As a result of that, A*v (wrongly) infers to Any, making everything down the road slower. Since you’re only (presumably) storing a small subset of all possible concrete types of Real at any one point, having the container be Real instead of e.g. Union{Float64,Int64} like it may be as a result of using the [ f(...) for _ in 1:n ] approach limits optimization to things that any Real can be. If it were Union{Float64,Int64}, the compiler would know that each element was 64 bits long and could make smarter decisions about how to make it fast. Otherwise it would have to take care of the cases where the element is e.g. Int8, which is much shorter with only 8 bits.
I’ve subscribed to the issue you opened, I’m interested to see where this goes. But AFAI understand, what you’re suggesting totally makes sense - get the inference right so that we can continue writing code that will be automatically type stable, instead of becoming Any.
julia> methods(+, (Int, Real))
# 7 methods for generic function "+":
[1] +(x::T, y::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at int.jl:87
[2] +(c::Union{Int16, Int32, Int64, Int8}, x::BigInt) in Base.GMP at gmp.jl:534
[3] +(c::Union{Int16, Int32, Int64, Int8}, x::BigFloat) in Base.MPFR at mpfr.jl:385
[4] +(a::Integer, b::Integer) in Base at int.jl:919
[5] +(x::T, y::T) where T<:Number in Base at promotion.jl:394
[6] +(y::Integer, x::Rational) in Base at rational.jl:295
[7] +(x::Number, y::Number) in Base at promotion.jl:321
julia> Base.return_types(+, (Int, Real))
7-element Vector{Any}:
Int64
BigInt
BigFloat
Any
Union{}
Rational
Any
Here there are 5 concrete and 2 abstract applicable methods, well over the limit of 3.
Intuitively that’s what I initially thought would happen, with the philosophy that it would/should “just work”. But I also accept that there’s more going on within the internals, and will leave it to you guys to figure out a viable solution.
v = Vector{Real}(undef, 3)
for n = 1:length(v)
v[n] = rand(Float64)
end
is
v = [rand(Float64) for n in 1:3]
(We’re assuming that rand(Float64) is a stand-in for some more complicated expression, possibly a function of n, etc.)
The latter syntax not only saves a bit of typing, it will also give a more narrow element type for v. In fact, even if you do use the first syntax to generate v, you could re-generate it afterwards using
(as long as the generating function is somewhat type stable, of course)
If you change your mind or the return type of the function called in your comprehension changes, your code will recompile to fit those types and still be fast. If you keep it abstractly typed, you’ll incur performance overhead from having to actually check the runtime type of the elements.