Improved parallel scaling of Davidson

plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f

@@ -38,7 +38,8 @@ subroutine run_wf
   integer, external              :: zmq_get_dvector, zmq_get_N_det_generators 
   integer, external              :: zmq_get8_dvector
   integer, external              :: zmq_get_ivector
-  integer, external              :: zmq_get_psi, zmq_get_N_det_selectors
+  integer, external              :: zmq_get_psi, zmq_get_N_det_selectors, zmq_get_psi_bilinear
+  integer, external              :: zmq_get_psi_notouch
   integer, external              :: zmq_get_N_states_diag
 
   zmq_context = f77_zmq_ctx_new ()
@@ -131,25 +132,17 @@ subroutine run_wf
       ! --------
 
       call wall_time(t0)
-      if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
+      if (zmq_get_psi_notouch(zmq_to_qp_run_socket,1) == -1) cycle
       if (zmq_get_N_states_diag(zmq_to_qp_run_socket,1) == -1) cycle
       if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states_diag) == -1) cycle
+      if (zmq_get_psi_bilinear(zmq_to_qp_run_socket,1) == -1) cycle
+      SOFT_TOUCH psi_det psi_coef psi_det_size N_det N_states psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns psi_bilinear_matrix_order
 
       call wall_time(t1)
       if (mpi_master) then
         call write_double(6,(t1-t0),'Broadcast time')
       endif
 
-      call wall_time(t0)
-      if (.True.) then
-        PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
-        PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
-        PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
-        PROVIDE psi_bilinear_matrix_transp_order
-      endif
-      call wall_time(t1)
-      call write_double(6,(t1-t0),'Sort time')
-
       call omp_set_nested(.True.)
       call davidson_slave_tcp(0)
       call omp_set_nested(.False.)

src/Davidson/davidson_parallel.irp.f

@@ -304,6 +304,9 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then
     stop 'Unable to put psi on ZMQ server'
   endif
+  if (zmq_put_psi_bilinear(zmq_to_qp_run_socket,1) == -1) then
+    stop 'Unable to put psi on ZMQ server'
+  endif
   if (zmq_put_N_states_diag(zmq_to_qp_run_socket, 1) == -1) then
     stop 'Unable to put N_states_diag on ZMQ server'
   endif
@@ -381,7 +384,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   integer*8 :: rc8
   double precision :: energy(N_st)
 
-  integer, external :: zmq_put_dvector, zmq_put_psi, zmq_put_N_states_diag
+  integer, external :: zmq_put_dvector, zmq_put_psi, zmq_put_N_states_diag, zmq_put_psi_bilinear
   integer, external :: zmq_put_dmatrix
 
   if (size(u_t) < 8388608) then

src/Determinants/zmq.irp.f

@@ -40,6 +40,174 @@ end
 
 
 
