Help me understand vector initialization

General Usage Performance
1

Here is a MWE of a problem that I am facing. In method 1, the array “arr” gets initialized correctly, but it seems very inelegant. In the second method, the array “brr” often does not get initialized correctly.

for k = 1:100
    local i = 3
    local arr = Vector{Float64}(undef, i)
    local brr = Vector{Float64}(undef, i)
    for j = 1:i
        arr[j] = 0
    end
    brr .*= 0 #<- this should have worked but doesn't
    arr[i] = 1
    brr[i] = 1
    if all(isfinite, arr) == false || all(isfinite, brr) == false
        print("$arr v/s $brr \n")
    end
end

On average, about 4 times out of 100, the output is like so:

[0.0, 0.0, 1.0] v/s [NaN, 0.0, 1.0]

What’s going on here?

2

you want brr .= 0 NaN*0==NaN.

1 Like
3

Thanks very much!

4

While there is nothing wrong with the other answer, I think it is rare that you’d want to really use the undef initializer. Bugs coming from a slightly wrong use of undef initializers can occur infrequently (because more often than not, the memory will be full of zeros), they can fail very silently and are in general hard to reproduce and debug.
Perhaps, the zeros function is not suitable for your actual use can but perhaps, fill(0,i) is?. 0 can of course be interchanges with basically anything.
This way of initialization, along with array comprehensions have completely replaces my use of undef initializers.
Not saying that there can never be a use case, but it always eases my mind when I manage to avoid such pitfalls.

2 Likes
5

One of the issues with the 0 * NaN behaviour is that BLAS methods don’t tend to respect this, i.e.
LinearAlgebra.BLAS.scal!(0., x)
eliminates possible NaN values in x. I thought it used to be the case that this BLAS method was called by more intrinsic Julia methods like LinearAlgebra.lmul!, but this seems to be no longer the case, as this method does preserve the NaNs correctly. However, doing something like LinearAlgebra.axpy!(0., x, y) where x contains NaNs will still discard those and not propagate them into y.

1 Like
6

Thanks for your note. Yes, fill works; specifically, fill!(brr,0.0) does the job on an existing array and brr = fill(0.0,i) to initialize a new array. Not sure about the difference in performance, though. Yes, the zeros function doesn’t work for me because of how I use the array next (passing it as an argument to the ChebyshevT function) which wants a vector, not Matrix, as input.

7

Why would zeros not be able to create a Vector?

julia> zeros(5)
5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0

julia> zeros(Int, 5)
5-element Vector{Int64}:
 0
 0
 0
 0
 0
1 Like
8

Ah yes, that works too. As you can tell, I am new to Juila, coming from a MATLAB background, where zeros(n) gives a square matrix of size n x n. Thanks for pointing this out.