10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-05 05:33:56 +01:00
quantum_package/src/Utils/map_module.f90
2017-04-12 20:23:04 +02:00

860 lines
23 KiB
Fortran

module map_module
! A map is an array of maps (cache_maps)
! A cache map is an array of keys and values sorted by keys
! A cache map has its own OpenMP lock
! To access a (key,value) pair in the map, the
! index of the cache_map in the map array is obtained
! by removing the first 15 bits of the key.
! The key in the cache_map is composed of the first
! 15 bits of the key. Therefore, it can be stored
! as integer*2 and is found by applying the map_mask
! to the initial key. The element are found in the
! cache_map using a binary search
!
! When using the map_update subroutine to build the map,
! the map_unique subroutine
! should be called before getting data from the map.
use omp_lib
integer, parameter :: integral_kind = 8
integer, parameter :: cache_key_kind = 2
integer, parameter :: cache_map_size_kind = 4
integer, parameter :: key_kind = 8
integer, parameter :: map_size_kind = 8
integer, parameter :: map_shift = -15
integer*8, parameter :: map_mask = ibset(0_8,15)-1_8
type cache_map_type
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
integer(omp_lock_kind) :: lock
end type cache_map_type
type map_type
type(cache_map_type), allocatable :: 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
integer(omp_lock_kind) :: lock
end type map_type
end module map_module
double precision function map_mb(map)
use map_module
use omp_lib
implicit none
type (map_type), intent(in) :: map
integer(map_size_kind) :: i
map_mb = dble(8+map_size_kind+map_size_kind+omp_lock_kind+4)
do i=0,map%map_size
map_mb = map_mb + dble(map%map(i)%map_size*(cache_key_kind+integral_kind) +&
8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind)
enddo
map_mb = map_mb / (1024.d0*1024.d0)
end
subroutine cache_map_init(map,sze)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_map_size_kind) :: sze
call omp_set_lock(map%lock)
map%n_elements = 0_8
map%map_size = 0_8
map%sorted = .True.
NULLIFY(map%value, map%key)
call cache_map_reallocate(map,sze)
call omp_unset_lock(map%lock)
end
subroutine map_init(map,keymax)
use map_module
implicit none
integer*8, intent(in) :: keymax
type (map_type), intent(inout) :: map
integer(map_size_kind) :: i
integer(cache_map_size_kind) :: sze
integer :: err
call omp_init_lock(map%lock)
call omp_set_lock(map%lock)
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
print *, 'Unable to allocate map'
stop 5
endif
sze = 2
do i=0_8,map%map_size
call omp_init_lock(map%map(i)%lock)
enddo
!$OMP PARALLEL DEFAULT(NONE) SHARED(map,sze) PRIVATE(i)
!$OMP DO SCHEDULE(STATIC,512)
do i=0_8,map%map_size
call cache_map_init(map%map(i),sze)
enddo
!$OMP ENDDO
!$OMP END PARALLEL
map%sorted = .True.
call omp_unset_lock(map%lock)
end
subroutine cache_map_reallocate(map,sze)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
type (cache_map_type), intent(inout) :: map
integer(cache_key_kind), pointer :: key_new(:)
real(integral_kind), pointer :: value_new(:)
integer(map_size_kind) :: i
integer :: err
!DIR$ ATTRIBUTES ALIGN : 64 :: key_new, value_new
if (sze < map%n_elements) then
print *, 'Unable to resize map : map too large'
stop 3
endif
! Resize keys
allocate( key_new(sze), stat=err )
if (err /= 0) then
print *, 'Unable to allocate map', sze
stop 1
endif
if (associated(map%key)) then
do i=1_8,min(size(map%key),map%n_elements)
key_new(i) = map%key(i)
enddo
deallocate(map%key)
endif
! Resize values
allocate( value_new(sze), stat=err )
if (err /= 0) then
print *, 'Unable to allocate map', sze
stop 2
endif
if (associated(map%value)) then
do i=1_8,min(size(map%key),map%n_elements)
value_new(i) = map%value(i)
enddo
deallocate(map%value)
endif
! Set new pointers
map%key => key_new
map%value => value_new
map%map_size = sze
end
subroutine cache_map_deinit(map)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer :: err
if (associated( map % value )) then
deallocate( map % value, stat=err )
if (err /= 0) then
print *, 'Unable to deallocate map'
stop 2
endif
NULLIFY(map%value)
endif
if (associated( map % key )) then
deallocate( map % key, stat=err )
if (err /= 0) then
print *, 'Unable to deallocate map'
stop 4
endif
NULLIFY(map%key)
endif
map%n_elements = 0_8
map%map_size = 0_8
call omp_destroy_lock(map%lock)
end
subroutine map_deinit(map)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer :: err
integer(map_size_kind) :: i
if (allocated( map % map )) then
do i=0_8, map%map_size
call cache_map_deinit(map%map(i))
enddo
deallocate( map % map, stat=err )
if (err /= 0) then
print *, 'Unable to deallocate map'
stop 6
endif
endif
map%n_elements = 0_8
map%map_size = 0_8
call omp_destroy_lock(map%lock)
end
subroutine cache_map_sort(map)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_map_size_kind), allocatable :: iorder(:)
integer(cache_map_size_kind) :: i
!DIR$ ATTRIBUTES ALIGN : 64 :: iorder
if (.not.map%sorted) then
allocate(iorder(map%n_elements))
do i=1,map%n_elements
iorder(i) = i
enddo
if (cache_key_kind == 2) then
call i2radix_sort(map%key,iorder,map%n_elements,-1)
else if (cache_key_kind == 4) then
call iradix_sort(map%key,iorder,map%n_elements,-1)
else if (cache_key_kind == 8) then
call i8radix_sort(map%key,iorder,map%n_elements,-1)
endif
if (integral_kind == 4) then
call set_order(map%value,iorder,map%n_elements)
else if (integral_kind == 8) then
call dset_order(map%value,iorder,map%n_elements)
endif
deallocate(iorder)
map%sorted = .True.
endif
end
subroutine map_sort(map)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer(map_size_kind) :: i
if (.not.map%sorted) then
!$OMP PARALLEL DO SCHEDULE(static,1024) DEFAULT(SHARED) PRIVATE(i)
do i=0_8,map%map_size
call omp_set_lock(map%map(i)%lock)
call cache_map_sort(map%map(i))
call omp_unset_lock(map%map(i)%lock)
enddo
!$OMP END PARALLEL DO
map%sorted = .True.
endif
end
subroutine cache_map_unique(map)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_key_kind) :: prev_key
integer(cache_map_size_kind) :: i, j
call cache_map_sort(map)
prev_key = -1_8
j=0
do i=1,map%n_elements
if (map%key(i) /= prev_key) then
j = j+1
map%value(j) = map%value(i)
map%key(j) = map%key(i)
prev_key = map%key(i)
else
map%value(j) = map%value(j)+map%value(i)
endif
enddo
map%n_elements = j
end
subroutine cache_map_shrink(map,thr)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
real(integral_kind) , intent(in) :: thr
integer(cache_map_size_kind) :: i,j
j=0
do i=1,map%n_elements
if (abs(map%value(i)) > thr) then
j = j+1
map%value(j) = map%value(i)
map%key(j) = map%key(i)
endif
enddo
map%n_elements = j
end
subroutine map_unique(map)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer(map_size_kind) :: i
integer(map_size_kind) :: icount
icount = 0_8
!$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
!$OMP REDUCTION(+:icount)
do i=0_8,map%map_size
call omp_set_lock(map%map(i)%lock)
call cache_map_unique(map%map(i))
call omp_unset_lock(map%map(i)%lock)
icount = icount + map%map(i)%n_elements
enddo
!$OMP END PARALLEL DO
map%n_elements = icount
end
subroutine map_shrink(map,thr)
use map_module
implicit none
type (map_type), intent(inout) :: map
real(integral_kind), intent(in) :: thr
integer(map_size_kind) :: i
integer(map_size_kind) :: icount
icount = 0_8
!$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
!$OMP REDUCTION(+:icount)
do i=0_8,map%map_size
call omp_set_lock(map%map(i)%lock)
call cache_map_shrink(map%map(i),thr)
call omp_unset_lock(map%map(i)%lock)
icount = icount + map%map(i)%n_elements
enddo
!$OMP END PARALLEL DO
map%n_elements = icount
end
subroutine map_update(map, key, value, sze, thr)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(inout) :: key(sze)
real(integral_kind), intent(inout) :: value(sze)
real(integral_kind), intent(in) :: thr
integer :: i
integer(map_size_kind) :: idx_cache, idx_cache_new
integer(cache_map_size_kind) :: idx
integer :: sze2
integer(cache_key_kind) :: cache_key
integer(map_size_kind) :: n_elements_temp
type (cache_map_type) :: local_map
logical :: map_sorted
sze2 = sze
map_sorted = .True.
n_elements_temp = 0_8
n_elements_temp = n_elements_temp + 1_8
do while (sze2>0)
i=1
do while (i<=sze)
if (key(i) /= 0_8) then
idx_cache = ishft(key(i),map_shift)
if (omp_test_lock(map%map(idx_cache)%lock)) then
local_map%key => map%map(idx_cache)%key
local_map%value => map%map(idx_cache)%value
local_map%sorted = map%map(idx_cache)%sorted
local_map%map_size = map%map(idx_cache)%map_size
local_map%n_elements = map%map(idx_cache)%n_elements
do
!DIR$ FORCEINLINE
call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements)
if (idx > 0_8) then
local_map%value(idx) = local_map%value(idx) + value(i)
else
! Assert that the map has a proper size
if (local_map%n_elements == local_map%map_size) then
call cache_map_unique(local_map)
call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
call cache_map_shrink(local_map,thr)
endif
cache_key = int(iand(key(i),map_mask),2)
local_map%n_elements = local_map%n_elements + 1
local_map%value(local_map%n_elements) = value(i)
local_map%key(local_map%n_elements) = cache_key
local_map%sorted = .False.
n_elements_temp = n_elements_temp + 1_8
endif ! idx > 0
key(i) = 0_8
i = i+1
sze2 = sze2-1
if (i>sze) then
i=1
endif
if ( (ishft(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then
exit
endif
enddo
map%map(idx_cache)%key => local_map%key
map%map(idx_cache)%value => local_map%value
map%map(idx_cache)%sorted = local_map%sorted
map%map(idx_cache)%n_elements = local_map%n_elements
map%map(idx_cache)%map_size = local_map%map_size
map_sorted = map_sorted .and. local_map%sorted
call omp_unset_lock(map%map(idx_cache)%lock)
endif ! omp_test_lock
else
i=i+1
endif ! key = 0
enddo ! i
enddo ! sze2 > 0
call omp_set_lock(map%lock)
map%n_elements = map%n_elements + n_elements_temp
map%sorted = map%sorted .and. map_sorted
call omp_unset_lock(map%lock)
end
subroutine map_append(map, key, value, sze)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(inout) :: key(sze)
real(integral_kind), intent(inout) :: value(sze)
integer :: i
integer(cache_map_size_kind) :: n_elements
integer(map_size_kind) :: idx_cache
integer(cache_key_kind) :: cache_key
do i=1,sze
idx_cache = ishft(key(i),map_shift)
call omp_set_lock(map%map(idx_cache)%lock)
n_elements = map%map(idx_cache)%n_elements + 1
! Assert that the map has a proper size
if (n_elements == map%map(idx_cache)%map_size) then
call cache_map_reallocate(map%map(idx_cache), n_elements+ ishft(n_elements,-1))
endif
cache_key = int(iand(key(i),map_mask),2)
map%map(idx_cache)%value(n_elements) = value(i)
map%map(idx_cache)%key(n_elements) = cache_key
map%map(idx_cache)%n_elements = n_elements
if (map%map(idx_cache)%sorted.and.n_elements > 1) then
map%map(idx_cache)%sorted = (map%map(idx_cache)%key(n_elements-1) <= cache_key)
map%sorted = map%sorted .and. map%map(idx_cache)%sorted
endif
call omp_unset_lock(map%map(idx_cache)%lock)
enddo
call omp_set_lock(map%lock)
map%n_elements = map%n_elements + sze
call omp_unset_lock(map%lock)
end
subroutine map_get(map, key, value)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer(key_kind), intent(in) :: key
real(integral_kind), intent(out) :: value
integer(map_size_kind) :: idx_cache
integer(cache_map_size_kind) :: idx
! index in tha pointers array
idx_cache = ishft(key,map_shift)
!DIR$ FORCEINLINE
call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx)
end
subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(key_kind), intent(in) :: key
integer(cache_map_size_kind), intent(in) :: ibegin, iend
real(integral_kind), intent(out) :: value
integer(cache_map_size_kind), intent(inout) :: idx
double precision, pointer :: v(:)
integer :: i
! call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend)
! if (idx > 0) then
! value = v(idx)
! else
! value = 0._integral_kind
! endif
end
subroutine map_get_many(map, key, value, sze)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(in) :: key(sze)
real(integral_kind), intent(out) :: value(sze)
integer :: i
integer(map_size_kind) :: idx_cache
integer(cache_map_size_kind) :: ibegin, iend
integer(cache_map_size_kind), allocatable :: idx(:)
!DIR$ ATTRIBUTES ALIGN : 64 :: idx
allocate(idx(sze))
do i=1,sze
idx_cache = ishft(key(i),map_shift)
iend = map%map(idx_cache)%n_elements
!DIR$ FORCEINLINE
call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend)
enddo
do i=1,sze
idx_cache = ishft(key(i),map_shift)
if (idx(i) > 0) then
value(i) = map%map(idx_cache)%value(idx(i))
else
value(i) = 0.
endif
enddo
deallocate(idx)
end
subroutine map_exists_many(map, key, sze)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(inout) :: key(sze)
integer :: i
integer(map_size_kind) :: idx_cache, idx_cache_prev
integer(cache_map_size_kind) :: ibegin, iend
integer(cache_map_size_kind), allocatable :: idx(:)
!DIR$ ATTRIBUTES ALIGN : 64 :: idx
idx_cache_prev = -1_map_size_kind
allocate(idx(sze))
do i=1,sze
idx_cache = ishft(key(i),map_shift)
iend = map%map(idx_cache)%n_elements
if (idx_cache == idx_cache_prev) then
if ((idx(i-1) > 0_cache_map_size_kind).and.(idx(i-1) < iend)) then
if ((key(i) == key(i-1)+1).and.(map%map(idx_cache)%key(idx(i-1))+1) == key(i)) then
idx(i) = idx(i-1)+1
cycle
endif
endif
endif
!DIR$ FORCEINLINE
call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend)
idx_cache_prev = idx_cache
enddo
do i=1,sze
idx_cache = ishft(key(i),map_shift)
if (idx(i) <= 0) then
key(i) = 0_key_kind
endif
enddo
deallocate(idx)
end
subroutine search_key_big(key,X,sze,idx)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
integer(key_kind) , intent(in) :: key
integer(cache_key_kind) , intent(in) :: X(sze)
integer(cache_map_size_kind), intent(out) :: idx
call search_key_big_interval(key,X,sze,idx,1,sze)
end
subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
integer(key_kind) , intent(in) :: key
integer(cache_key_kind) , intent(in) :: X(sze)
integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in
integer(cache_map_size_kind), intent(out) :: idx
integer(cache_map_size_kind) :: istep, ibegin, iend, i
integer(cache_key_kind) :: cache_key
if (sze /= 0) then
continue
else
idx = -1
return
endif
cache_key = int(iand(key,map_mask),2)
ibegin = min(ibegin_in,sze)
iend = min(iend_in,sze)
if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
istep = ishft(iend-ibegin,-1)
idx = ibegin + istep
do while (istep > 64)
idx = ibegin + istep
! TODO : Cache misses
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
idx = ibegin + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
cycle
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
cycle
else
return
endif
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
idx = idx + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
cycle
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
cycle
else
return
endif
else
return
endif
enddo
idx = ibegin
if (min(iend_in,sze) > ibegin+64) then
iend = ibegin+64
do while (cache_key > X(idx))
idx = idx+1
end do
else
do while (cache_key > X(idx))
idx = idx+1
if (idx /= iend) then
cycle
else
exit
endif
end do
endif
if (cache_key /= X(idx)) then
idx = 1-idx
endif
return
else
if (cache_key < X(ibegin)) then
idx = -ibegin
return
endif
if (cache_key > X(iend)) then
idx = -iend
return
endif
if (cache_key == X(ibegin)) then
idx = ibegin
return
endif
if (cache_key == X(iend)) then
idx = iend
return
endif
endif
end
subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
integer(key_kind) , intent(in) :: key
real(integral_kind) , intent(out) :: value
integer(cache_key_kind) , intent(in) :: X(sze)
real(integral_kind) , intent(in) :: Y(sze)
integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in
integer(cache_map_size_kind), intent(out) :: idx
integer(cache_map_size_kind) :: istep, ibegin, iend, i
integer(cache_key_kind) :: cache_key
if (sze /= 0) then
continue
else
idx = -1
value = 0.d0
return
endif
cache_key = int(iand(key,map_mask),2)
ibegin = min(ibegin_in,sze)
iend = min(iend_in,sze)
if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
istep = ishft(iend-ibegin,-1)
idx = ibegin + istep
do while (istep > 64)
idx = ibegin + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
idx = ibegin + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
cycle
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
cycle
else
value = Y(idx)
return
endif
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
idx = idx + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
cycle
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
cycle
else
value = Y(idx)
return
endif
else
value = Y(idx)
return
endif
enddo
idx = ibegin
value = Y(idx)
if (min(iend_in,sze) > ibegin+64) then
iend = ibegin+64
do while (cache_key > X(idx))
idx = idx+1
value = Y(idx)
end do
else
do while (cache_key > X(idx))
idx = idx+1
value = Y(idx)
if (idx /= iend) then
cycle
else
exit
endif
end do
endif
if (cache_key /= X(idx)) then
idx = 1-idx
value = 0.d0
endif
return
else
if (cache_key < X(ibegin)) then
idx = -ibegin
value = 0.d0
return
endif
if (cache_key > X(iend)) then
idx = -iend
value = 0.d0
return
endif
if (cache_key == X(ibegin)) then
idx = ibegin
value = Y(idx)
return
endif
if (cache_key == X(iend)) then
idx = iend
value = Y(idx)
return
endif
endif
end
subroutine get_cache_map_n_elements_max(map,n_elements_max)
use map_module
implicit none
! Returns the size of the largest cache_map
type (map_type), intent(in) :: map
integer(cache_map_size_kind), intent(out) :: n_elements_max
integer(map_size_kind) :: i
n_elements_max = 0_cache_map_size_kind
do i=0_8,map%map_size
n_elements_max = max(n_elements_max, map%map(i)%n_elements)
enddo
end
subroutine get_cache_map(map,map_idx,keys,values,n_elements)
use map_module
implicit none
type (map_type), intent(in) :: map
integer(map_size_kind), intent(in) :: map_idx
integer(cache_map_size_kind), intent(inout) :: n_elements
integer(key_kind), intent(out) :: keys(n_elements)
double precision, intent(out) :: values(n_elements)
integer(cache_map_size_kind) :: i
integer(key_kind) :: shift
shift = ishft(map_idx,-map_shift)
n_elements = map%map(map_idx)%n_elements
do i=1,n_elements
keys(i) = map%map(map_idx)%key(i) + shift
values(i) = map%map(map_idx)%value(i)
enddo
end