From 1ac36ab762f0f42834a915a5cfd72764e6052d18 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 19 Apr 2017 15:31:12 +0200
Subject: [PATCH] Accelerated selection

---
 plugins/Full_CI_ZMQ/run_selection_slave.irp.f |  3 +-
 plugins/Full_CI_ZMQ/selection_buffer.irp.f    | 19 ++----
 src/Utils/sort.irp.f                          | 58 ++++++++++++++-----
 3 files changed, 50 insertions(+), 30 deletions(-)

diff --git a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
index 85b52c30..9ea942a5 100644
--- a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
+++ b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
@@ -45,7 +45,7 @@ subroutine run_selection_slave(thread,iproc,energy)
       if(buf%N == 0) then
         ! Only first time 
         call create_selection_buffer(N, N*2, buf)
-        call create_selection_buffer(N, N*3, buf2)
+        call create_selection_buffer(N, N*2, buf2)
       else
         if(N /= buf%N) stop "N changed... wtf man??"
       end if
@@ -62,7 +62,6 @@ subroutine run_selection_slave(thread,iproc,energy)
         do i=1,buf%cur
           call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i))
         enddo
-        call sort_selection_buffer(buf2)
         buf%mini = buf2%mini
         pt2 = 0d0
         buf%cur = 0
diff --git a/plugins/Full_CI_ZMQ/selection_buffer.irp.f b/plugins/Full_CI_ZMQ/selection_buffer.irp.f
index 84992449..28ceaae3 100644
--- a/plugins/Full_CI_ZMQ/selection_buffer.irp.f
+++ b/plugins/Full_CI_ZMQ/selection_buffer.irp.f
@@ -41,7 +41,6 @@ subroutine sort_selection_buffer(b)
   implicit none
 
   type(selection_buffer), intent(inout) :: b
-  double precision, allocatable::  absval(:)
   integer, allocatable :: iorder(:)
   double precision, pointer :: vals(:)
   integer(bit_kind), pointer :: detmp(:,:,:)
@@ -49,29 +48,23 @@ subroutine sort_selection_buffer(b)
   logical, external :: detEq
   nmwen = min(b%N, b%cur)
 
-
-  allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val)))
-  absval = -dabs(b%val(:b%cur))
+  allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), vals(size(b%val)))
   do i=1,b%cur
     iorder(i) = i
   end do
-  ! Optimal for almost sorted data
-!  call sorted_dnumber(absval, b%cur, i)
-!  if (b%cur/i > 
-!  call insertion_dsort(absval, iorder, b%cur)
-  call dsort(absval, iorder, b%cur)
+  call dsort(b%val, iorder, b%cur)
   do i=1, nmwen
     detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
     detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
     vals(i) = b%val(iorder(i))
   end do
-  do i=nmwen+1, size(vals)
-    vals(i) = 0.d0
-  enddo
+  if (nmwen < b%N) then
+    vals(nmwen+1) = 0.d0
+  endif
   deallocate(b%det, b%val)
   b%det => detmp
   b%val => vals
-  b%mini = max(b%mini,dabs(b%val(b%N)))
+  b%mini = min(b%mini,b%val(1))
   b%cur = nmwen
 end subroutine
 
diff --git a/src/Utils/sort.irp.f b/src/Utils/sort.irp.f
index ba27c0f7..1e271fc0 100644
--- a/src/Utils/sort.irp.f
+++ b/src/Utils/sort.irp.f
@@ -12,23 +12,20 @@ BEGIN_TEMPLATE
   $type                          :: xtmp
   integer                        :: i, i0, j, jmax
   
-  do i=1,isize
+  do i=2,isize
     xtmp = x(i)
     i0 = iorder(i)
-    j = i-1
-    do j=i-1,1,-1
-      if ( x(j) > xtmp ) then
-        x(j+1) = x(j)
-        iorder(j+1) = iorder(j)
-      else
-        exit
-      endif
+    do j = i-1,1,-1
+      if ((x(j) <= xtmp)) exit
+      x(j+1) = x(j)
+      iorder(j+1) = iorder(j)
     enddo
     x(j+1) = xtmp
     iorder(j+1) = i0
   enddo
  end subroutine insertion_$Xsort
 
+
  subroutine heap_$Xsort(x,iorder,isize)
   implicit none
   BEGIN_DOC
@@ -179,7 +176,7 @@ BEGIN_TEMPLATE
   endif
 
   do i=2,isize
-    if (x(i-1) >= x(i)) then
+    if (x(i-1) <= x(i)) then
       n=n+1
     endif
   enddo
@@ -194,6 +191,31 @@ SUBST [ X, type ]
  i2 ; integer*2 ;;
 END_TEMPLATE
 
+!BEGIN_TEMPLATE
+! subroutine $Xsort(x,iorder,isize)
+!  implicit none
+!  BEGIN_DOC
+!  ! Sort array x(isize).
+!  ! iorder in input should be (1,2,3,...,isize), and in output
+!  ! contains the new order of the elements.
+!  END_DOC
+!  integer,intent(in)             :: isize
+!  $type,intent(inout)            :: x(isize)
+!  integer,intent(inout)          :: iorder(isize)
+!  integer                        :: n
+!  call sorted_$Xnumber(x,isize,n)
+!  if ( isize-n < 1000) then
+!    call insertion_$Xsort(x,iorder,isize)
+!  else
+!    call heap_$Xsort(x,iorder,isize)
+!  endif
+! end subroutine $Xsort
+!
+!SUBST [ X, type ]
+!   ; real ;;
+! d ; double precision ;;
+!END_TEMPLATE
+
 BEGIN_TEMPLATE
  subroutine $Xsort(x,iorder,isize)
   implicit none
@@ -207,16 +229,19 @@ BEGIN_TEMPLATE
   integer,intent(inout)          :: iorder(isize)
   integer                        :: n
   call sorted_$Xnumber(x,isize,n)
-  if ( isize-n < 1000) then
+  if (isize == n) then
+    return
+  endif
+  if ( isize < 512+n) then
     call insertion_$Xsort(x,iorder,isize)
   else
-    call heap_$Xsort(x,iorder,isize)
+    call $Yradix_sort(x,iorder,isize,-1)
   endif
  end subroutine $Xsort
 
-SUBST [ X, type ]
-   ; real ;;
- d ; double precision ;;
+SUBST [ X, type, Y ]
+   ; real ; i ;;
+ d ; double precision ; i8 ;;
 END_TEMPLATE
 
 BEGIN_TEMPLATE
@@ -422,6 +447,9 @@ BEGIN_TEMPLATE
     
     i0 = 0_$int_type
     i4 = maxval(x)
+    if (i4 == 0_$type) then
+      return
+    endif
     
     iradix_new = $integer_size-1-leadz(i4)
     mask = ibset(0_$type,iradix_new)