diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index f835f11f..3b737f5d 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -371,13 +371,16 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) ! i, j for k,l fixed. END_DOC integer, intent(in) :: k,l, sze - logical, intent(out) :: out_array(sze,sze) + double precision, intent(out) :: out_array(sze,sze) type(map_type), intent(inout) :: map integer :: i,j,kk,ll,m integer(key_kind),allocatable :: hash(:) integer ,allocatable :: pairs(:,:), iorder(:) + real(integral_kind), allocatable :: tmp_val(:) + PROVIDE mo_bielec_integrals_in_map - allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze)) + allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), & + tmp_val(sze*sze)) kk=0 out_array = 0.d0 @@ -401,16 +404,16 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) call i2radix_sort(hash,iorder,kk,-1) endif - call map_exists_many(mo_integrals_map, hash, kk) + call map_get_many(mo_integrals_map, hash, tmp_val, kk) do ll=1,kk m = iorder(ll) i=pairs(1,m) j=pairs(2,m) - out_array(i,j) = (hash(ll) /= 0_8) + out_array(i,j) = tmp_val(ll) enddo - deallocate(pairs,hash,iorder) + deallocate(pairs,hash,iorder,tmp_val) end integer*8 function get_mo_map_size() diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index 1c89355d..24f5a328 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -496,13 +496,16 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx) 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) - if (idx > 0) then - value = map%value(idx) - else - value = 0._integral_kind - endif +! 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 @@ -612,7 +615,7 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) istep = ishft(iend-ibegin,-1) idx = ibegin + istep - do while (istep > 32) + do while (istep > 16) idx = ibegin + istep if (cache_key < X(idx)) then iend = idx @@ -649,8 +652,8 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) endif enddo idx = ibegin - if (min(iend_in,sze) > ibegin+64) then - iend = ibegin+64 + if (min(iend_in,sze) > ibegin+16) then + iend = ibegin+16 !DIR$ VECTOR ALIGNED do while (cache_key > X(idx)) idx = idx+1 @@ -693,6 +696,126 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) 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 = iand(key,map_mask) + 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 > 16) + 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+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) + 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