Possible troubles when mul!(Y, A , transpose(B) )?

New to Julia
1

Hi, everyone.

I’m new to this language, I’ve been loving it so far, despite having encountered some setbacks. One of them could be a memory leak when multiplying a matrix with its transpose. I use the mul!(Y,A,B) function, and the manual makes it clear that Y should not be related to either A or B:

Blockquote mul!(Y, A, B) → Y. Note that Y must not be aliased with either A or B

However, I have experienced nasty memory leaks (i.e. the program gave correct numerical output, but its CPU usage kept increasing) when A and B are related, as e.g.:

B =transpose(A)
mul!( Y , A , B)

My solution has been to use the transpose! function, as in

transpose(!B , A)
mul!( Y , A , B)

Does this make sense? Unfortunately, I am not able to produce a minimal working code, since the leaks go away when simplifying the full code. The later, I am also unable to provide due to confidentiality reasons! :face_with_diagonal_mouth:

[ BTW, multiplying a matrix by its transpose, thus producing a positive semidefinite symmetric matrix, M =A’ A, which later appears in a linear system M x = y, to be solved as x = M \ y, is an extremely common situation. I appreciate general comments about this, but perhaps it merits a separate topic… ]

2

Is it both a memory leak and CPU going out of control, or was CPU supposed to say memory? If CPU, does it keep increasing even after the computation, or does it stop then?

Can you verify that the leak happens when you have mul!(Y, A, transpose(A)) but not for mul!(Y, A, copy(transpose(A))) in the original program, without changing other things?

3

In general, for dense A, mul!(Y, A, transpose(A)) should call BLAS.syrk!, which shouldn’t do anything very surprising.

1 Like
4

Sorry for the confusion: I meant the RAM kept increasing, as seen e.g. in “free” or “top” in linux.

Will try the copy() thing, but, as you know (and I learned the hard way, but that’s another story) A=copy(B) make still link the two matrices, as transpose() does.

I attach the piece of code that made this clear for me
Sin título

5

Well, I don’t think it links the matrix A and B in any way, but if B contains references it will only copy the reference, thus still referring to the same memory.

What you create in the example code is not a matrix, but a vector of vectors. In julia these are two very different things, where the outer vector contains references to the inner vector. Thus a copy will create a new vector with the same references, and updating the values in the inner vectors will be reflected in both a and b since they point to the same memory.

6

Yes, you are completely right. copy() should not link two matrices in any way. It was just my bad luck that I happened to define a vector of vectors in my code and run into some problems (it is a particle simulation, so each particle has vectors assigned to it, its velocity for example).

7

If you want to preserve that vector-of-vectors structure, I would recommend switching the “inner” vectors to static vectors from StaticArrays.jl. Prerrquisite is them alle being the same size (which should be the case for particle simulations).

8

I can’t reproduce this behavior. As shown below, there are no allocations from mul!. An MWE would help to figure out what the underlying problem is.

julia> A = [1 2; 3 4];

julia> B = transpose(A);

julia> Y = similar(A);

julia> @btime mul!(Y, A, B)
  61.312 ns (0 allocations: 0 bytes)
2×2 Matrix{Int64}:
  5  11
 11  25
2 Likes
11

Good suggestion, in fact I have been using static vectors for the positions, but not for velocities or other vectors (forces, accelerations…). Will definitely do that.

12

I know but, as I said, I have been unable to produce an MWE. Trimming the program makes the problem go away somehow, and the full program is not to be made public yet. Perhaps in a few months, when the research is published.