Are you using that function or the in-place version? This one
I think what I am doing is equivalent to the one you pointed out, this is what I am doing:
cholB1KtB1ff = let Ktu = makeKt(Φ, elems, u, d) B1KtB1ff = Symmetric(transpose(B1[:,ifree])*(λT/Δt^2*MM+λV*Ktu)*B1[:,ifree]) cholesky(B1KtB1ff) end
this should mean that neither
B1KtB1ff are kept, and memory should be freed after the
end, is that correct ?
would this be any different?
cholB1KtB1ff = let Ktu = makeKt(Φ, elems, u, d) cholesky!(Symmetric(transpose(B1[:,ifree])*(λT/Δt^2*MM+λV*Ktu)*B1[:,ifree])) end
MM has the same dimensions of
Ktu, is assembled once and never changes, and the other are scalars.
If you want specific help it would probably be best if you included some self-contained code. As far as general tips go, check out the Performance tips section in the manual. In particular, try doing allocation profiling (very convenient in the (VS)Code editor, just do
The latest Julia release also has a feature, seemingly still not documented, to trigger frequent garbage collection after a certain threshold. For example:
If I remember correctly, sometimes it was necessary to invoke this call multiple times in a row.
No, Julia isn’t C++. The memory gets freed when the garbage collector decides so.
I think I am following all of the guidelines of Performance tips, especially when it comes to enclosing code in functions, reusing vectors, and not having globals, I will look at
@profview_allocs some_call()) to see if I can get a better understanding of what is happening
GC.gc() multiple times, at the moment it is done at every increment of the boundary conditions, which is more often than the stiffness matrix factorization gets updated. But still the memory footprint of the program is from 50 to 100 times the expected
Just remember to run with
@profview_allocs an additional time before doing it for real, to get everything to compile.
The factorization objects (FO) returned by
cholesky contain large arrays managed by the solver library (so not directly visible to Julia tooling). When an FO goes out of scope and is garbage-collected, the GC should call a finalizer which frees the underlying storage. Are you sure that your obsolete FO’s are going out of scope?
I was thinking something like that was happening, but I was not aware of this peculiarity with the garbage collection of the FO objects,
in fact the memory occupancy of
Base.summarysize is 8 bytes, while the size in memory of the matrix is \approx 150 mb, that suggested that the decomposition was handled and kept somewhere out of the reach of julia, but thought it was still handled by
can you suggest how to do it, or where could I find information about it?
do you think it would make any difference if I switched to
MKL instead of
OpenBLAS that comes with the standard Julia installation?
Sorry if I wasn’t clear; the
gc call should automatically free the memory by calling the finalizer without extra effort from you, as long as you don’t keep the
cholB1KtB1ff variable (which is just a wrapped pointer) in scope. You might have been misled about the need to drop references to the FO’s by the small size visible to the Julia tools.
Finalizers are described in the Julia manual, and they are set up for the sparse array factorizations in the
SharedArrays interface to the
SuiteSparse library. Which BLAS library you use should not matter here.
I am keeping
cholB1KtB1ff in the main function because I need it multiple times, and I overwrite it when I update it, what should I do to make sure that
gc() takes care of disposing the old one when a new one is made?
cholB1KtB1ff was a pointer to something living elsewhere, but I thought
gc could reach it even if it was not visible through
I meant reusing the input buffer, with:
cholesky!(A::AbstractMatrix, NoPivot(); check = true) -> Cholesky
The same as [`cholesky`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky), but saves space by overwriting the input `A`, instead of creating a copy. An [`InexactError`](https://docs.julialang.org/en/v1/base/base/#Core.InexactError) exception is thrown if the factorization produces a number not representable by the element type of `A`, e.g. for integer types.
Also note that
A + B will create a new Matrix as a result. If you are doing this iteratively, you probably should compute
C .= A .+ B and reuse
C. Matrix multiplications can be done in place with
B[:,ifree] allocates a new vector. Use, perhaps,
I don’t know how much each of these things might make a difference there for you.
I was aware of that, but the old copies should be discharged by the
gc(), right? I was not too much worried about that because it takes place only a limited number of times, and it is not where the bulk of the computations lie
Although what you describe seems reasonable to me, there may be some subtlety about scope and object references in play here. Instead of doing so much in your main function, perhaps you could have an intermediate function which does the factorization of a given matrix, applies it to all the vectors associated with that problem, and then returns, thus isolating each
cholB1KtB1ff. Such function barriers sometimes help the memory management.
thanks, I’ll try that and I’ll post the results
Note that OP is using sparse matrices, for which different conventions apply. You can do this:
F = cholesky(A) # use F cholesky!(F, B) # use F
From the original post, it sounds like this is what OP is doing (“I also store the cholesky factorization for re-use”), but this is not what’s happening in the example code provided later.
I’m not an expert here, but in general, I think such reuse is only advisable if
B have the same sparsity pattern. However, if that’s the case, it’s highly recommended, as it not only saves allocations but also saves you from recomputing the symbolic factorization.
BLAS libraries don’t matter here. Sparse linear algebra is handled by different libraries altogether. Cholesky factorization and solving linear systems using the factorization object are both done by SuiteSparse.CHOLMOD, and matrix-vector multiplication is handled by a Julia native function within the SparseArray stdlib.
I can’t really help with the memory issue. I’ve also observed that memory will sometimes grow indefinitely in long-running code with lots of sparse arrays, so I’ve been suspecting that the wrapped memory underlying
SparseMatrixCSC isn’t always properly finalized, but I haven’t made any attempt to find the bug. This issue might be relevant: Memory leak in solving ODE many times · Issue #1946 · SciML/OrdinaryDiffEq.jl · GitHub.
Is your Julia program parallel, i.e. with two tasks (via @spawn or @threads)? Julia’s gc is not very good at managing memory in such cases. Actually, the only times I’ve had problems with memory usage growing without bounds have been with allocation in multiple threads.
yes, the code is parallel and uses
@threads in a
for loop for assembling the global stiffness matrix from the elemental stiffness matrices, but no memory is allocated within the
for loop, all the required arrays are allocated in advance. I will try to remove that and report what happens.
what I meant by re-using is that the same factorization is used for a few different problems before a new stiffness matrix is assembled and factorized. That is because the cost for factorization is the highest and it makes sense to keep using the same factorized matrix at the cost of taking more iterations rather than assembling and factorizing at each step (that would take one iteration only)
reusing the factorization in the sense you showed in your post proved not to be effective in my case
I’m facing the same problem, i.e. memory usage growing without bounds and I need to use multiple cores with allocations. How did you solve the problem in that case? @Andrea_Vigliotti did you test what happened when you removed parallelization?
In my case, parallelization is done via
pmap. (and no
if Sys.free_memory()/2^30 < 6.0 GC.gc() end
to large functions that are called by pmap…
I think there is a bug in the garbage collector, this should not be needed…
Created a bug report: Garbage collector not working · Issue #50658 · JuliaLang/julia · GitHub
that very much helped! Thanks a lot