I wrote a small test code that changes matrix elements and stores them in the “old” matrix, so no new allocations should be made.
I then compiled the code with @allocated to measure the allocations made:
 function test(N)
80080 M = Array{BigFloat}(N,N)
0 for i = 1:N
0 for j = 1:N
0 M[i,j] = big(pi) #spawn matrix elements
 end
 end

 #Now change matrix elements
0 for i = 1:N
0 for j = 1:N
0 M[i,j] = BigFloat(2)*M[i,j]
 end
 end
 end

As you can see, the code behaves as expected with only allocations made when the matrix is initialized.
However, when I add multi threading, new allocations are being made:
 function test(N)
80096 M = Array{BigFloat}(N,N)
0 Threads.@threads for i = 1:N
1176 Threads.@threads for j = 1:N
174744 M[i,j] = big(pi) #spawn matrix elements
 end
 end

 #Now change matrix elements
0 Threads.@threads for i = 1:N
1640 Threads.@threads for j = 1:N
773132 M[i,j] = BigFloat(2)*M[i,j]
 end
 end
 end

Can anyone explain this behavior and tell me how to bypass it?
Thanks in advance
EDIT:
To give you a quick update:
I tried to make my algorithm work by only using Arbs, to avoid getting new allocations by the BigFloats.
Unfortunately, this does not work because Arbs seem to not be made for this kind of usage.
They lose numerically stability during the matrix operations and get unprecise results (unless one puts the precision much higher than for BigFloats, but this makes the performance even worse than the BigFloatallocations).
I would like to know if someone knows another way how I could make this type of algorithm work.
Maybe there is another multiprecision floating type or a way to preallocate the BigFloat calculations?