# Closest-pair problem/Fortran

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Closest-pair problem/Fortran is part of Closest pair problem. You may find other members of Closest pair problem at Category:Closest pair problem.
Works with: Fortran version 95 and later

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.