dbcsr_1d_d_sort_low Subroutine

      recursive subroutine dbcsr_1d_d_sort_low(arr, indices, tmp_arr, tmp_idx)
      !! The actual sort routine.
      !! Only dbcsr_1d_d_sort and itself should call this.

         REAL(kind=real_8), dimension(:), intent(inout) :: arr
         !! the array to sort
         integer, dimension(size(arr)), intent(inout) :: indices
         !! elements-indices before the sort
         REAL(kind=real_8), dimension((size(arr) + 1)/2), intent(inout) :: tmp_arr
         !! scratch space
         integer, dimension((size(arr) + 1)/2), intent(inout) :: tmp_idx
         !! scratch space
         REAL(kind=real_8) :: a
         integer :: t, m, i, j, k
         LOGICAL :: swapped
         ! a,t:  used during swapping of elements in arr and indices

         swapped = .TRUE.

         ! If only a few elements are left we switch to bubble-sort for efficiency.
         if (size(arr) <= 7) then ! 7 seems to be a good choice for the moment
            DO j = size(arr) - 1, 1, -1
               swapped = .FALSE.
               DO i = 1, j
                  IF (arr(i+1) < arr(i)) THEN
                     ! swap arr(i) with arr(i+1)
                     a = arr(i)
                     arr(i) = arr(i + 1)
                     arr(i + 1) = a
                     ! swap indices(i) with indices(i+1)
                     t = indices(i)
                     indices(i) = indices(i + 1)
                     indices(i + 1) = t
                     swapped = .true.
                  END IF
               END DO
               IF (.NOT. swapped) EXIT
            END DO
            return
         end if

         ! split list in half and recursively sort both sublists
         m = (size(arr) + 1)/2 ! index where we going to divide the list in two
         call dbcsr_1d_d_sort_low(arr(1:m), indices(1:m), tmp_arr, tmp_idx)
         call dbcsr_1d_d_sort_low(arr(m + 1:), indices(m + 1:), tmp_arr, tmp_idx)

         ! Check for a special case: Can we just concatenate the two sorted sublists?
         ! This leads to O(n) scaling if the input is already sorted.
         if (arr(m+1) < arr(m)) then
            ! ...no - let's merge the two sorted sublists arr(:m) and arr(m+1:)
            ! Merge will be performed directly in arr. Need backup of first sublist.
            tmp_arr(1:m) = arr(1:m)
            tmp_idx(1:m) = indices(1:m)
            i = 1; ! number of elements consumed from 1st sublist
            j = 1; ! number of elements consumed from 2nd sublist
            k = 1; ! number of elements already merged

            do while (i <= m .and. j <= size(arr) - m)
            if (arr(m+j) < tmp_arr(i)) then
               arr(k) = arr(m + j)
               indices(k) = indices(m + j)
               j = j + 1
            else
               arr(k) = tmp_arr(i)
               indices(k) = tmp_idx(i)
               i = i + 1
            end if
            k = k + 1
            end do

            ! One of the two sublist is now empty.
            ! Copy possibly remaining tail of 1st sublist
            do while (i <= m)
               arr(k) = tmp_arr(i)
               indices(k) = tmp_idx(i)
               i = i + 1
               k = k + 1
            end do

            ! The possibly remaining tail of 2nd sublist is already at the right spot.

         end if

      end subroutine dbcsr_1d_d_sort_low