!>
!! \file  bucketSortMod.f90
!!
!! \brief バケットソート（bucket sort）
!!
!! \reference
!! https://ja.wikipedia.org/wiki/バケットソート
!!
!! \author W. Izumida <izumida@cmpt.phys.tohoku.ac.jp>
!!
!! \date 
!! start of the project : 2018/07/21
!!                        2018/07/21 通常ルーチン & 並列化ルーチンの作成
!!                        2018/07/22 メモリ使用に関するテストrunの記述を追加
!!        last modified : 2018/07/23
!!
!<
module bucketSortMod

  implicit none

  private
  public :: bucketSort_parallel, bucketSort

contains

  subroutine bucketSort_parallel(A, order)
    use omp_lib
    implicit none
    integer, intent(in)  :: A(:)
    integer, intent(out) :: order(:)

    integer, parameter            :: iWriteTest = 0

    type bucketType
       integer, allocatable :: idx(:)
    end type bucketType
    type(bucketType), allocatable :: bucket(:)     ! 各バケツ
    integer, allocatable          :: nidx(:,:)     ! 各バケツ内の配列数（各並列毎）
    integer, allocatable          :: pointidx(:,:) ! bucket%idxの引数のポインタ：
                                                   ! pointidx(ib,nPar+1) - 1 が バケツibにおける配列数となる。
    integer, allocatable          :: n(:)          ! 各並列内のAの数
    integer, allocatable          :: pointn(:)     ! 各並列分割に関するポインタ：

    integer                       :: nTot, iC, i, nb, ib
    integer                       :: Amin, Amax
    integer                       :: info
    integer                       :: nPar, iPar
    integer                       :: sum!, sum0

    nPar = omp_get_max_threads()
    if ( iWriteTest == 1 ) write(*,*) "並列数(nPar)=", nPar
    nTot = size(A)
    Amin = minval(A)
    Amax = maxval(A)
    nb = Amax - Amin + 1 ! バケツの個数
    allocate( n(nPar), pointn(nPar+1), nidx(nb,nPar), pointidx(nb,nPar+1), stat=info )
    if ( info /= 0 ) stop 'allocation error for n'
    ! TEST WRITE
    if ( iWriteTest == 1 ) then
       write(*,FMT='("各バケツ内に確保する配列（bucket(ib)%idx）の要素数を")',advance='no')
       write(*,FMT='("記憶するための配列（nidx, pointidx）に使用するメモリ")',advance='no')
       write(*,FMT='("（これは バケツの個数 nb x 並列数 nPar のオーダー）(1) (GB) = ",f20.10)') &
            & real(2 * 4 * nb*nPar, kind(0d0)) / real(1028*1028*1028, kind(0d0))
    end if

    ! 各並列が配列A数のどの箇所を担当するかを与える
    ! (1) 配列の要素数を数える
    n(:) = nTot / nPar
    do iPar = 1, mod(nTot, nPar)
       n(iPar) = n(iPar) + 1
    end do
    ! ポインタ：各並列の担当箇所は pointn(iPar) から pointn(iPar+1)-1 まで
    pointn(1) = 1
    do iPar = 1, nPar
       pointn(iPar+1) = pointn(iPar) + n(iPar)
    end do

    ! 並列処理
    !$OMP parallel do    &
    !$omp default(none), &
    !$OMP private(ib, info),   &
    !$OMP shared(nPar, A, nidx, n, Amin, pointn)
    do iPar = 1, nPar
       ! 各バケツ内の配列（bucket(ib)%idx）の要素数を数える
       nidx(:,iPar) = 0
       do i = 1, n(iPar)
          ib = A( i + pointn(iPar) - 1 ) - Amin + 1 ! 対応するバケツ
          nidx(ib,iPar) = nidx(ib,iPar) + 1
       end do
    end do
    !$omp end parallel do
    ! ポインタ(2)
    pointidx(:,1) = 1
    do iPar = 1, nPar
       pointidx(:,iPar+1) = pointidx(:,iPar) + nidx(:,iPar)
    end do

    ! 各バケツの容器を準備
    allocate( bucket(nb), stat=info )
    if ( info /= 0 ) stop 'allocation error for bucket'
    do ib = 1, nb
       if ( pointidx(ib,nPar+1)-1 >= 1 ) then
          allocate( bucket(ib)%idx(pointidx(ib,nPar+1)-1), stat=info )
          if ( info /= 0 ) stop 'allocation error for bucket'
       end if
    end do
    ! TEST WRITE
    if ( iWriteTest == 1 ) then
       !sum0 = 0
       sum  = 0
       do ib = 1, nb
          !if ( pointidx(ib,nPar+1)-1 == 0 ) sum0 = sum0 + 1
          if ( pointidx(ib,nPar+1)-1 >= 1 ) then
             sum = sum + pointidx(ib,nPar+1)-1
          end if
       end do
       write(*,FMT='("各バケツ内の配列（bucket(ib)%idx）のために使用するメモリ")',advance='no')
       write(*,FMT='("（これは配列 A のためのメモリに等しいはず）(2) (GB) = ",f20.10)') &
            & real(4 * sum, kind(0d0)) / real(1028*1028*1028, kind(0d0))
       !write(*,*) "sum0=", sum0
    end if

    ! 並列処理
    ! 各バケツのidxにAの引数iを格納していく
    !$OMP parallel do    &
    !$omp default(none), &
    !$OMP private(ib),   &
    !$OMP shared(nPar, nidx, n, A, Amin, bucket, pointidx, pointn)
    do iPar = 1, nPar
       nidx(:,iPar) = 0
       do i = 1, n(iPar)
          ib = A( i + pointn(iPar) - 1 ) - Amin + 1 ! 対応するバケツ
          nidx(ib,iPar) = nidx(ib,iPar) + 1
          ! バケツに元々のindexを格納
          bucket(ib)%idx( nidx(ib,iPar) + pointidx(ib,iPar) - 1 ) = i + pointn(iPar) - 1
       end do
    end do
    !$omp end parallel do

    deallocate( nidx )
    deallocate( n )
    deallocate( pointn )

    iC = 0
    do ib = 1, nb
       do i = 1, pointidx(ib,nPar+1) - 1
          iC = iC + 1
          order(iC) = bucket(ib)%idx(i)
       end do
    end do

    deallocate( bucket )
    deallocate( pointidx )

    return
  end subroutine bucketSort_parallel




  subroutine bucketSort(A, order)
    implicit none
    integer, intent(in)  :: A(:)
    integer, intent(out) :: order(:)

    integer, parameter            :: iWriteTest = 1

    type bucketType
       integer, allocatable :: idx(:)
    end type bucketType
    type(bucketType), allocatable :: bucket(:) ! 各バケツ
    integer, allocatable          :: nidx(:)   ! 各バケツ内の配列数
    integer                       :: n, iC, i, nb, ib
    integer                       :: Amin, Amax
    integer                       :: info
    integer                       :: sum

    n = size(A)
    Amin = minval(A)
    Amax = maxval(A)
    nb = Amax - Amin + 1 ! バケツの個数
    allocate( bucket(nb), nidx(nb), stat=info )
    if ( info /= 0 ) stop 'allocation error for bucket'
    ! TEST WRITE
    if ( iWriteTest == 1 ) then
       write(*,FMT='("各バケツ内に確保する配列（bucket(ib)%idx）の要素数を")',advance='no')
       write(*,FMT='("記憶するための配列（nidx）に使用するメモリ(1) (GB) = ",f20.10)') &
            & real(4 * nb, kind(0d0)) / real(1028*1028*1028, kind(0d0))
    end if

    ! 各バケツの配列数を数える
    nidx(:) = 0
    do i = 1, n
       ib = A(i) - Amin + 1     ! 対応するバケツ
       nidx(ib) = nidx(ib) + 1
    end do
    do ib = 1, nb
       if ( nidx(ib) >= 1 ) then
          allocate( bucket(ib)%idx(nidx(ib)), stat=info )
          if ( info /= 0 ) stop 'allocation error for bucket(ib)%idx'
       end if
    end do
    ! TEST WRITE
    if ( iWriteTest == 1 ) then
       sum  = 0
       do ib = 1, nb
          if ( nidx(ib) >= 1 ) then
             sum = sum + nidx(ib)
          end if
       end do
       write(*,FMT='("各バケツ内の配列（bucket(ib)%idx）のために使用するメモリ")',advance='no')
       write(*,FMT='("（これは配列 A のためのメモリに等しいはず）(2) (GB) = ",f20.10)') &
            & real(4 * sum, kind(0d0)) / real(1028*1028*1028, kind(0d0))
    end if

    ! 各バケツに格納していく
    nidx(:) = 0
    do i = 1, n
       ib = A(i) - Amin + 1     ! 対応するバケツ
       nidx(ib) = nidx(ib) + 1
       bucket(ib)%idx(nidx(ib)) = i
    end do

    iC = 0
    do ib = 1, nb
       do i = 1, nidx(ib)
          iC = iC + 1
          order(iC) = bucket(ib)%idx(i)
       end do
    end do

    deallocate( bucket )
    deallocate( nidx )

    return
  end subroutine bucketSort




end module bucketSortMod