From 4b578d9849df7fa548a0b8627f714df6248a8440 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 20 Jun 2024 13:43:46 +0200
Subject: [PATCH] mmap is now shared in cholesky

---
 src/ao_two_e_ints/cholesky.irp.f | 36 ++++++++++++++++++--------------
 src/utils/fortran_mmap.c         |  3 +++
 2 files changed, 23 insertions(+), 16 deletions(-)

diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f
index a680e7ee..063cc365 100644
--- a/src/ao_two_e_ints/cholesky.irp.f
+++ b/src/ao_two_e_ints/cholesky.irp.f
@@ -31,12 +31,11 @@ END_PROVIDER
    integer*8                      :: ndim8
    integer                        :: rank
    double precision               :: tau, tau2
-   double precision, pointer      :: L(:,:), Delta(:,:)
+   double precision, pointer      :: L(:,:)
 
    double precision               :: s
-   double precision               :: dscale, dscale_tmp
 
-   double precision, allocatable  :: D(:), Ltmp_p(:,:), Ltmp_q(:,:), D_sorted(:), Delta_col(:)
+   double precision, allocatable  :: D(:), Ltmp_p(:,:), Ltmp_q(:,:), D_sorted(:), Delta_col(:), Delta(:,:)
    integer, allocatable           :: addr1(:), addr2(:)
    integer*8, allocatable         :: Lset(:), Dset(:)
    logical, allocatable           :: computed(:)
@@ -102,7 +101,7 @@ END_PROVIDER
      rank = 0
 
      allocate( D(ndim8), Lset(ndim8), Dset(ndim8), D_sorted(ndim8))
-     allocate( addr1(ndim8), addr2(ndim8), Delta_col(ndim8) )
+     allocate( addr1(ndim8), addr2(ndim8), Delta_col(ndim8), computed(ndim8) )
 
      call resident_memory(mem0)
 
@@ -149,11 +148,9 @@ END_PROVIDER
      Dmax = D_sorted(1)
 
      ! 2.
-     dscale = 1.d0
-     dscale_tmp = dscale*dscale*Dmax
      np8=0_8
      do p8=1,ndim8
-       if ( dscale_tmp*D(p8) >= tau2 ) then
+       if ( Dmax*D(p8) >= tau2 ) then
          np8 = np8+1_8
          Lset(np8) = p8
        endif
@@ -203,16 +200,23 @@ END_PROVIDER
 
          mem = mem0                                 &
              + np*memory_of_double(nq)              & ! Delta(np,nq)
-             + (np+nq)*memory_of_double(block_size) & ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
-             + memory_of_int(nq)                      ! computed(nq)
+             + (np+nq)*memory_of_double(block_size)   ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
 
          if (mem > qp_max_mem*0.5d0) then
-           nq = nq/2
+           Dmin = D_sorted(nq/2)
+           do ii=nq/2,np-1
+             if (D_sorted(ii) < Dmin) then
+               nq = ii
+               exit
+             endif
+           enddo
          else
            exit
          endif
 
        enddo
+!call print_memory_usage
+!print *, 'np, nq, Predicted memory: ', np, nq, mem
 
        if (nq <= 0) then
          print *, nq
@@ -247,8 +251,7 @@ END_PROVIDER
        endif
 
 
-       allocate(computed(nq))
-       computed(:) = .False.
+       computed(1:nq) = .False.
 
 
        !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q)
@@ -406,7 +409,6 @@ END_PROVIDER
 
        deallocate(Ltmp_p)
        deallocate(Ltmp_q)
-       deallocate(computed)
        deallocate(Delta)
 
        ! i.
@@ -419,11 +421,9 @@ END_PROVIDER
 
        Dmax = D_sorted(1)
 
-       dscale = 1.d0
-       dscale_tmp = dscale*dscale*Dmax
        np8=0_8
        do p8=1,ndim8
-         if ( dscale_tmp*D(p8) >= tau2 ) then
+         if ( Dmax*D(p8) >= tau2 ) then
            np8 = np8+1_8
            Lset(np8) = p8
          endif
@@ -436,6 +436,10 @@ END_PROVIDER
      print *,  '============ ============='
      print *,  ''
 
+     deallocate( D, Lset, Dset, D_sorted )
+     deallocate( addr1, addr2, Delta_col, computed )
+
+
      allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
 
      if (ierr /= 0) then
diff --git a/src/utils/fortran_mmap.c b/src/utils/fortran_mmap.c
index 2dbe42b8..0306f64f 100644
--- a/src/utils/fortran_mmap.c
+++ b/src/utils/fortran_mmap.c
@@ -49,10 +49,13 @@ void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only,
         }
 
         if (single_node == 1) {
+          map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
+/*
           map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_POPULATE | MAP_NONBLOCK | MAP_NORESERVE, fd, 0);
           if (map == MAP_FAILED) {
              map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_PRIVATE, fd, 0);
           }
+*/
         } else {
           map = mmap(NULL, bytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
         }