UndefVarError in for loop

General Usage
22

It is an indicator to users that the function modifies its input arguments. Take a look at the Julia Style Guide for information on this and other stylistic conventions. For instance, NodeVoltage written that way indicates that it is a type. node_voltage would be preferred. None of this affects how the code functions though, so you can ignore it if you aren’t sharing your code.

1 Like
23

@Nathan_Boyer ! is an indicator to users that the function modifies at least one of its input arguments (since A_Mat is not being modified), correct?

24

Correct.

25

@Nathan_Boyer Thank you very much!
Could you please provide the version of your MWE code with mutation to let me fully understand the difference between the them?

26

Well, I meant something like below, but the number of allocations didn’t go down like I expected. You are reaching the limit of my understanding, so someone else may have to help from here.

function march_mutate!(NodeVoltage, A_Mat, I_Mat, V_Mat)
  dt = 1e-6
  for t in 0:dt:0.03 
    I_Mat[:] = curMatRLC(NodeVoltage, dt)
    V_Mat[:] = lu(A_Mat) \ I_Mat
    NodeVoltage[:,:] = reshape(V_Mat, 3, :)
  end
  return NodeVoltage, A_Mat, I_Mat, V_Mat
end

NodeVoltage  = zeros(3,5)
LNV = length(NodeVoltage)
A_Mat = rand(LNV, LNV)
I_Mat = zeros(LNV)
V_Mat = zeros(LNV)

julia> march!(NodeVoltage, A_Mat) == march_mutate!(NodeVoltage, A_Mat, I_Mat, V_Mat)
true

julia> @btime march!($NodeVoltage, $A_Mat)
  70.630 ms (210008 allocations: 84.69 MiB)

julia> @btime march_mutate!($NodeVoltage, $A_Mat, $I_Mat, $V_Mat)
  70.984 ms (210007 allocations: 84.69 MiB)
1 Like
27

A version which saves your time history could look something like below. There are ways to make it faster by pre-allocating everything.

function march_and_save(NodeVoltage, A_Mat, I_Mat, V_Mat)
  NodeVoltage_Vec = [NodeVoltage]
  A_Mat_Vec = [A_Mat]
  I_Mat_Vec = [I_Mat]
  V_Mat_Vec = [V_Mat]
  dt = 1e-6
  for t in 0:dt:0.03
    push!(I_Mat_Vec, curMatRLC(NodeVoltage, dt))
    push!(V_Mat_Vec, lu(A_Mat) \ I_Mat)
    push!(NodeVoltage_Vec, reshape(V_Mat, 3, :))
    push!(A_Mat_Vec, A_Mat)
  end
  return NodeVoltage_Vec, A_Mat_Vec, I_Mat_Vec, V_Mat_Vec
end

julia> @btime march_and_save($NodeVoltage, $A_Mat, $I_Mat, $V_Mat)
  74.990 ms (210067 allocations: 86.69 MiB)
2 Likes
28

@Nathan_Boyer thank you very much! this is what I was planning to do so. And yeah it is interesting to know the reason for the same number of allocations
@Jeff_Emanuel @jlperla @stevengj Could you please share with us the reason?

29

It doesn’t seem very complicated. Your 0:dt:03 loop has 30001 iterations. On each iteration, you are calling lu(A_Mat) (2 allocations: one for the resulting factorization and one for a workspace vector needed by LAPACK), and \ I_mat allocates another matrix for the result. Your curMatRLC function allocates two matrices (ones and the result of .±). Plus I think reshape has to allocate a small heap object (a “view”-like wrapper). This is at least 6 allocations per iteration, which corresponds to 6 * 30001 == 210007 allocations, accounting for 85% of the allocations. (I think there must be one more allocation in your inner loop that I’m missing, because 7 allocations per iteration would account for 99.97% of your allocations, with the remaining 60 allocations attributable to push! and the initial allocations of A_Mat_Vec etc.)

For example,

V_Mat[:] = lu(A_Mat) \ I_Mat

overwrites V_Mat, but it is not a fully in-place operation. lu(A_Mat) allocates as I noted above, and moreover lu(A_Mat) \ I_Mat allocates a new matrix on the right-hand side and then copies it over to V_Mat[:]. If you wanted to write the result directly into V_Mat you could use ldiv! instead.

1 Like
30

Hi Nathan,

good job!

31

For example, the following version of your march_mutate! function has only 2 allocations (from the lu call):

function march_mutate2!(NodeVoltage, A_Mat, V_Mat)
  dt = 1e-6
  LU = lu(A_Mat)
  for t in 0:dt:0.03
    V_Mat .= 1 + (NodeVoltage[1,1]==11 ? -dt : dt)
    ldiv!(LU, V_Mat)
    NodeVoltage[:] = V_Mat
  end
  return NodeVoltage, A_Mat, V_Mat
end

Note that (a) we construct the LU factorization only once and re-use it, (b) there is no need to construct an explicit ones array when you have .=, (c) we use ldiv! to compute the result in-place, and (d) there is no need for reshape.

(Since the right-hand-sides differ only by constant factors, you could skip the ldiv! entirely and just rescale the results, but I’m assuming this is a toy version of a real problem where the right-hand sides are more interesting.)

2 Likes
32

Thank you for your explanation.

Right, this is a toy version and in the real most of the matrices are changing.

33

I know the right hand sides are allocating a lot and there are many ways to improve the code. I didn’t want to change too much for fear of breaking the toy. I just thought f!(x) would allocate more than g!(x) when transferring to the left hand side.

function f!(x)
  for i in 1:10000
    x = ones(length(x))
  end
end

function g!(x)
  for i in 1:10000
    x[:] = ones(length(x))
  end
end

julia> x = fill(2, 1000)
1000-element Vector{Int64}:
 2
 2
 ⋮
 2
 2

julia> @btime f!(x)
  7.853 ms (10000 allocations: 77.51 MiB)

julia> @btime g!(x)
  20.816 ms (10000 allocations: 77.51 MiB)

But I guess there are no extra copies in f!. The difference is just where in memory the function places the result. The mutating version actually takes longer.

1 Like
34

@Nathan_Boyer
Yes, our point is about the allocations by mutation and re-assignment. I tried also in my real code but actually I was always getting the same number of allocations.

35

Nathan solved the original problem AFAIU. Shouldn’t you start a new thread for discussing reducing allocs?