topological_sort Function

  pure function topological_sort(dag) result(order)
    !! Provide array of vertex numbers ordered in a way that respects dependencies
    type(dag_t), intent(in) :: dag
    integer, allocatable :: order(:)

    call assert(all(dag%vertices(:)%edges_allocated()), "dag_s topological_sort: all(dag%vertices(:)%edges_allocated())")

    block
      type(searched_and_ordered_t) searched_and_ordered
      integer v

      searched_and_ordered = searched_and_ordered_t(s = [integer::], o = [integer::])

      do v = 1, size(dag%vertices)
        if (.not. any(v == searched_and_ordered%s)) &
          searched_and_ordered = depth_first_search(v, [integer::], searched_and_ordered%o)
      end do
      order = searched_and_ordered%o
    end block

  contains

    pure recursive function depth_first_search(v, d, o) result(hybrid)
      integer, intent(in) :: v
      integer, intent(in), dimension(:) :: d, o
      type(searched_and_ordered_t) hybrid

      call assert(.not. any(v == d), "depth_first_search: cycle detected", intrinsic_array_t([v,d]))

      hybrid = searched_and_ordered_t(s = [integer::], o = o)

      associate(dependencies => dag%depends_on(v))
        block
          integer w
          do w = 1, size(dependencies)
            associate(w_dependencies => dependencies(w))
              if (.not. any(w_dependencies == hybrid%s)) hybrid = depth_first_search(w_dependencies, [d, v], hybrid%o)
            end associate
          end do
        end block
      end associate

      if (.not. any(v == hybrid%o)) hybrid%o = [v, hybrid%o]
      hybrid = searched_and_ordered_t(s = [v, hybrid%s], o = hybrid%o)

    end function

  end function topological_sort