Closest-pair problem/Fortran
(Redirected from Closest pair problem/Fortran)
Closest-pair problem/Fortran is part of Closest pair problem. You may find other members of Closest pair problem at Category:Closest pair problem.
This module defines the point type and implements only two operations on "vectors" ("difference" and "length")
module Points_Module
implicit none
type point
real :: x, y
end type point
interface operator (-)
module procedure pt_sub
end interface
interface len
module procedure pt_len
end interface
public :: point
private :: pt_sub, pt_len
contains
function pt_sub(a, b) result(c)
type(point), intent(in) :: a, b
type(point) :: c
c = point(a%x - b%x, a%y - b%y)
end function pt_sub
function pt_len(a) result(l)
type(point), intent(in) :: a
real :: l
l = sqrt((a%x)**2 + (a%y)**2)
end function pt_len
end module Points_Module
Then we need a modified version of the fortran quicksort able to handle "points" and that can use a custom comparator.
module qsort_module
use Points_Module
implicit none
contains
recursive subroutine qsort(a, comparator)
type(point), intent(inout) :: a(:)
interface
integer function comparator(aa, bb)
use Points_Module
type(point), intent(in) :: aa, bb
end function comparator
end interface
integer :: split
if(size(a) > 1) then
call partition(a, split, comparator)
call qsort(a(:split-1), comparator)
call qsort(a(split:), comparator)
end if
end subroutine qsort
subroutine partition(a, marker, comparator)
type(point), intent(inout) :: a(:)
integer, intent(out) :: marker
interface
integer function comparator(aa, bb)
use Points_Module
type(point), intent(in) :: aa, bb
end function comparator
end interface
type(point) :: pivot, temp
integer :: left, right
left = 0
right = size(a) + 1
pivot = a(size(a)/2)
do while (left < right)
right = right - 1
do while (comparator(a(right), pivot) > 0)
right = right - 1
end do
left = left + 1
do while (comparator(a(left), pivot) < 0)
left = left + 1
end do
if ( left < right ) then
temp = a(left)
a(left) = a(right)
a(right) = temp
end if
end do
if (left == right) then
marker = left + 1
else
marker = left
end if
end subroutine partition
end module qsort_module
The module containing the custom comparators.
module Order_By_XY
use Points_Module
implicit none
contains
integer function order_by_x(aa, bb)
type(point), intent(in) :: aa, bb
if ( aa%x < bb%x ) then
order_by_x = -1
elseif (aa%x > bb%x) then
order_by_x = 1
else
order_by_x = 0
end if
end function order_by_x
integer function order_by_y(aa, bb)
type(point), intent(in) :: aa, bb
if ( aa%y < bb%y ) then
order_by_y = -1
elseif (aa%y > bb%y) then
order_by_y = 1
else
order_by_y = 0
end if
end function order_by_y
end module Order_By_XY
The closest pair functions' module.
module ClosestPair
use Points_Module
use Order_By_XY
use qsort_module
implicit none
private :: closest_pair_real
contains
function closest_pair_simple(p, pair) result(distance)
type(point), dimension(:), intent(in) :: p
type(point), dimension(:), intent(out), optional :: pair
real :: distance
type(point), dimension(2) :: cp
integer :: i, j
real :: d
if ( size(P) < 2 ) then
distance = huge(0.0)
else
distance = len(p(1) - p(2))
cp = (/ p(1), p(2) /)
do i = 1, size(p) - 1
do j = i+1, size(p)
d = len(p(i) - p(j))
if ( d < distance ) then
distance = d
cp = (/ p(i), p(j) /)
end if
end do
end do
if ( present(pair) ) pair = cp
end if
end function closest_pair_simple
function closest_pair(p, pair) result(distance)
type(point), dimension(:), intent(in) :: p
type(point), dimension(:), intent(out), optional :: pair
real :: distance
type(point), dimension(2) :: dp
type(point), dimension(size(p)) :: xp, yp
integer :: i
xp = p
yp = p
call qsort(xp, order_by_x)
call qsort(yp, order_by_y)
distance = closest_pair_real(xp, yp, dp)
if ( present(pair) ) pair = dp
end function closest_pair
recursive function closest_pair_real(xp, yp, pair) result(distance)
type(point), dimension(:), intent(in) :: xp, yp
type(point), dimension(:), intent(out) :: pair
real :: distance
type(point), dimension(2) :: pairl, pairr
type(point), dimension(:), allocatable :: ys
type(point), dimension(:), allocatable :: pl, yl
type(point), dimension(:), allocatable :: pr, yr
real :: dl, dr, dmin, xm, d
integer :: ns, i, k, j, midx
if ( size(xp) <= 3 ) then
distance = closest_pair_simple(xp, pair)
else
midx = ceiling(size(xp)/2.0)
allocate(ys(size(xp)))
allocate(pl(midx), yl(midx))
allocate(pr(size(xp)-midx), yr(size(xp)-midx))
pl = xp(1:midx)
pr = xp((midx+1):size(xp))
xm = xp(midx)%x
k = 1; j = 1
do i = 1, size(yp)
if ( yp(i)%x > xm ) then
! guard ring; this is an "idiosyncrasy" that should be fixed in a
! smarter way
if ( k <= size(yr) ) yr(k) = yp(i)
k = k + 1
else
! guard ring (see above)
if ( j <= size(yl) ) yl(j) = yp(i)
j = j + 1
end if
end do
dl = closest_pair_real(pl, yl, pairl)
dr = closest_pair_real(pr, yr, pairr)
pair = pairr
dmin = dr
if ( dl < dr ) then
pair = pairl
dmin = dl
end if
ns = 0
do i = 1, size(yp)
if ( abs(yp(i)%x - xm) < dmin ) then
ns = ns + 1
ys(ns) = yp(i)
end if
end do
distance = dmin
do i = 1, ns-1
k = i + 1
do while ( k <= ns .and. abs(ys(k)%y - ys(i)%y) < dmin )
d = len(ys(k) - ys(i))
if ( d < distance ) then
distance = d
pair = (/ ys(k), ys(i) /)
end if
k = k + 1
end do
end do
deallocate(ys)
deallocate(pl, yl)
deallocate(pr, yr)
end if
end function closest_pair_real
end module ClosestPair
Testing:
program TestClosestPair
use ClosestPair
implicit none
integer, parameter :: n = 10000
integer :: i
real :: x, y
type(point), dimension(n) :: points
type(point), dimension(2) :: p
real :: dp, dr
! init the random generator here if needed
do i = 1, n
call random_number(x)
call random_number(y)
points(i) = point(x*20.0-10.0, y*20.0-10.0)
end do
dp = closest_pair_simple(points, p)
print *, "sim ", dp
dr = closest_pair(points, p)
print *, "rec ", dr
end program TestClosestPair
Time gave 2.92user 0.00system 0:02.94elapsed for brute force, and 0.02user 0.00system 0:00.03elapsed for the other one.