From cb1227a9a98060a136728962b5e86e47119a314b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Sep 2017 20:23:50 +0200 Subject: [PATCH] OK --- plugins/FourIdx/four_idx.irp.f | 44 ++++++ plugins/FourIdx/four_index.irp.f | 147 +++++++++++++++++++++ src/Determinants/two_body_dm_map.irp.f | 2 +- src/Integrals_Bielec/mo_bi_integrals.irp.f | 10 +- src/Utils/map_module.f90 | 49 ++++++- 5 files changed, 243 insertions(+), 9 deletions(-) create mode 100644 plugins/FourIdx/four_idx.irp.f create mode 100644 plugins/FourIdx/four_index.irp.f diff --git a/plugins/FourIdx/four_idx.irp.f b/plugins/FourIdx/four_idx.irp.f new file mode 100644 index 00000000..de5927bf --- /dev/null +++ b/plugins/FourIdx/four_idx.irp.f @@ -0,0 +1,44 @@ +program FourIdx + use map_module + implicit none + BEGIN_DOC +! Performs a four index transformation of the two-electron integrals + END_DOC + + type(map_type) :: test_map + integer(key_kind) :: key_max + integer(map_size_kind) :: sze + + call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max) + sze = key_max + call map_init(test_map,sze) + + call four_index_transform(ao_integrals_map,test_map, & + mo_coef, size(mo_coef,1), & + 1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num, & + 1, 1, 1, 1, mo_tot_num, mo_tot_num, mo_tot_num, mo_tot_num) + + integer :: i,j,k,l + real(integral_kind) :: integral1, integral2 + + provide mo_bielec_integrals_in_map + + do i=1,mo_tot_num + do j=1,mo_tot_num + do k=1,mo_tot_num + do l=1,mo_tot_num + call bielec_integrals_index(i,j,k,l,key_max) + call map_get(test_map,key_max,integral1) + call map_get(mo_integrals_map,key_max,integral2) + if (dabs(integral2) >=1.d-10 ) then + if (dabs(integral1 / integral2 -1.d0) > .001d0) then + print *, i,j,k,l + print *, integral1, integral2 + print *, '' + endif + endif + enddo + enddo + enddo + enddo +end diff --git a/plugins/FourIdx/four_index.irp.f b/plugins/FourIdx/four_index.irp.f new file mode 100644 index 00000000..eba99f2c --- /dev/null +++ b/plugins/FourIdx/four_index.irp.f @@ -0,0 +1,147 @@ +subroutine four_index_transform(map_a,map_c,matrix_B,LDB, & + i_start, j_start, k_start, l_start, & + i_end , j_end , k_end , l_end , & + a_start, b_start, c_start, d_start, & + a_end , b_end , c_end , d_end ) + implicit none + use map_module + BEGIN_DOC +! Performs a four-index transformation of map_a(N^4) into map_c(M^4) using b(NxM) +! C_{abcd} = \sum_{ijkl} A_{ijkl}.B_{ia}.B_{jb}.B_{kc}.B_{ld} +! Loops run over *_start->*_end + END_DOC + type(map_type), intent(in) :: map_a + type(map_type), intent(inout) :: map_c + integer, intent(in) :: LDB + double precision, intent(in) :: matrix_B(LDB,*) + integer, intent(in) :: i_start, j_start, k_start, l_start + integer, intent(in) :: i_end , j_end , k_end , l_end + integer, intent(in) :: a_start, b_start, c_start, d_start + integer, intent(in) :: a_end , b_end , c_end , d_end + + double precision, allocatable :: T(:,:,:), U(:,:,:), V(:,:,:) + integer :: i_max, j_max, k_max, l_max + integer :: i_min, j_min, k_min, l_min + integer :: i, j, k, l + integer :: a, b, c, d + double precision, external :: get_ao_bielec_integral + integer(key_kind) :: idx + real(integral_kind) :: tmp + integer(key_kind), allocatable :: key(:) + real(integral_kind), allocatable :: value(:) + + + + i_min = min(i_start,a_start) + i_max = max(i_end ,a_end ) + j_min = min(j_start,b_start) + j_max = max(j_end ,b_end ) + k_min = min(k_start,c_start) + k_max = max(k_end ,c_end ) + l_min = min(l_start,d_start) + l_max = max(l_end ,d_end ) + + ASSERT (0 < i_max) + ASSERT (0 < j_max) + ASSERT (0 < k_max) + ASSERT (0 < l_max) + ASSERT (LDB >= i_max) + ASSERT (LDB >= j_max) + ASSERT (LDB >= k_max) + ASSERT (LDB >= l_max) + + allocate( T(i_min:i_max,j_min:j_max,k_min:k_max), & + U(i_min:i_max,j_min:j_max,k_min:k_max), & + V(i_min:i_max,j_min:j_max,k_min:k_max), & + key(i_max*j_max*k_max), & + value(i_max*j_max*k_max) ) + + do d=d_start,d_end + U = 0.d0 + print *, d + do l=l_start,l_end + if (dabs(matrix_B(l,d)) < 1.d-10) then + cycle + endif + do k=k_start,k_end + do j=j_start,j_end + do i=i_start,i_end + call bielec_integrals_index(i,j,k,l,idx) + call map_get(map_a,idx,tmp) + T(i,j,k) = tmp + enddo + enddo + enddo + + V = 0.d0 + do a=a_start,a_end + do k=k_start,k_end + do j=j_start,j_end + do i=i_start,i_end + V(j,k,a) = V(j,k,a) + T(i,j,k)*matrix_B(i,a) + enddo + enddo + enddo + enddo +! call DGEMM('T','N', (j_end-j_start+1),(k_end-k_start+1), & +! (i_end-i_start+1), 1.d0, & +! T, size(T,1)* + + T = 0.d0 + do b=b_start,b_end + do a=a_start,a_end + do k=k_start,k_end + do j=j_start,j_end + T(k,a,b) = T(k,a,b) + V(j,k,a)*matrix_B(j,b) + enddo + enddo + enddo + enddo + + V = 0.d0 + do c=c_start,c_end + do b=b_start,b_end + do a=a_start,a_end + do k=k_start,k_end + V(a,b,c) = V(a,b,c) + T(k,a,b)*matrix_B(k,c) + enddo + enddo + enddo + enddo + + do c=c_start,c_end + do b=b_start,b_end + do a=a_start,a_end +! do c=c_start,c_end +! do b=b_start,d +! do a=a_start,min(b,c) + U(a,b,c) = U(a,b,c) + V(a,b,c) * matrix_B(l,d) + enddo + enddo + enddo + + enddo + + idx = 0_8 + do c=c_start,c_end + do b=b_start,b_end + do a=a_start,a_end +! do c=c_start,c_end +! do b=b_start,d +! do a=a_start,min(b,c) + if (dabs(U(a,b,c)) < 1.d-15) then + cycle + endif + idx = idx+1_8 + call bielec_integrals_index(a,b,c,d,key(idx)) + value(idx) = U(a,b,c) + enddo + enddo + enddo + call map_append(map_c, key, value, idx) + call map_sort(map_c) + call map_unique(map_c) + + enddo + +end diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f index aa8f630b..2228b1b5 100644 --- a/src/Determinants/two_body_dm_map.irp.f +++ b/src/Determinants/two_body_dm_map.irp.f @@ -187,7 +187,7 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl) print*,'n_elements = ',n_elements call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& real(mo_integrals_threshold,integral_kind)) - call map_unique(two_body_dm_ab_map) + call map_merge(two_body_dm_ab_map) deallocate(buffer_i,buffer_value) diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 05eb8dff..84cfd228 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -146,7 +146,7 @@ subroutine set_integrals_jj_into_map enddo call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& real(mo_integrals_threshold,integral_kind)) - call map_unique(mo_integrals_map) + call map_merge(mo_integrals_map) end subroutine set_integrals_exchange_jj_into_map @@ -167,7 +167,7 @@ subroutine set_integrals_exchange_jj_into_map enddo call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& real(mo_integrals_threshold,integral_kind)) - call map_unique(mo_integrals_map) + call map_merge(mo_integrals_map) end @@ -458,7 +458,7 @@ subroutine add_integrals_to_map(mask_ijkl) real(mo_integrals_threshold,integral_kind)) deallocate(buffer_i, buffer_value) !$OMP END PARALLEL - call map_unique(mo_integrals_map) + call map_merge(mo_integrals_map) call wall_time(wall_2) call cpu_time(cpu_2) @@ -773,7 +773,7 @@ subroutine add_integrals_to_map_three_indices(mask_ijk) real(mo_integrals_threshold,integral_kind)) deallocate(buffer_i, buffer_value) !$OMP END PARALLEL - call map_unique(mo_integrals_map) + call map_merge(mo_integrals_map) call wall_time(wall_2) call cpu_time(cpu_2) @@ -1035,7 +1035,7 @@ subroutine add_integrals_to_map_no_exit_34(mask_ijkl) ! print*, 'Communicating the map' ! call communicate_mo_integrals() !IRP_ENDIF - call map_unique(mo_integrals_map) + call map_merge(mo_integrals_map) call wall_time(wall_2) call cpu_time(cpu_2) diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index ac16f97e..29f7440c 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -13,7 +13,7 @@ module map_module ! cache_map using a binary search ! ! When using the map_update subroutine to build the map, -! the map_unique subroutine +! the map_merge subroutine ! should be called before getting data from the map. use omp_lib @@ -274,7 +274,7 @@ subroutine map_sort(map) end -subroutine cache_map_unique(map) +subroutine cache_map_merge(map) use map_module implicit none type (cache_map_type), intent(inout) :: map @@ -298,6 +298,28 @@ subroutine cache_map_unique(map) 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) + endif + enddo + map%n_elements = j + +end + subroutine cache_map_shrink(map,thr) use map_module implicit none @@ -338,6 +360,27 @@ subroutine map_unique(map) end +subroutine map_merge(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_merge(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 @@ -402,7 +445,7 @@ subroutine map_update(map, key, value, sze, thr) 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_merge(local_map) call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements) call cache_map_shrink(local_map,thr) endif