Compiler Explorer

Source code

!> Float point kind
module nnps_kinds

implicit none

private
    public :: rk

integer, parameter :: rk = kind(0.0)

end module nnps_kinds

!> int_vector integer vector
module nnps_int_vector

implicit none

private
    public :: int_vector

!> int_vector integer vector
    type int_vector
        integer :: len = 0  !! true vector length
        integer, allocatable :: items(:)  !! integer vector
    contains
        procedure :: push_back
        procedure, private :: extend
    end type int_vector

contains

!> push back 1 item
    pure subroutine push_back(self, item)
        class(int_vector), intent(inout) :: self
        integer, intent(in) :: item

if (allocated(self%items)) then
            if (self%len == size(self%items)) call self%extend()
            self%len = self%len + 1
            self%items(self%len) = item
        else
            allocate (self%items(1), source=item)
            self%len = 1
        end if

end subroutine push_back

!> extend
    pure subroutine extend(self)
        class(int_vector), intent(inout) :: self
        integer, allocatable :: tmp(:)

allocate (tmp(size(self%items)))
        self%items = [self%items, tmp]  ! address of self%items may not be changed, which is good

end subroutine extend

end module nnps_int_vector

!> key-value
module nnps_key_value

use nnps_int_vector, only: int_vector
    implicit none

private
    public :: key_values

!> key-value pair
    type key_value
        integer :: key(3)  !! key
        type(int_vector) :: value  !! value
    end type key_value

!> key-value pairs
    type key_values
        type(key_value), allocatable :: items(:)  !! key-value pair
    contains
        procedure :: push_back, get_value
    end type key_values

contains

!> constructor
    pure type(key_value) function key_value_constructor(key, value) result(self)
        integer, intent(in) :: key(3)
        integer, intent(in) :: value

self%key = key
        call self%value%push_back(value)

end function key_value_constructor

!> push back key-value pair
    pure subroutine push_back(self, key, value, stat)
        class(key_values), intent(inout) :: self
        integer, intent(in) :: key(3)
        integer, intent(in) :: value
        logical, intent(out) :: stat
        integer :: i

if (allocated(self%items)) then
            do i = 1, size(self%items)
                if (all(self%items(i)%key == key)) then
                    call self%items(i)%value%push_back(value)
                    if (self%items(i)%value%len == 1) then
                        stat = .true.
                    else
                        stat = .false.
                    end if
                    return
                end if
            end do
            self%items = [self%items, key_value_constructor(key, value)]
            stat = .true.
        else
            allocate (self%items(1), source=key_value_constructor(key, value))
            stat = .true.
        end if

end subroutine push_back

!> get value @todo improve performance
    subroutine get_value(self, key, ptr)
        class(key_values), intent(in), target :: self
        integer, intent(in) :: key(3)
        integer, pointer, intent(inout) :: ptr(:)
        integer :: i

if (.not. allocated(self%items)) return
        do i = 1, size(self%items)
            if (all(self%items(i)%key == key)) then
                ptr => self%items(i)%value%items(1:self%items(i)%value%len)
                return
            end if
        end do

end subroutine get_value

end module nnps_key_value

!> spatial Hashing
module nnps_spatial_hashing

use nnps_key_value, only: key_values
    implicit none

private
    public :: shash_tbl

!> spatial Hashing
    type shash_tbl
        type(key_values), allocatable :: buckets(:)  !! Buckets
    contains
        procedure :: allocate => shash_tbl_allocate, set, hash
    end type shash_tbl

contains

!> Hashing function
    pure integer function hash(self, key)
        class(shash_tbl), intent(in) :: self
        integer, intent(in) :: key(3)

hash = modulo(ieor(ieor(73856093*key(1), 19349663*key(2)), 83492791*key(3)), &
                      size(self%buckets))

end function hash

!> Allocate
    pure subroutine shash_tbl_allocate(self, m)
        class(shash_tbl), intent(inout) :: self
        integer, intent(in) :: m

allocate (self%buckets(0:m - 1))

end subroutine shash_tbl_allocate

!> Push
    pure subroutine set(self, key, value, stat)
        class(shash_tbl), intent(inout) :: self
        integer, intent(in) :: key(3), value
        logical, intent(out) :: stat

call self%buckets(self%hash(key))%push_back(key, value, stat)

end subroutine set

end module nnps_spatial_hashing

program main

use nnps_kinds, only: rk
    use nnps_spatial_hashing, only: shash_tbl
    use nnps_int_vector, only: int_vector
    implicit none
    real(rk), dimension(:, :), allocatable :: loc
    integer :: m, i, ik(3)
    real(rk) :: min(3)
    logical :: lstat
    type(shash_tbl) :: tbl
    type(int_vector) :: iks
    integer, pointer :: value(:)

m = 2000
    allocate (loc(3, m))
    call tbl%allocate(2*m)
    call randu(loc, -100.0_rk, 100.0_rk)

open (1, file='debug.txt')
    write (*, *) "indexing..."
    write (1, *) "indexing..."

min = minval(loc, 2) - sqrt(epsilon(1.0_rk))
    do i = 1, m
        ik = ceiling((loc(:, i) - min)/0.6_rk)
        call tbl%set(key=ik, value=i, stat=lstat)
        if (lstat) then
            call iks%push_back(ik(1))
            call iks%push_back(ik(2))
            call iks%push_back(ik(3))
        end if
        write (1, *) "(key)i,j,k=", ik, "(value)index=", i
    end do

write (*, *) "searching..."
    write (1, *) "searching..."
    do i = 1, iks%len, 3
        nullify (value)
        call tbl%buckets(tbl%hash(iks%items(i:i + 2)))%get_value(iks%items(i:i + 2), value)
        if (associated(value)) then
            if (any(value > m) .or. any(value < 1)) then
                write (1, *) "(key)i,j,k=", iks%items(i:i + 2), "(value)index=", value(:)
                write (*, *) "(key)i,j,k=", iks%items(i:i + 2), "(value)index=", value(:)
                error stop "index out of range, bug occurred, using `/libs:dll` option would &
                &solve this problem, output file is `debug.txt`"
            end if
        end if
    end do

contains

!> uniform distribution
    impure elemental subroutine randu(x, min, max)
        real(rk), intent(in) :: min, max
        real(rk), intent(out) :: x

call random_number(x)
        x = min + x*(max - min)

end subroutine randu

end program main