module mod_weighted_graph
  implicit none
  integer(8), parameter :: long_infty = lshift(1_8, 62)
  integer, parameter :: int_infty = lshift(1, 30)

  type t_edge
    integer :: from = 0
    integer :: to = 0
    integer :: capacity = 0
    integer :: reverse = 0
    integer(8) :: weight = 0
  end type

  type t_list
    type(t_edge), allocatable :: at(:)
  end type

  type weighted_graph
    type(t_list), allocatable :: list(:)
  contains
    procedure :: dijkstra => dijkstra
    procedure :: ford_fulkerson => ford_fulkerson
    procedure :: prim => prim
  end type

  type t_node
    type(t_node), private, pointer :: left => null(), right => null()
    integer :: rank = 1
    integer :: vertex, from
    integer(8) :: cost
  end type

  type leftist_heap
    integer :: size = 0
    type(t_node), private, pointer :: root => null()
  end type

  interface weighted_graph
    module procedure :: newwg0
  end interface weighted_graph

contains
  
  type(t_edge) function newe0(from, to, weight, capacity, reverse) result(res)
    integer, intent(in) :: from, to
    integer(8), intent(in) :: weight
    integer, optional, intent(in) :: capacity, reverse
    res%from = from
    res%to = to
    res%weight = weight
    if (present(capacity)) res%capacity = capacity
    if (present(reverse)) res%reverse = reverse
  end

  type(weighted_graph) function newwg0(n, m, a, b, weight, capacity, undirected) result(res)
    integer, intent(in) :: n, m, a(:), b(:)
    integer(8), intent(in) :: weight(:)
    integer, optional, intent(in) :: capacity(:)
    logical, optional, intent(in) :: undirected
    logical :: u, cap
    integer :: c(n), i
    u = merge(undirected, .true., present(undirected))
    cap = present(capacity)
    c = 0
    do i = 1, m
      c(a(i)) = c(a(i)) + 1
      if (cap) c(b(i)) = c(b(i)) + 1
      if (u) then
        c(b(i)) = c(b(i)) + 1
        if (cap) c(a(i)) = c(a(i)) + 1
      end if
    end do
    allocate(res%list(n))
    do i = 1, n
      allocate(res%list(i)%at(c(i)))
    end do
    do i = m, 1, -1
      if (cap) then
        res%list(a(i))%at(c(a(i))) = newe0(a(i), b(i), weight(i), capacity(i), c(b(i)))
        res%list(b(i))%at(c(b(i))) = newe0(b(i), a(i), weight(i), 0, c(a(i)))
        c(a(i)) = c(a(i)) - 1
        c(b(i)) = c(b(i)) - 1
      else
        res%list(a(i))%at(c(a(i))) = newe0(a(i), b(i), weight(i))
        c(a(i)) = c(a(i)) - 1
      end if
      if (u) then
        if (cap) then
          res%list(b(i))%at(c(b(i))) = newe0(b(i), a(i), weight(i), capacity(i), c(a(i)))
          res%list(a(i))%at(c(a(i))) = newe0(a(i), b(i), weight(i), 0, c(b(i)))
          c(a(i)) = c(a(i)) - 1
          c(b(i)) = c(b(i)) - 1
        else
          res%list(b(i))%at(c(b(i))) = newe0(b(i), a(i), weight(i))
          c(b(i)) = c(b(i)) - 1
        end if
      end if
    end do
  end

  logical function less(a, b) result(res)
    type(t_node), pointer, intent(in) :: a, b
    res = a%cost < b%cost
  end

  subroutine swap(a, b)
    type(t_node), pointer, intent(inout) :: a, b
    type(t_node), pointer :: c
    c => a
    a => b
    b => c
  end

  function newn0(vertex, cost) result(res)
    integer, intent(in) :: vertex
    integer(8), intent(in) :: cost
    type(t_node), pointer :: res
    res => null()
    allocate(res)
    res%vertex = vertex
    res%cost = cost
  end

  recursive function meld_node(a, b) result(res)
    type(t_node), pointer, intent(in) :: a, b
    type(t_node), pointer :: res
    logical :: f
    if (.not.associated(a)) then
      res => b; return
    end if
    if (.not.associated(b)) then
      res => a; return
    end if
    if (less(a, b)) then
      res => a
      res%right => meld_node(res%right, b)
    else
      res => b
      res%right => meld_node(res%right, a)
    end if
    f = .not.associated(res%left)
    if (.not.f) f = res%left%rank < res%right%rank
    if (f) call swap(res%left, res%right)
    res%rank = merge(res%right%rank, 0, associated(res%right)) + 1
  end

  function poll(this) result(res)
    class(leftist_heap), intent(inout) :: this
    type(t_node), pointer :: res
    res => this%root
    this%root => meld_node(this%root%left, this%root%right)
    this%size = this%size - 1
  end

  subroutine offer(this, vertex, cost, from)
    class(leftist_heap), intent(inout) :: this
    integer, intent(in) :: vertex
    integer(8), intent(in) :: cost
    integer, intent(in), optional :: from
    type(t_node), pointer :: node
    node => newn0(vertex, cost)
    if (present(from)) node%from = from
    this%root => meld_node(this%root, node)
    this%size = this%size + 1
  end

  function dijkstra(this, start) result(res)
    class(weighted_graph), intent(in) :: this
    integer, intent(in) :: start
    integer :: i, v
    integer(8) :: res(size(this%list)), cost
    type(t_node), pointer :: node
    type(leftist_heap) :: pq
    type(t_edge) :: e
    res = long_infty
    res(start) = 0
    call offer(pq, start, res(start))
    do while (pq%size > 0)
      node => poll(pq)
      v = node%vertex
      cost = node%cost
      deallocate(node)
      if (cost > res(v)) cycle
      do i = 1, size(this%list(v)%at)
        e = this%list(v)%at(i)
        if (res(e%to) > res(v) + e%weight) then
          res(e%to) = res(v) + e%weight
          call offer(pq, e%to, res(e%to))
        end if
      end do
    end do
  end

  integer function ford_fulkerson(this, start, goal) result(res)
    class(weighted_graph), intent(inout) :: this
    integer, intent(in) :: start, goal
    integer :: dflow
    logical :: used(size(this%list))
    res = 0
    do
      used = .false.
      dflow = ffdfs(this, used, start, goal, int_infty)
      if (dflow == 0) return
      res = res + dflow
    end do
  end

  recursive integer function ffdfs(this, used, u, goal, dflow) result(res)
    class(weighted_graph), intent(inout) :: this
    logical, intent(inout) :: used(:)
    integer, intent(in) :: u, goal
    integer, intent(in) :: dflow
    integer :: i, v, capacity, r
    if (u == goal) then
      res = dflow
      return
    end if
    used(u) = .true.
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)%to
      capacity = this%list(u)%at(i)%capacity
      if (used(v) .or. capacity <= 0) cycle
      res = ffdfs(this, used, v, goal, min(dflow, capacity))
      r = this%list(u)%at(i)%reverse
      if (res > 0) then
        this%list(u)%at(i)%capacity = this%list(u)%at(i)%capacity - res
        this%list(v)%at(r)%capacity = this%list(v)%at(r)%capacity + res
        return
      end if
    end do
    res = 0
  end

  integer(8) function bellman_ford(n, m, a, b, weight, start, goal) result(res)
    integer, intent(in) :: n, m, a(:), b(:), start, goal
    integer(8), intent(in) :: weight(:)
    type(weighted_graph) :: gs, gg
    logical :: useds(n), usedg(n), changed
    integer :: i, u, v, cnt, x(m), y(m), step
    integer(8) :: w(m), cost(n)
    gs = newwg0(n, m, a, b, weight, undirected = .false.)
    gg = newwg0(n, m, b, a, weight, undirected = .false.)
    useds = .false.
    usedg = .false.
    call bfdfs(gs, useds, start)
    call bfdfs(gg, usedg, goal)
    cnt = 0
    do i = 1, m
      u = a(i)
      v = b(i)
      if (useds(u) .and. useds(v) .and. usedg(u) .and. usedg(v)) then
        cnt = cnt + 1
        x(cnt) = u
        y(cnt) = v
        w(cnt) = weight(i)
      end if
    end do
    cost = long_infty
    cost(start) = 0
    res = long_infty
    step = 0
    changed = .true.
    do while (changed)
      changed = .false.
      do i = 1, cnt
        if (cost(y(i)) > cost(x(i)) + w(i)) then
          cost(y(i)) = cost(x(i)) + w(i)
          changed = .true.
        end if
      end do
      step = step + 1
      if (step > n) then
        res = -long_infty; exit
      end if
    end do
    if (res == long_infty) res = cost(goal)
  end

  recursive subroutine bfdfs(this, used, u)
    type(weighted_graph), intent(in) :: this
    logical, intent(inout) :: used(:)
    integer, intent(in) :: u
    integer :: i, v
    used(u) = .true.
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)%to
      if (used(v)) cycle
      call bfdfs(this, used, v)
    end do
  end

  type(weighted_graph) function prim(this) result(res)
    class(weighted_graph), intent(in) :: this
    logical :: used(size(this%list))
    integer :: i, v, from, a(size(this%list) - 1), b(size(this%list) - 1), cnt
    integer(8) :: cost, weight(size(this%list) - 1)
    type(t_node), pointer :: node
    type(leftist_heap) :: pq
    type(t_edge) :: e
    do i = 1, size(this%list(1)%at)
      e = this%list(1)%at(i)
      call offer(pq, e%to, e%weight, 1)
    end do
    used = .false.
    used(1) = .true.
    cnt = 0
    do while (pq%size > 0)
      node => poll(pq)
      v = node%vertex
      from = node%from
      cost = node%cost
      deallocate(node)
      if (used(v)) cycle
      used(v) = .true.
      cnt = cnt + 1
      a(cnt) = from
      b(cnt) = v
      weight(cnt) = cost
      do i = 1, size(this%list(v)%at)
        e = this%list(v)%at(i)
        if (used(e%to)) cycle
        call offer(pq, e%to, e%weight, v)
      end do
    end do
    res = newwg0(size(this%list), size(this%list) - 1, a, b, weight)
  end

end module mod_weighted_graph
program two_d_plane_two_n_points
  use mod_weighted_graph
  integer :: n, m = 0, i, j
  integer, dimension(100) :: a = 0, b = 0, c = 0, d = 0
  integer, dimension(20000) :: x = 0, y = 0
  type(weighted_graph) :: g
  read(*,*) n
  do i = 1, n
    read(*,*) a(i), b(i)
  end do
  do i = 1, n
    read(*,*) c(i), d(i)
  end do
  do i = 1, n
    m = m + 1
    x(m) = 1
    y(m) = i + 1
    m = m + 1
    x(m) = i + n + 1
    y(m) = 2 * n + 2
  end do
  do i = 1, n
    do j = 1, n
      if (a(i) < c(j) .and. b(i) < d(j)) then
        m = m + 1
        x(m) = i + 1
        y(m) = j + n + 1
      end if
    end do
  end do
  g = weighted_graph(2 * n + 2, m, x, y, &
    [(0_8, i = 1, m)], [(1, i = 1, m)], .false.)
  write(*,'(i0)') g%ford_fulkerson(1, 2 * n + 2)
end program two_d_plane_two_n_points

Measured memory_bytes: 3116