LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-04-16 13:39:32 +02:00
Code Issues Releases Wiki Activity

AO integrals with ZeroMQ

Browse Source
This commit is contained in:
Anthony Scemama 2015-12-01 16:32:01 +01:00
parent 665638ebc0
commit 94cf7c2b32
8 changed files with 235 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Integrals_Bielec MOGuess
Integrals_Bielec MOGuess

2
src/Integrals_Bielec/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Pseudo Bitmask
Pseudo Bitmask ZMQ

219
src/Integrals_Bielec/ao_bi_integrals.irp.f
View File

@ -301,7 +301,7 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
double precision :: thresh
thresh = ao_integrals_threshold
integer :: n_centers, i
integer :: i
if (ao_overlap_abs(j,l) < thresh) then
buffer_value = 0._integral_kind
@ -329,6 +329,7 @@ end
BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
implicit none
use f77_zmq
use map_module
BEGIN_DOC
! Map of Atomic integrals
@ -345,9 +346,8 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
integer(key_kind),allocatable :: buffer_i(:)
integer,parameter :: size_buffer = 1024*64
real(integral_kind),allocatable :: buffer_value(:)
integer(omp_lock_kind) :: lock
integer :: n_integrals, n_centers, thread_num
integer :: n_integrals, rc
integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax
integral = ao_bielec_integral(1,1,1,1)
@ -363,120 +363,61 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
endif
endif
kk=1
do l = 1, ao_num ! r2
do j = 1, l ! r2
jl_pairs(1,kk) = j
jl_pairs(2,kk) = l
kk += 1
enddo
enddo
PROVIDE progress_bar
call omp_init_lock(lock)
lmax = ao_num*(ao_num+1)/2
print*, 'Providing the AO integrals'
call wall_time(wall_0)
call wall_time(wall_1)
call cpu_time(cpu_1)
call start_progress(lmax,'AO integrals (MB)',0.d0)
!$OMP PARALLEL PRIVATE(i,j,k,l,kk, &
!$OMP integral,buffer_i,buffer_value,n_integrals, &
!$OMP cpu_2,wall_2,i1,j1,thread_num) &
!$OMP DEFAULT(NONE) &
!$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, &
!$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, &
!$OMP ao_overlap_abs,ao_overlap,abort_here, &
!$OMP wall_0,progress_bar,progress_value, &
!$OMP ao_bielec_integral_schwartz)
allocate(buffer_i(size_buffer))
allocate(buffer_value(size_buffer))
n_integrals = 0
!$ thread_num = omp_get_thread_num()
!$OMP DO SCHEDULE(dynamic)
do kk=1,lmax
IRP_IF COARRAY
if (mod(kk-this_image(),num_images()) /= 0) then
cycle
endif
IRP_ENDIF
if (abort_here) then
cycle
endif
if (thread_num == 0) then
progress_bar(1) = kk
endif
j = jl_pairs(1,kk)
l = jl_pairs(2,kk)
j1 = j+ishft(l*l-l,-1)
if (ao_overlap_abs(j,l) < thresh) then
cycle
endif
do k = 1, ao_num ! r1
i1 = ishft(k*k-k,-1)
if (i1 > j1) then
exit
endif
do i = 1, k
i1 += 1
if (i1 > j1) then
exit
endif
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
cycle
endif
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
cycle
endif
!DIR$ FORCEINLINE
integral = ao_bielec_integral(i,k,j,l)
if (abs(integral) < thresh) then
cycle
endif
n_integrals += 1
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
buffer_value(n_integrals) = integral
if (n_integrals > 1024 ) then
if (omp_test_lock(lock)) then
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
n_integrals = 0
call omp_unset_lock(lock)
endif
endif
if (n_integrals == size(buffer_i)) then
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
n_integrals = 0
endif
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(kk)/float(lmax), '% in ', &
wall_2-wall_1, 's', map_mb(ao_integrals_map) ,'MB'
progress_value = dble(map_mb(ao_integrals_map))
endif
endif
enddo
!$OMP END DO NOWAIT
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
deallocate(buffer_i)
deallocate(buffer_value)
!$OMP END PARALLEL
call omp_destroy_lock(lock)
call stop_progress
if (abort_here) then
stop 'Aborting in AO integrals calculation'
integer(ZMQ_PTR) :: zmq_socket_rep_inproc, zmq_socket_push_inproc
zmq_socket_rep_inproc = f77_zmq_socket(zmq_context, ZMQ_REP)
rc = f77_zmq_bind(zmq_socket_rep_inproc, 'inproc://req_rep')
if (rc /= 0) then
stop 'Unable to connect zmq_socket_rep_inproc'
endif
IRP_IF COARRAY
print*, 'Communicating the map'
call communicate_ao_integrals()
IRP_ENDIF COARRAY
integer(ZMQ_PTR) :: thread(0:nproc)
external :: ao_bielec_integrals_in_map_slave, ao_bielec_integrals_in_map_collector
rc = pthread_create( thread(0), ao_bielec_integrals_in_map_collector )
! Create client threads
do i=1,nproc
rc = pthread_create( thread(i), ao_bielec_integrals_in_map_slave )
enddo
character*(64) :: message_string
do l = ao_num, 1, -1
rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0)
print *, l
! TODO : error handling
ASSERT (rc >= 0)
ASSERT (message == 'get_ao_integrals')
rc = f77_zmq_send( zmq_socket_rep_inproc, l, 4, 0)
enddo
do i=1,nproc
rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0)
! TODO : error handling
ASSERT (rc >= 0)
ASSERT (message == 'get_ao_integrals')
rc = f77_zmq_send( zmq_socket_rep_inproc, 0, 4, 0)
enddo
! TODO terminer thread(0)
rc = f77_zmq_unbind(zmq_socket_rep_inproc, 'inproc://req_rep')
do i=1,nproc
rc = pthread_join( thread(i) )
enddo
zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH)
rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull')
if (rc /= 0) then
stop 'Unable to connect zmq_socket_push_inproc'
endif
rc = f77_zmq_send( zmq_socket_push_inproc, -1, 4, ZMQ_SNDMORE)
rc = f77_zmq_send( zmq_socket_push_inproc, 0_key_kind, key_kind, ZMQ_SNDMORE)
rc = f77_zmq_send( zmq_socket_push_inproc, 0_integral_kind, integral_kind, 0)
rc = pthread_join( thread(0) )
print*, 'Sorting the map'
call map_sort(ao_integrals_map)
call cpu_time(cpu_2)
@ -1256,3 +1197,57 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
end
subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
implicit none
use map_module
BEGIN_DOC
! Parallel client for AO integrals
END_DOC
integer, intent(in) :: j,l
integer,intent(out) :: n_integrals
integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num)
real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
integer :: i,k
double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2
double precision :: integral, wall_0
double precision :: thresh
integer :: kk, m, j1, i1
thresh = ao_integrals_threshold
n_integrals = 0
j1 = j+ishft(l*l-l,-1)
do k = 1, ao_num ! r1
i1 = ishft(k*k-k,-1)
if (i1 > j1) then
exit
endif
do i = 1, k
i1 += 1
if (i1 > j1) then
exit
endif
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
cycle
endif
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
cycle
endif
!DIR$ FORCEINLINE
integral = ao_bielec_integral(i,k,j,l)
if (abs(integral) < thresh) then
cycle
endif
n_integrals += 1
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
buffer_value(n_integrals) = integral
enddo
enddo
end

