Fix OpenMP memory leak and segfault in fragment Hessian diagonalization #1366
+37 −22
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
Replaced OpenMP reduction on large arrays with manual packing/unpacking and critical section updates. Removed large automatic array 'hess_mask' that caused stack overflow. Signed-off-by: Daniel Wiczew <daniel.wiczew@univ-lorraine.fr>
Contributor Author
| Okay, now DCO and CI went through |
foxtran suggested changes Nov 25, 2025
Comment on lines +381 to +396
| do i = 1, nspinsyst(isystem) | ||
| at_i = ispinsyst(i, isystem) | ||
| ! Eigenvalues | ||
| do k = 1, 3 | ||
| eig_calc(3*(at_i-1)+k) = eig(3*(i-1)+k) | ||
| end do | ||
| ! Eigenvectors | ||
| do j = 1, nspinsyst(isystem) | ||
| at_j = ispinsyst(j, isystem) | ||
| do k = 1, 3 | ||
| do l = 1, 3 | ||
| ev_calc(3*(at_i-1)+k, 3*(at_j-1)+l) = mini_hess(3*(i-1)+k, 3*(j-1)+l) | ||
| end do | ||
| end do | ||
| end do | ||
| end do |
Contributor
There was a problem hiding this comment.
Please, split these into two cycles: one for eigval, another for eigvec.
Order of update for eigvec should be:
do j do l do i do k for continuous access for mini_hess
Comment on lines +362 to +372
| do i = 1, nspinsyst(isystem) | ||
| at_i = ispinsyst(i, isystem) | ||
| do j = 1, nspinsyst(isystem) | ||
| at_j = ispinsyst(j, isystem) | ||
| do k = 1, 3 | ||
| do l = 1, 3 | ||
| mini_hess(3*(i-1)+k, 3*(j-1)+l) = hess(3*(at_i-1)+k, 3*(at_j-1)+l) | ||
| end do | ||
| end do | ||
| end do | ||
| end do |
Contributor
There was a problem hiding this comment.
Please, fill mini_hess continuously
Signed-off-by: Daniel Wiczew <daniel.wiczew@univ-lorraine.fr>
foxtran approved these changes Nov 26, 2025
Signed-off-by: Daniel Wiczew <daniel.wiczew@univ-lorraine.fr>
Contributor Author
| I think last my commit is unnecessary |
Member
| No, your last commit (reordering the indices in the hessian as @foxtran suggested), was right. In Fortran, one wants to order the loops such that the fastest changing index (so the innermost loop) addresses the array in its first index (this lies continuously in memory). So please revert this that you have again j, l, i, k as the order of the loops. |
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit. This suggestion is invalid because no changes were made to the code. Suggestions cannot be applied while the pull request is closed. Suggestions cannot be applied while viewing a subset of changes. Only one suggestion per line can be applied in a batch. Add this suggestion to a batch that can be applied as a single commit. Applying suggestions on deleted lines is not supported. You must change the existing code in this line in order to create a valid suggestion. Outdated suggestions cannot be applied. This suggestion has been applied or marked resolved. Suggestions cannot be applied from pending reviews. Suggestions cannot be applied on multi-line comments. Suggestions cannot be applied while the pull request is queued to merge. Suggestion cannot be applied right now. Please check back later.
This commit resolves a segmentation fault and massive memory consumption in the GFN-FF fragment Hessian diagonalization (
frag_hess.f90).The original implementation used a large automatic array
hess_mask(size 3N x 3N) inside an OpenMP parallel region, leading to stack overflows on systems with moderate atom counts (~2000+). Additionally, an OpenMPreduction(+:ev_calc)on the eigenvector matrix caused each thread to allocate a private copy of the full matrix, resulting in O(N_threads * N_atoms^2) memory usage and OOM crashes.Changes:
hess_maskautomatic array.!$omp critical.pack/unpack) with explicit loops to handle fragment data, improving efficiency and reducing memory overhead.Outcomes:
Less technical story:
Generally I had a problem with the frequency calculations when I did use --bhess with 3k+ atom systems. Setting OMP_STACKSIZE and unlimit wasn't helping (the program did run until it ran out of the memory).
So I did some debugging why by compiling with
-g -traceback -check boundsand ran through gdb to get:Things to review: What about the
!$omp do schedule(dynamic)- better to use here static or dynamic ?