+integer function zmq_get_psi_notouch(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get the wave function from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+
+  integer, external              :: zmq_get_N_states
+  integer, external              :: zmq_get_N_det
+  integer, external              :: zmq_get_psi_det_size
+  integer*8, external            :: zmq_get_psi_det
+  integer*8, external            :: zmq_get_psi_coef
+
+  zmq_get_psi_notouch = 0
+
+  if (zmq_get_N_states(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  if (zmq_get_N_det(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  if (zmq_get_psi_det_size(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  
+  if (size(psi_det,kind=8) /= N_int*2_8*psi_det_size*bit_kind) then
+    deallocate(psi_det)
+    allocate(psi_det(N_int,2,psi_det_size))
+  endif
+
+  if (size(psi_coef,kind=8) /= psi_det_size*N_states) then
+    deallocate(psi_coef)
+    allocate(psi_coef(psi_det_size,N_states))
+  endif
+
+  if (zmq_get_psi_det(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  if (zmq_get_psi_coef(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+
+end
+
+
+integer function zmq_get_psi(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get the wave function from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer, external :: zmq_get_psi_notouch
+  zmq_get_psi = zmq_get_psi_notouch(zmq_to_qp_run_socket, worker_id)
+  SOFT_TOUCH psi_det psi_coef psi_det_size N_det N_states
+
+end
+
+
+
+
+
+integer function zmq_put_psi_bilinear(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put the wave function on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  character*(256)                :: msg
+
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_values
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_rows
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_columns
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_order
+
+  zmq_put_psi_bilinear = 0
+  if (zmq_put_psi_bilinear_matrix_rows(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+  if (zmq_put_psi_bilinear_matrix_columns(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+  if (zmq_put_psi_bilinear_matrix_order(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+  if (zmq_put_psi_bilinear_matrix_values(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+end
+
+
+integer function zmq_get_psi_bilinear(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get the wave function from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_values
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_rows
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_columns
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_order
+
+  zmq_get_psi_bilinear= 0
+
+  if (size(psi_bilinear_matrix_values,kind=8) /= N_det*N_states) then
+    deallocate(psi_bilinear_matrix_values)
+    allocate(psi_bilinear_matrix_values(N_det,N_states))
+  endif
+
+  if (size(psi_bilinear_matrix_rows,kind=8) /= N_det) then
+    deallocate(psi_bilinear_matrix_rows)
+    allocate(psi_bilinear_matrix_rows(N_det))
+  endif
+
+  if (size(psi_bilinear_matrix_columns,kind=8) /= N_det) then
+    deallocate(psi_bilinear_matrix_columns)
+    allocate(psi_bilinear_matrix_columns(N_det))
+  endif
+
+  if (size(psi_bilinear_matrix_order,kind=8) /= N_det) then
+    deallocate(psi_bilinear_matrix_order)
+    allocate(psi_bilinear_matrix_order(N_det))
+  endif
+
+  if (zmq_get_psi_bilinear_matrix_values(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+  if (zmq_get_psi_bilinear_matrix_rows(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+  if (zmq_get_psi_bilinear_matrix_columns(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+  if (zmq_get_psi_bilinear_matrix_order(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+  SOFT_TOUCH psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns psi_bilinear_matrix_order
+
+end
+
+
+
+
+
+
+
 BEGIN_TEMPLATE 
 
 integer function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
@@ -140,11 +308,14 @@ psi_det_size ;;
 
 END_TEMPLATE
 
-integer*8 function zmq_put_psi_det(zmq_to_qp_run_socket,worker_id)
+
+BEGIN_TEMPLATE
+
+integer*8 function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none
   BEGIN_DOC
-! Put psi_det on the qp_run scheduler
+! Put $X on the qp_run scheduler
   END_DOC
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer, intent(in)            :: worker_id
@@ -154,49 +325,22 @@ integer*8 function zmq_put_psi_det(zmq_to_qp_run_socket,worker_id)
   integer*8 :: zmq_put_i8matrix
   integer   :: ni, nj
 
-  if (size(psi_det,kind=8) <= 8388608_8) then
+  if (size($X,kind=8) <= 8388608_8) then
-    ni = size(psi_det,kind=4)
+    ni = size($X,kind=4)
     nj = 1
   else
     ni = 8388608_8
-    nj = int(size(psi_det,kind=8)/8388608_8,4) + 1
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
   endif
-  rc8 = zmq_put_i8matrix(zmq_to_qp_run_socket, 1, 'psi_det', psi_det, ni, nj, size(psi_det,kind=8))
+  rc8 = zmq_put_i8matrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
-  zmq_put_psi_det = rc8
+  zmq_put_$X = rc8
 end
 
-integer*8 function zmq_put_psi_coef(zmq_to_qp_run_socket,worker_id)
+integer*8 function zmq_get_$X(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none
   BEGIN_DOC
-! Put psi_coef on the qp_run scheduler
+! Get $X on the qp_run scheduler
-  END_DOC
-  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
-  integer, intent(in)            :: worker_id
-  integer*8                      :: rc8
-  character*(256)                :: msg
-
-  zmq_put_psi_coef = 0
-
-  integer*8 :: zmq_put_dmatrix
-  integer   :: ni, nj
-
-  if (size(psi_coef,kind=8) <= 8388608_8) then
-    ni = size(psi_coef,kind=4)
-    nj = 1
-  else
-    ni = 8388608
-    nj = int(size(psi_coef,kind=8)/8388608_8,4) + 1
-  endif
-  rc8 = zmq_put_dmatrix(zmq_to_qp_run_socket, 1, 'psi_coef', psi_coef, ni, nj, size(psi_coef,kind=8) )
-  zmq_put_psi_coef = rc8
-end
-
-integer*8 function zmq_get_psi_det(zmq_to_qp_run_socket,worker_id)
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Get psi_det on the qp_run scheduler
   END_DOC
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer, intent(in)            :: worker_id
@@ -206,99 +350,148 @@ integer*8 function zmq_get_psi_det(zmq_to_qp_run_socket,worker_id)
   integer*8 :: zmq_get_i8matrix
   integer   :: ni, nj
 
-  if (size(psi_det,kind=8) <= 8388608_8) then
+  if (size($X,kind=8) <= 8388608_8) then
-    ni = size(psi_det,kind=4)
+    ni = size($X,kind=4)
     nj = 1
   else
     ni = 8388608
-    nj = int(size(psi_det,kind=8)/8388608_8,4) + 1
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
   endif
-  rc8 = zmq_get_i8matrix(zmq_to_qp_run_socket, 1, 'psi_det', psi_det, ni, nj, size(psi_det,kind=8))
+  rc8 = zmq_get_i8matrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
-  zmq_get_psi_det = rc8
+  zmq_get_$X = rc8
 end
 
-integer*8 function zmq_get_psi_coef(zmq_to_qp_run_socket,worker_id)
+SUBST [ X ]
+
+psi_det ;;
+
+END_TEMPLATE
+
+BEGIN_TEMPLATE
+
+integer*8 function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none
   BEGIN_DOC
-! get psi_coef on the qp_run scheduler
+! Put $X on the qp_run scheduler
   END_DOC
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer, intent(in)            :: worker_id
   integer*8                      :: rc8
   character*(256)                :: msg
 
-  zmq_get_psi_coef = 0_8
+  integer*8 :: zmq_put_imatrix
+  integer   :: ni, nj
+
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
+    nj = 1
+  else
+    ni = 8388608_8
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
+  endif
+  rc8 = zmq_put_imatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
+  zmq_put_$X = rc8
+end
+
+integer*8 function zmq_get_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer*8                      :: rc8
+  character*(256)                :: msg
+
+  integer*8 :: zmq_get_imatrix
+  integer   :: ni, nj
+
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
+    nj = 1
+  else
+    ni = 8388608
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
+  endif
+  rc8 = zmq_get_imatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
+  zmq_get_$X = rc8
+end
+
+SUBST [ X ]
+
+psi_bilinear_matrix_rows ;;
+psi_bilinear_matrix_columns ;;
+psi_bilinear_matrix_order ;;
+
+END_TEMPLATE
+
+
+BEGIN_TEMPLATE
+
+integer*8 function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer*8                      :: rc8
+  character*(256)                :: msg
+
+  zmq_put_$X = 0
+
+  integer*8 :: zmq_put_dmatrix
+  integer   :: ni, nj
+
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
+    nj = 1
+  else
+    ni = 8388608
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
+  endif
+  rc8 = zmq_put_dmatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8) )
+  zmq_put_$X = rc8
+end
+
+integer*8 function zmq_get_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! get $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer*8                      :: rc8
+  character*(256)                :: msg
+
+  zmq_get_$X = 0_8
 
   integer*8 :: zmq_get_dmatrix
   integer   :: ni, nj
 
-  if (size(psi_coef,kind=8) <= 8388608_8) then
+  if (size($X,kind=8) <= 8388608_8) then
-    ni = size(psi_coef,kind=4)
+    ni = size($X,kind=4)
     nj = 1
   else
     ni = 8388608
-    nj = int(size(psi_coef,kind=8)/8388608_8,4) + 1
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
   endif
-  rc8 = zmq_get_dmatrix(zmq_to_qp_run_socket, 1, 'psi_coef', psi_coef, ni, nj, size(psi_coef,kind=8) )
+  rc8 = zmq_get_dmatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8) )
-  zmq_get_psi_coef = rc8
+  zmq_get_$X = rc8
 end
 
+SUBST [ X ]
+
+psi_coef ;;
+psi_bilinear_matrix_values ;;
+
+END_TEMPLATE
+
+
 !---------------------------------------------------------------------------
 
 
-integer function zmq_get_psi(zmq_to_qp_run_socket, worker_id)
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Get the wave function from the qp_run scheduler
-  END_DOC
-  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
-  integer, intent(in)            :: worker_id
-
-  integer, external              :: zmq_get_N_states
-  integer, external              :: zmq_get_N_det
-  integer, external              :: zmq_get_psi_det_size
-  integer*8, external            :: zmq_get_psi_det
-  integer*8, external            :: zmq_get_psi_coef
-
-  zmq_get_psi = 0
-
-  if (zmq_get_N_states(zmq_to_qp_run_socket, worker_id) == -1) then
-    zmq_get_psi = -1
-    return
-  endif
-  if (zmq_get_N_det(zmq_to_qp_run_socket, worker_id) == -1) then
-    zmq_get_psi = -1
-    return
-  endif
-  if (zmq_get_psi_det_size(zmq_to_qp_run_socket, worker_id) == -1) then
-    zmq_get_psi = -1
-    return
-  endif
-  
-  if (size(psi_det,kind=8) /= N_int*2_8*psi_det_size*bit_kind) then
-    deallocate(psi_det)
-    allocate(psi_det(N_int,2,psi_det_size))
-  endif
-
-  if (size(psi_coef,kind=8) /= psi_det_size*N_states) then
-    deallocate(psi_coef)
-    allocate(psi_coef(psi_det_size,N_states))
-  endif
-
-  if (zmq_get_psi_det(zmq_to_qp_run_socket, worker_id) == -1_8) then
-    zmq_get_psi = -1
-    return
-  endif
-  if (zmq_get_psi_coef(zmq_to_qp_run_socket, worker_id) == -1_8) then
-    zmq_get_psi = -1
-    return
-  endif
-  SOFT_TOUCH psi_det psi_coef psi_det_size N_det N_states
-
-end
-
-
-
-

src/ZMQ/put_get.irp.f

@@ -608,13 +608,13 @@ integer function zmq_get_int(zmq_to_qp_run_socket, worker_id, name, x)
     include 'mpif.h'
     call MPI_BCAST (zmq_get_int, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
     if (ierr /= MPI_SUCCESS) then
-      print *,  irp_here//': Unable to broadcast zmq_get_i8matrix'
+      print *,  irp_here//': Unable to broadcast zmq_get_int'
       stop -1
     endif
     call MPI_BARRIER(MPI_COMM_WORLD,ierr)
     call MPI_BCAST (x, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
     if (ierr /= MPI_SUCCESS) then
-      print *,  irp_here//': Unable to broadcast zmq_get_i8matrix'
+      print *,  irp_here//': Unable to broadcast zmq_get_int'
       stop -1
     endif
   IRP_ENDIF
@@ -723,7 +723,7 @@ integer function zmq_get_i8matrix(zmq_to_qp_run_socket, worker_id, name, x, size
   integer, intent(in)            :: size_x1, size_x2
   integer*8, intent(in)          :: sze
   character*(*), intent(in)      :: name
-  double precision, intent(out)  :: x(size_x1,size_x2)
+  integer*8, intent(out)         :: x(size_x1,size_x2)
   integer*8                      :: rc, ni
   integer*8                      :: j
   character*(256)                :: msg
@@ -777,7 +777,7 @@ integer function zmq_get_i8matrix(zmq_to_qp_run_socket, worker_id, name, x, size
       stop -1
     endif
     call MPI_BARRIER(MPI_COMM_WORLD,ierr)
-    call broadcast_chunks_double(x, sze)
+    call broadcast_chunks_integer8(x, sze)
   IRP_ENDIF
 
 end
@@ -786,3 +786,124 @@ end
 
 
 
+integer function zmq_put_imatrix(zmq_to_qp_run_socket, worker_id, name, x, size_x1, size_x2, sze)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put a float vector on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  character*(*)                  :: name
+  integer, intent(in)            :: size_x1, size_x2
+  integer*8, intent(in)          :: sze
+  integer, intent(in)            :: x(size_x1, size_x2)
+  integer*8                      :: rc, ni
+  integer*8                      :: j
+  character*(256)                :: msg
+
+  zmq_put_imatrix = 0
+
+  ni = size_x1
+  do j=1,size_x2
+    if (j == size_x2) then
+      ni = sze - (j-1_8)*int(size_x1,8)
+    endif
+    write(msg,'(A,1X,I8,1X,A,I8.8)') 'put_data '//trim(zmq_state), worker_id, trim(name), j
+    rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),ZMQ_SNDMORE)
+    if (rc /= len(trim(msg))) then
+      zmq_put_imatrix = -1
+      print *,  irp_here, 'Failed in put_data', rc, j
+      return 
+    endif
+
+    rc = f77_zmq_send8(zmq_to_qp_run_socket,x(1,j),ni*4_8,0)
+    if (rc /= ni*4_8) then
+      print *,  irp_here, 'Failed in send ', rc, j
+      zmq_put_imatrix = -1
+      return
+    endif
+
+    rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
+    if (msg(1:rc) /= 'put_data_reply ok') then
+      print *,  irp_here, 'Failed in recv ',  rc, j
+      zmq_put_imatrix = -1
+      return
+    endif
+  enddo
+
+end
+
+
+integer function zmq_get_imatrix(zmq_to_qp_run_socket, worker_id, name, x, size_x1, size_x2, sze)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get a float vector from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer, intent(in)            :: size_x1, size_x2
+  integer*8, intent(in)          :: sze
+  character*(*), intent(in)      :: name
+  integer, intent(out)           :: x(size_x1,size_x2)
+  integer*8                      :: rc, ni
+  integer*8                      :: j
+  character*(256)                :: msg
+
+  PROVIDE zmq_state
+  ! Success
+  zmq_get_imatrix = 0
+
+  if (mpi_master) then
+    ni = size_x1
+    do j=1, size_x2
+      if (j == size_x2) then
+        ni = sze - (j-1)*size_x1
+      endif
+      write(msg,'(A,1X,I8,1X,A,I8.8)') 'get_data '//trim(zmq_state), worker_id, trim(name),j
+      rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),0)
+      if (rc /= len(trim(msg))) then
+        zmq_get_imatrix = -1
+        print *,  irp_here, 'rc /= len(trim(msg))', rc, len(trim(msg))
+        go to 10
+      endif
+
+      rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
+      if (msg(1:14) /= 'get_data_reply') then
+        print *,  irp_here, 'msg(1:14) /= get_data_reply', msg(1:14)
+        zmq_get_imatrix = -1
+        go to 10
+      endif
+
+      rc = f77_zmq_recv8(zmq_to_qp_run_socket,x(1,j),ni*4_8,0)
+      if (rc /= ni*4_8) then
+        print *,  irp_here, 'rc /= ni*8', rc, ni*4_8
+        zmq_get_imatrix = -1
+        go to 10
+      endif
+    enddo
+  endif
+
+  10 continue
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    integer :: ierr
+    include 'mpif.h'
+    call MPI_BCAST (zmq_get_imatrix, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      print *,  irp_here//': Unable to broadcast zmq_get_imatrix'
+      stop -1
+    endif
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
+    call broadcast_chunks_integer(x, sze)
+  IRP_ENDIF
+
+end
+
+
+