3
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -247,8 +247,7 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
print*, 'MO map initialized'
END_PROVIDER
subroutine insert_into_ao_integrals_map(n_integrals, &
buffer_i, buffer_values)
subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
use map_module
implicit none
BEGIN_DOC

1
src/ZMQ/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@

15
src/ZMQ/README.rst Normal file
View File

@ -0,0 +1,15 @@
===
ZMQ
===
Socket address : defined as an environment variable : QP_RUN_ADDRESS
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

4
src/ZMQ/f77_zmq_module.f90 Normal file
View File

@ -0,0 +1,4 @@
module f77_zmq
include 'f77_zmq.h'
end module

105
src/ZMQ/zmq.irp.f Normal file
View File

@ -0,0 +1,105 @@
use f77_zmq
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
implicit none
BEGIN_DOC
! Context for the ZeroMQ library
END_DOC
zmq_context = f77_zmq_ctx_new ()
END_PROVIDER
BEGIN_PROVIDER [ character*(128), qp_run_address ]
&BEGIN_PROVIDER [ integer, zmq_port_start ]
implicit none
BEGIN_DOC
! Address of the qp_run socket
! Example : tcp://130.120.229.139:12345
END_DOC
character*(128) :: buffer
call getenv('QP_RUN_ADDRESS',buffer)
if (trim(buffer) == '') then
stop 'QP_RUN_ADDRESS environment variable not defined'
endif
print *, trim(buffer)
integer :: i
do i=len(buffer),1,-1
if ( buffer(i:i) == ':') then
qp_run_address = trim(buffer(1:i-1))
read(buffer(i+1:), *) zmq_port_start
exit
endif
enddo
END_PROVIDER
function zmq_port(ishift)
implicit none
integer, intent(in) :: ishift
character*(8) :: zmq_port
write(zmq_port,'(I8)') zmq_port_start+ishift
zmq_port = adjustl(trim(zmq_port))
end
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_to_qp_run_socket ]
implicit none
BEGIN_DOC
! Socket on which the qp_run process replies
END_DOC
integer :: rc
zmq_to_qp_run_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
rc = f77_zmq_connect(zmq_to_qp_run_socket, trim(qp_run_address))
if (rc /= 0) then
stop 'Unable to connect zmq_to_qp_run_socket'
endif
integer :: i
i=4
rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_SNDTIMEO, 120000, i)
if (rc /= 0) then
stop 'Unable to set send timout in zmq_to_qp_run_socket'
endif
rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_RCVTIMEO, 120000, i)
if (rc /= 0) then
stop 'Unable to set recv timout in zmq_to_qp_run_socket'
endif
END_PROVIDER
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_push ]
implicit none
BEGIN_DOC
! Socket on which to push the results (1)
END_DOC
integer :: rc
character*(64) :: address
character*(8), external :: zmq_port
zmq_socket_push = f77_zmq_socket(zmq_context, ZMQ_PUSH)
address = trim(qp_run_address)//':'//zmq_port(1)
rc = f77_zmq_connect(zmq_socket_push, trim(address))
if (rc /= 0) then
stop 'Unable to connect zmq_socket_push'
endif
END_PROVIDER
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_pull ]
implicit none
BEGIN_DOC
! Socket which pulls the results (2)
END_DOC
integer :: rc
character*(64) :: address
character*(8), external :: zmq_port
zmq_socket_pull = f77_zmq_socket(zmq_context, ZMQ_PULL)
address = 'tcp://*:'//zmq_port(2)
rc = f77_zmq_bind(zmq_socket_pull, trim(address))
if (rc /= 0) then
stop 'Unable to connect zmq_socket_pull'
endif
END_PROVIDER