⚠️ Warning: This is a draft ⚠️

This means it might contain formatting issues, incorrect code, conceptual problems, or other severe issues.

If you want to help to improve and eventually enable this page, please fork RosettaGit's repository and open a merge request on GitHub.

{{collection|Closest pair problem}}

{{works with|Fortran|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 [[Quicksort#Fortran|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.