Memory leak for linear scaling ompGemm_m kernel on atom movement

[davidbowler]

I have identified a significant memory leak in O(N) geometry optimisation, though not the source. The system (which I will attach input files for when I have a better internet connection) is 4,096 atoms of silicon (8x8x8 fundamental cells) running on 64 processes (though this does not seem to matter) with a SZ basis and L range of 16 bohr. With the develop branch, at some point after the initial movement of atoms the memory use spikes (I've seen at least 16GB in one process) and even with 1.5TB of RAM the job fails. The same behaviour is not seen for v1.2.

[davidbowler]

Input files attached.
MemoryError.zip

[davidbowler]

The problem seems to lie only in the ompGemm_m version of multiply_kernel. After atoms are moved, the maxima found at the start of m_kern_max and m_kern_min become very large (and in some cases overflow) leading to the memory problem. I'm not yet sure why these maxima fail...

[tkoskela]

I created a new test in the testsuite test_004 using the files in MemoryError.zip. It seems like all the multiply kernels fail on it.

See https://github.com/OrderN/CONQUEST-release/actions/runs/7035615252/job/19146415009

Interesting enough, ompGemm_m fails in test_003, so maybe there is a different problem there

[tkoskela]

Regarding ompGemm_m, running test_003 on my laptop it looks like there can be uninitialized values in bndim2 and nbnab. These come from pointers to the prefetched array in multiply_module.

The problem is that we are assuming a size mx_part for the integer arrays nbnab,ibarddr, etc. but this is not always correct. Because we are using a pointer, we don't explicitly allocate a size for the array (although it looks like it might be possible to figure out from the code in the kernel). The maxval function then looks at all elements from 1 to mx_part and in case of some of the pointers finds uninitialized memory.

Personally I would be strongly in favour of getting rid of the pointers. One possible solution would be to just pass the integer indices to each sub-array instead, if the underlying data structure is too difficult to change.

[tkoskela]

A simple solution to this particular issue seems to be using assumed-shape arrays in the multiply kernel, ie. integer(integ) :: ibaddr(:) instead of integer(integ) :: ibaddr(mx_part). How exactly the correct shape gets propagated to the kernel is anyone's guess, maybe because the pointer is pointing to a chunk of memory of specified size, that information is contained in the pointer.