Julia intermediate array allocation with dots

user_2393329 0
Viewed 98

I have a PDE model in which I am trying to reduce allocations. Currently the following line is the biggest source of allocations:

part2 .= r.*u .+ deriv .+ q .+ u.^3

part2, deriv, and u are preallocated Float64 matrices. r and q are Float64 numbers.

My understanding was that the use of preallocated left hand side and . notation throughout should avoid allocation of intermediate temporary arrays, so why is this line still my biggest source of allocations?

If I replace this line with the following and add another preallocated matrix part1, I get zero allocations, but much less readable code.

broadcast!(*,part2,u,u,u)
broadcast!(*,part1,r,u)
broadcast!(+,part2,part2,part1,deriv,q)
julia
1 Answers

Unless your matrices are some non-standard objects that do some weird things the operation you have shown does not allocate:

julia> part2 = rand(1000, 1000);

julia> deriv = rand(1000, 1000);

julia> u = rand(1000, 1000);

julia> r, q = rand(2);

julia> using BenchmarkTools

julia> @btime $part2 .= $r .* $u .+ $deriv .+ $q .+ $u .^ 3;
  1.011 ms (0 allocations: 0 bytes)

(maybe you have used @time and seen the allocations that compilation causes when you run the code for the first time?)

user_12695670
Related