Added quicksort

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -258,7 +258,7 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove,
 
     time = omp_get_wtime()
   
-    if(time - timeLast > 1d1 .or. more /= 1) then
+    if(time - timeLast > 3d0 .or. more /= 1) then
       timeLast = time
       do i=1, first_det_of_teeth(1)-1
         if(.not.(actually_computed(i))) then
@@ -331,7 +331,7 @@ end function
 
 BEGIN_PROVIDER [ integer, comb_teeth ]
   implicit none
-  comb_teeth = 50
+  comb_teeth = 100
 END_PROVIDER
 
 
@@ -369,7 +369,7 @@ subroutine get_last_full_tooth(computed, last_tooth)
   
   last_tooth = 0
    combLoop : do i=comb_teeth, 1, -1
-     missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-12) ! /4096
+     missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-14) ! /16384
      do j=first_det_of_teeth(i), first_det_of_teeth(i+1)-1
        if(.not.computed(j)) then
          missing -= 1
@@ -543,7 +543,7 @@ end subroutine
   comb_step = 1d0/dfloat(comb_teeth)
   first_det_of_comb = 1
   do i=1,N_det_generators
-    if(pt2_weight(i)/norm_left < comb_step*.5d0) then
+    if(pt2_weight(i)/norm_left < comb_step*.25d0) then
       first_det_of_comb = i
       exit
     end if

plugins/Full_CI_ZMQ/selection.irp.f

@@ -509,7 +509,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   logical :: ok
   integer :: s1, s2, p1, p2, ib, j, istate
   integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
-  double precision :: e_pert, delta_E, val, Hii, max_e_pert,tmp
+  double precision :: e_pert, delta_E, val, Hii, min_e_pert,tmp
   double precision, external :: diag_H_mat_elem_fock
   
   logical, external :: detEq
@@ -536,7 +536,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
       
       Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
-      max_e_pert = 0d0
+      min_e_pert = 0d0
       
       do istate=1,N_states
         delta_E = E0(istate) - Hii
@@ -545,14 +545,14 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         if (delta_E < 0.d0) then
             tmp = -tmp
         endif
-        e_pert = 0.5d0 * ( tmp - delta_E)
+        e_pert = 0.5d0 * (tmp - delta_E)
         pt2(istate) = pt2(istate) + e_pert
-        max_e_pert = min(e_pert,max_e_pert)
+        min_e_pert = min(e_pert,min_e_pert)
 !        ci(istate) = e_pert / mat(istate, p1, p2)
       end do
       
-      if(dabs(max_e_pert) > buf%mini) then
-        call add_to_selection_buffer(buf, det, max_e_pert)
+      if(min_e_pert <= buf%mini) then
+        call add_to_selection_buffer(buf, det, min_e_pert)
       end if
     end do
   end do

plugins/Full_CI_ZMQ/selection_buffer.irp.f

@@ -25,7 +25,7 @@ subroutine add_to_selection_buffer(b, det, val)
   double precision, intent(in) :: val
   integer :: i
 
-  if(dabs(val) >= b%mini) then
+  if(val <= b%mini) then
     b%cur += 1
     b%det(1:N_int,1:2,b%cur) = det(1:N_int,1:2)
     b%val(b%cur) = val
@@ -41,33 +41,25 @@ 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(:,:,:)
   integer :: i, nmwen
   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 = b%val(:b%cur)
+  allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)))
   do i=1,b%cur
     iorder(i) = i
   end do
-  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
-  deallocate(b%det, b%val)
+  deallocate(b%det)
   b%det => detmp
-  b%val => vals
-  b%mini = max(b%mini,dabs(b%val(b%N)))
+  b%mini = min(b%mini,b%val(b%N))
   b%cur = nmwen
+  deallocate(iorder)
 end subroutine
 

src/Utils/sort.irp.f

@@ -27,6 +27,56 @@ BEGIN_TEMPLATE
   enddo
  end subroutine insertion_$Xsort
 
+ subroutine quick_$Xsort(x, iorder, isize)
+  implicit none
+  BEGIN_DOC
+  ! Sort array x(isize) using the quicksort algorithm.
+  ! 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)
+  call rec_$X_quicksort(x,iorder,isize,1,isize)
+ end
+
+ recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last)
+  implicit none
+  integer, intent(in)            :: isize, first, last
+  integer,intent(inout)          :: iorder(isize)
+  $type, intent(inout)           :: x(isize)
+  $type                          :: c, tmp
+  integer                        :: itmp
+  integer                        :: i, j
+   
+  c = x( ishft(first+last,-1) )
+  i = first
+  j = last
+  do
+    do while (x(i) < c)
+      i=i+1
+    end do
+    do while (c < x(j))
+      j=j-1
+    end do
+    if (i >= j) exit
+    tmp  = x(i)
+    x(i) = x(j)
+    x(j) = tmp
+    itmp      = iorder(i)
+    iorder(i) = iorder(j)
+    iorder(j) = itmp
+    i=i+1
+    j=j-1
+  enddo
+  if (first < i-1) then
+    call rec_$X_quicksort(x, iorder, isize, first, i-1)
+  endif
+  if (j+1 < last) then
+    call rec_$X_quicksort(x, iorder, isize, j+1, last)
+  endif
+ end
+
  subroutine heap_$Xsort(x,iorder,isize)
   implicit none
   BEGIN_DOC
@@ -206,10 +256,11 @@ BEGIN_TEMPLATE
 !  if (isize == n) then
 !    return
 !  endif
-  if ( isize < 16) then
+  if ( isize < 32) then
     call insertion_$Xsort(x,iorder,isize)
   else
-    call heap_$Xsort(x,iorder,isize)
+!    call heap_$Xsort(x,iorder,isize)
+    call quick_$Xsort(x,iorder,isize)
   endif
  end subroutine $Xsort