bucket_search.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module bucket_search
  use hecmw
  implicit none
 
  private
  public :: bucketdb
  public :: bucketdb_init
  public :: bucketdb_finalize
  public :: bucketdb_setup
  public :: bucketdb_getbucketid
  public :: bucketdb_registerpre
  public :: bucketdb_allocate
  public :: bucketdb_register
  public :: bucketdb_getnumcand
  public :: bucketdb_getcand
 
  integer(kind=kint), parameter :: DEBUG = 0
 
  type bucket
    integer(kind=kint) :: n
    integer(kind=kint) :: n_max
    integer(kind=kint), pointer :: member(:) => null()  
  end type bucket
 
  type bucketdb
    private
    real(kind=kreal) :: x_min(3)                           
    real(kind=kreal) :: x_max(3)                           
    real(kind=kreal) :: d(3)                               
    integer(kind=kint) :: ndiv(3)
    type(bucket), pointer :: buckets(:,:,:) => null()      
    integer(kind=kint) :: n_tot
    integer(kind=kint), pointer :: member_all(:) => null() 
  end type bucketdb
 
contains
 
  subroutine assert(cond, mesg)
    implicit none
    logical, intent(in) :: cond
    character(len=*) :: mesg
    if (debug > 0) then
      if (.not. cond) then
        write(0,*) 'ASSERTION FAILED: ',mesg
        call hecmw_abort( hecmw_comm_get_comm() )
      endif
    endif
  end subroutine assert
 
!!!
!!! routines for type(bucket)
!!!
 
  subroutine bucket_init(bkt)
    implicit none
    type(bucket), intent(inout) :: bkt
    bkt%n = 0
    bkt%n_max = 0
    nullify(bkt%member)
  end subroutine bucket_init
 
  subroutine bucket_finalize(bkt)
    implicit none
    type(bucket), intent(inout) :: bkt
    !if (bkt%n > 0) deallocate(bkt%member)
    nullify(bkt%member)
    bkt%n = 0
    bkt%n_max = 0
  end subroutine bucket_finalize
 
  subroutine bucket_incr_count(bkt)
    implicit none
    type(bucket), intent(inout) :: bkt
    !$omp atomic
    bkt%n = bkt%n + 1
  end subroutine bucket_incr_count
 
  subroutine bucket_assign(bkt, mem)
    implicit none
    type(bucket), intent(inout) :: bkt
    integer(kind=kint), pointer :: mem(:)
    bkt%member => mem
    bkt%n_max = bkt%n
    bkt%n = 0
  end subroutine bucket_assign
 
  subroutine bucket_register(bkt, sid)
    implicit none
    type(bucket), intent(inout) :: bkt
    integer(kind=kint), intent(in) :: sid
    integer(kind=kint) :: idx
    !$omp atomic capture
    bkt%n = bkt%n + 1
    idx = bkt%n
    !$omp end atomic
    call assert(idx <= bkt%n_max, 'bucket_register: too many members')
    bkt%member(idx) = sid
  end subroutine bucket_register
 
  function bucket_get_n(bkt)
    implicit none
    integer(kind=kint) :: bucket_get_n
    type(bucket), intent(in) :: bkt
    bucket_get_n = bkt%n
  end function bucket_get_n
 
  subroutine bucket_get_member(bkt, n, memb)
    implicit none
    type(bucket), intent(in) :: bkt
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(out) :: memb(n)
    call assert(n == bkt%n, 'bucket_get_member: wrong n')
    memb(1:n) = bkt%member(1:n)
  end subroutine bucket_get_member
 
