mmap works for integrals

plugins/Mmap/fortran_mmap.c

@@ -0,0 +1,72 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include <sys/mman.h>
+
+
+void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only)
+{
+    int i;
+    int fd;
+    int result;
+    void* map;
+
+    if (read_only == 1)
+    {
+        fd = open(filename, O_RDONLY, (mode_t)0600);
+        if (fd == -1) {
+            printf("%s:\n", filename);
+            perror("Error opening mmap file for reading");
+            exit(EXIT_FAILURE);
+        }
+        map = mmap(0, bytes, PROT_READ, MAP_SHARED, fd, 0);
+    }
+    else
+    {
+        fd = open(filename, O_RDWR | O_CREAT, (mode_t)0600);
+        if (fd == -1) {
+            printf("%s:\n", filename);
+            perror("Error opening mmap file for writing");
+            exit(EXIT_FAILURE);
+        }
+
+        result = lseek(fd, bytes, SEEK_SET);
+        if (result == -1) {
+            close(fd);
+            printf("%s:\n", filename);
+            perror("Error calling lseek() to stretch the file");
+            exit(EXIT_FAILURE);
+        }
+        
+        result = write(fd, "", 1);
+        if (result != 1) {
+            close(fd);
+            printf("%s:\n", filename);
+            perror("Error writing last byte of the file");
+            exit(EXIT_FAILURE);
+        }
+
+        map = mmap(0, bytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
+    }
+
+    if (map == MAP_FAILED) {
+        close(fd);
+        printf("%s:\n", filename);
+        perror("Error mmapping the file");
+        exit(EXIT_FAILURE);
+    }
+
+    *file_descr = fd;
+    return map;
+}
+
+void munmap_fortran(size_t bytes, int fd, void* map)
+{
+    if (munmap(map, bytes) == -1) {
+        perror("Error un-mmapping the file");
+    }
+    close(fd);
+}

plugins/Mmap/mmap.f90

@@ -0,0 +1,69 @@
+module mmap_module
+
+  use iso_c_binding
+
+  interface
+
+    ! File descriptors
+    ! ----------------
+
+    type(c_ptr) function c_mmap_fortran(filename, length, fd, read_only) bind(c,name='mmap_fortran')
+      use iso_c_binding
+      character(c_char), intent(in)  :: filename(*)
+      integer(c_size_t), intent(in), value :: length
+      integer(c_int), intent(out)    :: fd
+      integer(c_int), intent(in), value    :: read_only
+    end function
+
+    subroutine c_munmap(length, fd, map) bind(c,name='munmap_fortran')
+      use iso_c_binding
+      integer(c_size_t), intent(in), value :: length
+      integer(c_int), intent(in), value :: fd
+      type(c_ptr), intent(in), value    :: map
+    end subroutine
+
+  end interface
+
+  contains
+
+  subroutine mmap(filename, shape, bytes, fd, read_only, map)
+      use iso_c_binding
+      implicit none
+      character*(*), intent(in)      :: filename   ! Name of the mapped file
+      integer*8, intent(in)          :: shape(:)   ! Shape of the array to map
+      integer, intent(in)            :: bytes      ! Number of bytes per element
+      logical, intent(in)            :: read_only  ! If true, mmap is read-only
+      integer, intent(out)           :: fd         ! File descriptor
+      type(c_ptr), intent(out)       :: map        ! C Pointer
+     
+      integer(c_long)                :: length
+      integer(c_int)                 :: fd_
+
+      length = PRODUCT( shape(:) ) * bytes
+      if (read_only) then
+          map = c_mmap_fortran( trim(filename)//char(0), length, fd_, 1)
+      else
+          map = c_mmap_fortran( trim(filename)//char(0), length, fd_, 0)
+      endif
+      fd = fd_
+  end subroutine
+
+  subroutine munmap(shape, bytes, fd, map)
+      use iso_c_binding
+      implicit none
+      integer*8, intent(in)          :: shape(:)  ! Shape of the array to map
+      integer, intent(in)            :: bytes     ! Number of bytes per element
+      integer, intent(in)            :: fd        ! File descriptor
+      type(c_ptr), intent(in)        :: map       ! C pointer
+
+      integer(c_long)                :: length
+      integer(c_int)                 :: fd_
+
+      length = PRODUCT( shape(:) ) * bytes
+      fd_ = fd
+      call c_munmap( length, fd_, map)
+  end 
+
+end module mmap_module
+
+

src/Integrals_Bielec/NEEDED_CHILDREN_MODULES

@@ -1 +1,2 @@
-Pseudo Bitmask ZMQ
+Pseudo Bitmask ZMQ Mmap
+

src/Integrals_Bielec/ao_bi_integrals.irp.f

@@ -353,11 +353,10 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   if (read_ao_integrals) then
     integer                        :: load_ao_integrals
     print*,'Reading the AO integrals'
-    if (load_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') == 0) then
+      call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
       print*, 'AO integrals provided'
       ao_bielec_integrals_in_map = .True.
       return
-    endif
   endif
   
   print*, 'Providing the AO integrals'
@@ -402,8 +401,10 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   print*, ' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
   
   ao_bielec_integrals_in_map = .True.
+
   if (write_ao_integrals) then
-    call dump_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin')
+    call ezfio_set_work_empty(.False.)
+    call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
     call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
   endif
   

src/Integrals_Bielec/mo_bi_integrals.irp.f

@@ -30,10 +30,9 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
   if (read_mo_integrals) then
     integer                        :: load_mo_integrals
     print*,'Reading the MO integrals'
-    if (load_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') == 0) then
-      print*, 'MO integrals provided'
-      return
-    endif
+    call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
+    print*, 'MO integrals provided'
+    return
   endif
   
   call add_integrals_to_map(full_ijkl_bitmask_4)
@@ -299,7 +298,8 @@ subroutine add_integrals_to_map(mask_ijkl)
   print*,' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
   
   if (write_mo_integrals) then
-    call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin')
+    call ezfio_set_work_empty(.False.)
+    call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
     call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
   endif
   

src/Utils/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-
+Mmap

src/Utils/map_module.f90

@@ -30,8 +30,8 @@ module map_module
  integer*8, parameter           :: map_mask  = ibset(0_8,15)-1_8
 
  type cache_map_type
-  integer(cache_key_kind), pointer :: key(:)
   real(integral_kind), pointer   :: value(:)
+  integer(cache_key_kind), pointer :: key(:)
   logical                        :: sorted
   integer(cache_map_size_kind)   :: map_size
   integer(cache_map_size_kind)   :: n_elements
@@ -39,10 +39,14 @@ module map_module
  end type cache_map_type
 
  type map_type
-  type(cache_map_type), allocatable :: map(:)
+  type(cache_map_type), pointer  :: map(:)
+  real(integral_kind), pointer   :: consolidated_value(:)
+  integer(cache_key_kind), pointer :: consolidated_key(:)
+  integer*8, pointer             :: consolidated_idx(:)
+  logical                        :: sorted
+  logical                        :: consolidated
   integer(map_size_kind)         :: map_size
   integer(map_size_kind)         :: n_elements
-  logical                        :: sorted
   integer(omp_lock_kind)         :: lock
  end type map_type
 
@@ -92,6 +96,7 @@ subroutine map_init(map,keymax)
   
   map%n_elements = 0_8
   map%map_size = ishft(keymax,map_shift)
+  map%consolidated = .False.
   
   allocate(map%map(0_8:map%map_size),stat=err)
   if (err /= 0) then
@@ -655,12 +660,10 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
     idx = ibegin
     if (min(iend_in,sze) > ibegin+16) then
       iend = ibegin+16
-      !DIR$ VECTOR ALIGNED
       do while (cache_key > X(idx))
         idx = idx+1
       end do
     else
-      !DIR$ VECTOR ALIGNED
       do while (cache_key > X(idx))
         idx = idx+1
         if (idx /= iend) then
@@ -768,13 +771,11 @@ subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in
     value = Y(idx)
     if (min(iend_in,sze) > ibegin+16) then
       iend = ibegin+16
-      !DIR$ VECTOR ALIGNED
       do while (cache_key > X(idx))
         idx = idx+1
         value = Y(idx)
       end do
     else
-      !DIR$ VECTOR ALIGNED
       do while (cache_key > X(idx))
         idx = idx+1
         value = Y(idx)
@@ -853,3 +854,121 @@ subroutine get_cache_map(map,map_idx,keys,values,n_elements)
   enddo
   
 end
+
+
+subroutine map_save_to_disk(filename,map)
+  use map_module
+  use mmap_module
+  implicit none
+  character*(*), intent(in)      :: filename
+  type(map_type), intent(inout)  :: map
+  type(c_ptr)                    :: c_pointer(3)
+  integer                        :: fd(3)
+  integer*8                      :: i,k
+  integer                        :: j
+
+
+
+  if (map % consolidated) then
+    stop 'map already consolidated'
+  endif
+
+  call mmap(trim(filename)//'_consolidated_idx', (/ map % map_size + 2_8 /), 8, fd(1), .False., c_pointer(1))
+  call c_f_pointer(c_pointer(1),map % consolidated_idx, (/ map % map_size +2_8/))
+
+  call mmap(trim(filename)//'_consolidated_key', (/ map % n_elements /), cache_key_kind, fd(2), .False., c_pointer(2))
+  call c_f_pointer(c_pointer(2),map % consolidated_key, (/ map % n_elements /))
+
+  call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .False., c_pointer(3))
+  call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
+
+  if (.not.associated(map%consolidated_key)) then
+    stop 'cannot consolidate map : consolidated_key not associated'
+  endif
+
+  if (.not.associated(map%consolidated_value)) then
+    stop 'cannot consolidate map : consolidated_value not associated'
+  endif
+
+  if (.not.associated(map%consolidated_idx)) then
+    stop 'cannot consolidate map : consolidated_idx not associated'
+  endif
+
+  call map_sort(map)
+  k = 1_8
+  do i=0_8, map % map_size
+    map % consolidated_idx (i+1) = k
+    do j=1, map % map(i) % n_elements
+      map % consolidated_value(k) = map % map(i) % value(j)
+      map % consolidated_key  (k) = map % map(i) % key(j)
+      k = k+1_8
+    enddo
+    deallocate(map % map(i) % value)
+    deallocate(map % map(i) % key)
+    map % map(i) % value => map % consolidated_value ( map % consolidated_idx (i+1) :)
+    map % map(i) % key   => map % consolidated_key   ( map % consolidated_idx (i+1) :)
+  enddo
+  map % consolidated_idx (map % map_size + 2_8) = k
+  map % consolidated = .True.
+
+
+  call munmap( (/ map % map_size + 2_8 /), 8, fd(1), c_pointer(1))
+  call mmap(trim(filename)//'_consolidated_idx', (/ map % map_size + 2_8 /), 8, fd(1), .True., c_pointer(1))
+  call c_f_pointer(c_pointer(1),map % consolidated_idx, (/ map % map_size +2_8/))
+
+  call munmap( (/ map % n_elements /), cache_key_kind, fd(2), c_pointer(2))
+  call mmap(trim(filename)//'_consolidated_key', (/ map % n_elements /), cache_key_kind, fd(2), .True., c_pointer(2))
+  call c_f_pointer(c_pointer(2),map % consolidated_key, (/ map % n_elements /))
+
+  call munmap( (/ map % n_elements /), integral_kind, fd(3), c_pointer(3))
+  call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .True., c_pointer(3))
+  call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
+
+end
+
+subroutine map_load_from_disk(filename,map)
+  use map_module
+  use mmap_module
+  implicit none
+  character*(*), intent(in)      :: filename
+  type(map_type), intent(inout)  :: map
+  type(c_ptr)                    :: c_pointer(3)
+  integer                        :: fd(3)
+  integer*8                      :: i,k
+  integer                        :: n_elements
+
+
+
+  if (map % consolidated) then
+    stop 'map already consolidated'
+  endif
+
+  call mmap(trim(filename)//'_consolidated_idx', (/ map % map_size + 2_8 /), 8, fd(1), .True., c_pointer(1))
+  call c_f_pointer(c_pointer(1),map % consolidated_idx, (/ map % map_size + 2_8/))
+
+  map% n_elements = map % consolidated_idx (map % map_size+2_8)-1
+
+  call mmap(trim(filename)//'_consolidated_key', (/ map % n_elements /), cache_key_kind, fd(2), .True., c_pointer(2))
+  call c_f_pointer(c_pointer(2),map % consolidated_key, (/ map % n_elements /))
+
+  call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .True., c_pointer(3))
+  call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
+
+  k = 1_8
+  do i=0_8, map % map_size
+    deallocate(map % map(i) % value)
+    deallocate(map % map(i) % key)
+    map % map(i) % value => map % consolidated_value ( map % consolidated_idx (i+1) :)
+    map % map(i) % key   => map % consolidated_key   ( map % consolidated_idx (i+1) :)
+    map % map(i) % sorted = .True.
+    n_elements = map % consolidated_idx (i+2) - k
+    k = map % consolidated_idx (i+2)
+    map % map(i) % map_size = n_elements
+    map % map(i) % n_elements = n_elements
+  enddo
+  map % n_elements = k-1
+  map % sorted = .True. 
+  map % consolidated = .True.
+
+end
+