!!!
!!! routines for type(bucketDB)
!!!
 
  subroutine bucketdb_init(bktdb)
    implicit none
    type(bucketdb), intent(inout) :: bktdb
    bktdb%x_min(:) = 0.d0
    bktdb%x_max(:) = 0.d0
    bktdb%d(:) = 0.d0
    bktdb%ndiv(:) = 0
    nullify(bktdb%buckets)
    bktdb%n_tot = 0
    nullify(bktdb%member_all)
  end subroutine bucketdb_init
 
  subroutine bucketdb_finalize(bktdb)
    implicit none
    type(bucketdb), intent(inout) :: bktdb
    integer(kind=kint) :: i, j, k
    if (bktdb%n_tot > 0) then
      deallocate(bktdb%member_all)
      bktdb%n_tot = 0
    endif
    if (any(bktdb%ndiv == 0)) then
      bktdb%ndiv(:) = 0
      return
    endif
    do k = 1, bktdb%ndiv(3)
      do j = 1, bktdb%ndiv(2)
        do i = 1, bktdb%ndiv(1)
          call bucket_finalize(bktdb%buckets(i,j,k))
        enddo
      enddo
    enddo
    deallocate(bktdb%buckets)
    bktdb%ndiv(:) = 0
  end subroutine bucketdb_finalize
 
  subroutine bucketdb_setup(bktdb, x_min, x_max, dmin, n_tot)
    implicit none
    type(bucketdb), intent(inout) :: bktdb
    real(kind=kreal), intent(in) :: x_min(3)  
    real(kind=kreal), intent(in) :: x_max(3)  
    real(kind=kreal), intent(in) :: dmin      
    integer(kind=kint), intent(in) :: n_tot
    real(kind=kreal) :: xrange(3)
    integer(kind=kint) :: i, j, k
    real(kind=kreal), parameter :: eps = 1.d-6
    if (debug >= 1) write(0,*) 'DEBUG: bucketDB_setup', x_min, x_max, dmin, n_tot
    if (associated(bktdb%buckets)) deallocate(bktdb%buckets)
    bktdb%x_min(:) = x_min(:)
    bktdb%x_max(:) = x_max(:)
    xrange(:) = x_max(:) - x_min(:)
    call assert(all(xrange > 0.d0), 'bucketDB_setup: invalid x_min, x_max')
    do i = 1, 3
      bktdb%ndiv(i) = max(floor(xrange(i) / dmin), 1)
      bktdb%d(i) = xrange(i) / bktdb%ndiv(i) * (1.d0 + eps)
    enddo
    if (debug >= 1) write(0,*) 'DEBUG: bucketDB_setup: ndiv, d: ', bktdb%ndiv, bktdb%d
    call assert(all(bktdb%d > 0.d0), 'bucketDB_setup: invalid bktdb%d')
    allocate(bktdb%buckets(bktdb%ndiv(1), bktdb%ndiv(2), bktdb%ndiv(3)), stat=i)
    if( i /= 0 ) then
      write(*,*) 'Allocation error: bktdb%buckets', bktdb%ndiv
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    do k = 1, bktdb%ndiv(3)
      do j = 1, bktdb%ndiv(2)
        do i = 1, bktdb%ndiv(1)
          call bucket_init(bktdb%buckets(i,j,k))
        enddo
      enddo
    enddo
    if (bktdb%n_tot /= n_tot) then
      if (associated(bktdb%member_all)) deallocate(bktdb%member_all)
      allocate(bktdb%member_all(n_tot))
      bktdb%n_tot = n_tot
    endif
  end subroutine bucketdb_setup
 
  function encode_bid(bktdb, baddr)
    implicit none
    integer(kind=kint) :: encode_bid
    type(bucketdb), intent(in) :: bktdb
    integer(kind=kint), intent(in) :: baddr(3)
    if (any(baddr <= 0) .or. any(baddr > bktdb%ndiv)) then
      encode_bid = -1
    else
      encode_bid = &
           (baddr(3)-1) * bktdb%ndiv(1) * bktdb%ndiv(2) + (baddr(2)-1) * bktdb%ndiv(1) + baddr(1)
    endif
  end function encode_bid
 
  function decode_bid(bktdb, bid)
    implicit none
    integer(kind=kint) :: decode_bid(3)
    type(bucketdb), intent(in) :: bktdb
    integer(kind=kint), intent(in) :: bid
    call assert(bid <= bktdb%ndiv(1)*bktdb%ndiv(2)*bktdb%ndiv(3), 'decode_bid: out of range')
    if (bid < 0) then
      decode_bid(:) = -1
    else
      decode_bid(1) = mod(bid-1, bktdb%ndiv(1)) + 1
      decode_bid(2) = mod((bid-1)/bktdb%ndiv(1), bktdb%ndiv(2)) + 1
      decode_bid(3) = (bid-1)/(bktdb%ndiv(1) * bktdb%ndiv(2)) + 1
      call assert(encode_bid(bktdb, decode_bid) == bid, 'decode_bid')
    endif
  end function decode_bid
 
  function bucketdb_getbucketid(bktdb, x)
    implicit none
    integer(kind=kint) :: bucketdb_getbucketid
    type(bucketdb), intent(in) :: bktdb
    real(kind=kreal), intent(in) :: x(3)        
    integer(kind=kint) :: baddr(3)
    integer(kind=kint) :: i
    if (bktdb%n_tot == 0) then
      bucketdb_getbucketid = -1
      return
    endif
    do i = 1, 3
      call assert(bktdb%d(i) > 0.d0, 'bucketDB_getBucketID: bktdb%d(i) is zero')
      baddr(i) = floor((x(i) - bktdb%x_min(i)) / bktdb%d(i)) + 1
    enddo
    if (debug >= 2) write(0,*) '  DEBUG: bucketDB_getBucketID: ',x,baddr
    bucketdb_getbucketid = encode_bid(bktdb, baddr)
  end function bucketdb_getbucketid
 
  subroutine bucketdb_registerpre(bktdb, bid)
    implicit none
    type(bucketdb), intent(inout) :: bktdb
    integer(kind=kint), intent(in) :: bid
    integer(kind=kint) :: baddr(3)
    baddr = decode_bid(bktdb, bid)
    call assert(all(baddr > 0) .and. all(baddr <= bktdb%ndiv), 'bucketDB_register_pre: block ID out of range')
    call bucket_incr_count(bktdb%buckets(baddr(1),baddr(2),baddr(3)))
    if (debug >= 2) write(0,*) '  DEBUG: bucketDB_registerPre: ', baddr
  end subroutine bucketdb_registerpre
 
  subroutine bucketdb_allocate(bktdb)
    implicit none
    type(bucketdb), intent(inout) :: bktdb
    integer(kind=kint) :: i, j, k, count, n
    integer(kind=kint), pointer :: pmemb(:)
    count = 0
    do k = 1, bktdb%ndiv(3)
      do j = 1, bktdb%ndiv(2)
        do i = 1, bktdb%ndiv(1)
          !call bucket_allocate(bktdb%buckets(i,j,k))
          n = bucket_get_n(bktdb%buckets(i,j,k))
          pmemb => bktdb%member_all(count+1:count+n)
          call bucket_assign(bktdb%buckets(i,j,k), pmemb)
          count = count + n
        enddo
      enddo
    enddo
  end subroutine bucketdb_allocate
 
  subroutine bucketdb_register(bktdb, bid, sid)
    implicit none
    type(bucketdb), intent(inout) :: bktdb
    integer(kind=kint), intent(in) :: bid
    integer(kind=kint), intent(in) :: sid
    integer(kind=kint) :: baddr(3)
    baddr = decode_bid(bktdb, bid)
    call assert(all(baddr > 0) .and. all(baddr <= bktdb%ndiv), 'bucketDB_register: block ID our of range')
    call bucket_register(bktdb%buckets(baddr(1),baddr(2),baddr(3)), sid)
    if (debug >= 2) write(0,*) '  DEBUG: bucketDB_register: ', baddr, sid
  end subroutine bucketdb_register
 
  function bucketdb_getnumcand(bktdb, bid)
    implicit none
    integer(kind=kint) :: bucketdb_getnumcand
    type(bucketdb), intent(in) :: bktdb
    integer(kind=kint), intent(in) :: bid
    integer(kind=kint) :: baddr(3), ncand, i, j, k, is, ie, js, je, ks, ke
    if (bid < 0) then
      bucketdb_getnumcand = 0
      return
    endif
    baddr = decode_bid(bktdb, bid)
    ncand = 0
    is = max(baddr(1)-1, 1)
    ie = min(baddr(1)+1, bktdb%ndiv(1))
    js = max(baddr(2)-1, 1)
    je = min(baddr(2)+1, bktdb%ndiv(2))
    ks = max(baddr(3)-1, 1)
    ke = min(baddr(3)+1, bktdb%ndiv(3))
    do k = ks, ke
      do j = js, je
        do i = is, ie
          ncand = ncand + bucket_get_n(bktdb%buckets(i,j,k))
        enddo
      enddo
    enddo
    bucketdb_getnumcand = ncand
    if (debug >= 2) write(0,*) '  DEBUG: bucketDB_getNumCand: ',ncand
  end function bucketdb_getnumcand
 
  subroutine bucketdb_getcand(bktdb, bid, ncand, cand)
    implicit none
    type(bucketdb), intent(in) :: bktdb
    integer(kind=kint), intent(in) :: bid
    integer(kind=kint), intent(in) :: ncand
    integer(kind=kint), intent(out), target :: cand(ncand)
    integer(kind=kint) :: baddr(3), i, j, k, n, cnt, is, ie, js, je, ks, ke
    integer(kind=kint), pointer :: pcand(:)
    if (bid < 0 .or. ncand == 0) return
    baddr = decode_bid(bktdb, bid)
    cnt = 0
    is = max(baddr(1)-1, 1)
    ie = min(baddr(1)+1, bktdb%ndiv(1))
    js = max(baddr(2)-1, 1)
    je = min(baddr(2)+1, bktdb%ndiv(2))
    ks = max(baddr(3)-1, 1)
    ke = min(baddr(3)+1, bktdb%ndiv(3))
    do k = ks, ke
      do j = js, je
        do i = is, ie
          n = bucket_get_n(bktdb%buckets(i,j,k))
          pcand => cand(cnt+1:cnt+n)
          call bucket_get_member(bktdb%buckets(i,j,k), n, pcand)
          cnt = cnt + n
          call assert(cnt <= ncand, 'bucketDB_get_cand: array overflow')
        enddo
      enddo
    enddo
    call assert(cnt == ncand, 'bucketDB_get_cand: count mismatch')
    if (debug >= 3) write(0,*) '    DEBUG: bucketDB_getCand: ',cand
  end subroutine bucketdb_getcand
 
end module bucket_search