graph_coloring.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2013 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_graph_coloring_h
#  define dealii_graph_coloring_h


#  include <deal.II/base/config.h>

#  include <deal.II/base/thread_management.h>

#  include <deal.II/lac/sparsity_tools.h>

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#  include <boost/unordered_map.hpp>
#  include <boost/unordered_set.hpp>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

#  include <functional>
#  include <set>
#  include <vector>


DEAL_II_NAMESPACE_OPEN

namespace GraphColoring
{
  namespace internal
  {
    inline bool
    have_nonempty_intersection(
      const std::vector<types::global_dof_index> &indices1,
      const std::vector<types::global_dof_index> &indices2)
    {
      // we assume that both arrays are sorted, so we can walk
      // them in lockstep and see if we encounter an index that's
      // in both arrays. once we reach the end of either array,
      // we know that there is no intersection
      std::vector<types::global_dof_index>::const_iterator p = indices1.begin(),
                                                           q = indices2.begin();
      while ((p != indices1.end()) && (q != indices2.end()))
        {
          if (*p < *q)
            ++p;
          else if (*p > *q)
            ++q;
          else
            // conflict found!
            return true;
        }

      // no conflict found!
      return false;
    }


    template <typename Iterator>
    std::vector<std::vector<Iterator>>
    create_partitioning(
      const Iterator &                         begin,
      const typename identity<Iterator>::type &end,
      const std::function<std::vector<types::global_dof_index>(
        const Iterator &)> &                   get_conflict_indices)
    {
      // Number of iterators.
      unsigned int n_iterators = 0;

      // Create a map from conflict indices to iterators
      boost::unordered_map<types::global_dof_index, std::vector<Iterator>>
        indices_to_iterators;
      for (Iterator it = begin; it != end; ++it)
        {
          const std::vector<types::global_dof_index> conflict_indices =
            get_conflict_indices(it);
          const unsigned int n_conflict_indices = conflict_indices.size();
          for (unsigned int i = 0; i < n_conflict_indices; ++i)
            indices_to_iterators[conflict_indices[i]].push_back(it);
          ++n_iterators;
        }

      // create the very first zone which contains only the first
      // iterator. then create the other zones. keep track of all the
      // iterators that have already been assigned to a zone
      std::vector<std::vector<Iterator>> zones(1,
                                               std::vector<Iterator>(1, begin));
      std::set<Iterator>                 used_it;
      used_it.insert(begin);
      while (used_it.size() != n_iterators)
        {
          // loop over the elements of the previous zone. for each element of
          // the previous zone, get the conflict indices and from there get
          // those iterators that are conflicting with the current element
          typename std::vector<Iterator>::iterator previous_zone_it(
            zones.back().begin());
          typename std::vector<Iterator>::iterator previous_zone_end(
            zones.back().end());
          std::vector<Iterator> new_zone;
          for (; previous_zone_it != previous_zone_end; ++previous_zone_it)
            {
              const std::vector<types::global_dof_index> conflict_indices =
                get_conflict_indices(*previous_zone_it);

              const unsigned int n_conflict_indices(conflict_indices.size());
              for (unsigned int i = 0; i < n_conflict_indices; ++i)
                {
                  const std::vector<Iterator> &conflicting_elements =
                    indices_to_iterators[conflict_indices[i]];
                  for (unsigned int j = 0; j < conflicting_elements.size(); ++j)
                    {
                      // check that the iterator conflicting with the current
                      // one is not associated to a zone yet and if so, assign
                      // it to the current zone. mark it as used
                      //
                      // we can shortcut this test if the conflicting iterator
                      // is the current iterator
                      if ((conflicting_elements[j] != *previous_zone_it) &&
                          (used_it.count(conflicting_elements[j]) == 0))
                        {
                          new_zone.push_back(conflicting_elements[j]);
                          used_it.insert(conflicting_elements[j]);
                        }
                    }
                }
            }

          // If there are iterators in the new zone, then the zone is added to
          // the partition. Otherwise, the graph is disconnected and we need to
          // find an iterator on the other part of the graph. start the whole
          // process again with the first iterator that hasn't been assigned to
          // a zone yet
          if (new_zone.size() != 0)
            zones.push_back(new_zone);
          else
            for (Iterator it = begin; it != end; ++it)
              if (used_it.count(it) == 0)
                {
                  zones.push_back(std::vector<Iterator>(1, it));
                  used_it.insert(it);
                  break;
                }
        }

      return zones;
    }



    template <typename Iterator>
    void
    make_dsatur_coloring(
      std::vector<Iterator> &             partition,
      const std::function<std::vector<types::global_dof_index>(
        const Iterator &)> &              get_conflict_indices,
      std::vector<std::vector<Iterator>> &partition_coloring)
    {
      partition_coloring.clear();

      // Number of zones composing the partitioning.
      const unsigned int        partition_size(partition.size());
      std::vector<unsigned int> sorted_vertices(partition_size);
      std::vector<int>          degrees(partition_size);
      std::vector<std::vector<types::global_dof_index>> conflict_indices(
        partition_size);
      std::vector<std::vector<unsigned int>> graph(partition_size);

      // Get the conflict indices associated to each iterator. The
      // conflict_indices have to be sorted so we can more easily find conflicts
      // later on
      for (unsigned int i = 0; i < partition_size; ++i)
        {
          conflict_indices[i] = get_conflict_indices(partition[i]);
          std::sort(conflict_indices[i].begin(), conflict_indices[i].end());
        }

      // Compute the degree of each vertex of the graph using the
      // intersection of the conflict indices.
      for (unsigned int i = 0; i < partition_size; ++i)
        for (unsigned int j = i + 1; j < partition_size; ++j)
          // If the two iterators share indices then we increase the degree of
          // the vertices and create an ''edge'' in the graph.
          if (have_nonempty_intersection(conflict_indices[i],
                                         conflict_indices[j]))
            {
              ++degrees[i];
              ++degrees[j];
              graph[i].push_back(j);
              graph[j].push_back(i);
            }

      // Sort the vertices by decreasing degree.
      std::vector<int>::iterator degrees_it;
      for (unsigned int i = 0; i < partition_size; ++i)
        {
          // Find the largest element.
          degrees_it         = std::max_element(degrees.begin(), degrees.end());
          sorted_vertices[i] = degrees_it - degrees.begin();
          // Put the largest element to -1 so it cannot be chosen again.
          *degrees_it = -1;
        }

      // Color the graph.
      std::vector<boost::unordered_set<unsigned int>> colors_used;
      for (unsigned int i = 0; i < partition_size; ++i)
        {
          const unsigned int current_vertex(sorted_vertices[i]);
          bool               new_color(true);
          // Try to use an existing color, i.e., try to find a color which is
          // not associated to one of the vertices linked to current_vertex.
          // Loop over the color.
          for (unsigned int j = 0; j < partition_coloring.size(); ++j)
            {
              // Loop on the vertices linked to current_vertex. If one vertex
              // linked to current_vertex is already using the color j, this
              // color cannot be used anymore.
              bool unused_color(true);
              for (const auto adjacent_vertex : graph[current_vertex])
                if (colors_used[j].count(adjacent_vertex) == 1)
                  {
                    unused_color = false;
                    break;
                  }
              if (unused_color)
                {
                  partition_coloring[j].push_back(partition[current_vertex]);
                  colors_used[j].insert(current_vertex);
                  new_color = false;
                  break;
                }
            }
          // Add a new color.
          if (new_color)
            {
              partition_coloring.push_back(
                std::vector<Iterator>(1, partition[current_vertex]));
              boost::unordered_set<unsigned int> tmp;
              tmp.insert(current_vertex);
              colors_used.push_back(tmp);
            }
        }
    }



    template <typename Iterator>
    std::vector<std::vector<Iterator>>
    gather_colors(
      const std::vector<std::vector<std::vector<Iterator>>> &partition_coloring)
    {
      std::vector<std::vector<Iterator>> coloring;

      // Count the number of iterators in each color.
      const unsigned int partition_size(partition_coloring.size());
      std::vector<std::vector<unsigned int>> colors_counter(partition_size);
      for (unsigned int i = 0; i < partition_size; ++i)
        {
          const unsigned int n_colors(partition_coloring[i].size());
          colors_counter[i].resize(n_colors);
          for (unsigned int j = 0; j < n_colors; ++j)
            colors_counter[i][j] = partition_coloring[i][j].size();
        }

      // Find the partition with the largest number of colors for the even
      // partition.
      unsigned int       i_color(0);
      unsigned int       max_even_n_colors(0);
      const unsigned int colors_size(colors_counter.size());
      for (unsigned int i = 0; i < colors_size; i += 2)
        {
          if (max_even_n_colors < colors_counter[i].size())
            {
              max_even_n_colors = colors_counter[i].size();
              i_color           = i;
            }
        }
      coloring.resize(max_even_n_colors);
      for (unsigned int j = 0; j < colors_counter[i_color].size(); ++j)
        coloring[j] = partition_coloring[i_color][j];

      for (unsigned int i = 0; i < partition_size; i += 2)
        {
          if (i != i_color)
            {
              boost::unordered_set<unsigned int> used_k;
              for (unsigned int j = 0; j < colors_counter[i].size(); ++j)
                {
                  // Find the color in the current partition with the largest
                  // number of iterators.
                  std::vector<unsigned int>::iterator it;
                  it = std::max_element(colors_counter[i].begin(),
                                        colors_counter[i].end());
                  unsigned int min_iterators(static_cast<unsigned int>(-1));
                  unsigned int pos(0);
                  // Find the color of coloring with the least number of colors
                  // among the colors that have not been used yet.
                  for (unsigned int k = 0; k < max_even_n_colors; ++k)
                    if (used_k.count(k) == 0)
                      if (colors_counter[i_color][k] < min_iterators)
                        {
                          min_iterators = colors_counter[i_color][k];
                          pos           = k;
                        }
                  colors_counter[i_color][pos] += *it;
                  // Concatenate the current color with the existing coloring.
                  coloring[pos].insert(
                    coloring[pos].end(),
                    partition_coloring[i][it - colors_counter[i].begin()]
                      .begin(),
                    partition_coloring[i][it - colors_counter[i].begin()]
                      .end());
                  used_k.insert(pos);
                  // Put the number of iterators to the current color to zero.
                  *it = 0;
                }
            }
        }

      // If there is more than one partition, do the same thing that we did for
      // the even partitions to the odd partitions
      if (partition_size > 1)
        {
          unsigned int max_odd_n_colors(0);
          for (unsigned int i = 1; i < partition_size; i += 2)
            {
              if (max_odd_n_colors < colors_counter[i].size())
                {
                  max_odd_n_colors = colors_counter[i].size();
                  i_color          = i;
                }
            }
          coloring.resize(max_even_n_colors + max_odd_n_colors);
          for (unsigned int j = 0; j < colors_counter[i_color].size(); ++j)
            coloring[max_even_n_colors + j] = partition_coloring[i_color][j];

          for (unsigned int i = 1; i < partition_size; i += 2)
            {
              if (i != i_color)
                {
                  boost::unordered_set<unsigned int> used_k;
                  for (unsigned int j = 0; j < colors_counter[i].size(); ++j)
                    {
                      // Find the color in the current partition with the
                      // largest number of iterators.
                      std::vector<unsigned int>::iterator it;
                      it = std::max_element(colors_counter[i].begin(),
                                            colors_counter[i].end());
                      unsigned int min_iterators(static_cast<unsigned int>(-1));
                      unsigned int pos(0);
                      // Find the color of coloring with the least number of
                      // colors among the colors that have not been used yet.
                      for (unsigned int k = 0; k < max_odd_n_colors; ++k)
                        if (used_k.count(k) == 0)
                          if (colors_counter[i_color][k] < min_iterators)
                            {
                              min_iterators = colors_counter[i_color][k];
                              pos           = k;
                            }
                      colors_counter[i_color][pos] += *it;
                      // Concatenate the current color with the existing
                      // coloring.
                      coloring[max_even_n_colors + pos].insert(
                        coloring[max_even_n_colors + pos].end(),
                        partition_coloring[i][it - colors_counter[i].begin()]
                          .begin(),
                        partition_coloring[i][it - colors_counter[i].begin()]
                          .end());
                      used_k.insert(pos);
                      // Put the number of iterators to the current color to
                      // zero.
                      *it = 0;
                    }
                }
            }
        }

      return coloring;
    }
  } // namespace internal


  template <typename Iterator>
  std::vector<std::vector<Iterator>>
  make_graph_coloring(
    const Iterator &                               begin,
    const typename identity<Iterator>::type &      end,
    const std::function<std::vector<types::global_dof_index>(
      const typename identity<Iterator>::type &)> &get_conflict_indices)
  {
    Assert(begin != end,
           ExcMessage(
             "GraphColoring is not prepared to deal with empty ranges!"));

    // Create the partitioning.
    std::vector<std::vector<Iterator>> partitioning =
      internal::create_partitioning(begin, end, get_conflict_indices);

    // Color the iterators within each partition.
    // Run the coloring algorithm on each zone in parallel
    const unsigned int partitioning_size(partitioning.size());
    std::vector<std::vector<std::vector<Iterator>>> partition_coloring(
      partitioning_size);

    Threads::TaskGroup<> tasks;
    for (unsigned int i = 0; i < partitioning_size; ++i)
      tasks += Threads::new_task(&internal::make_dsatur_coloring<Iterator>,
                                 partitioning[i],
                                 get_conflict_indices,
                                 partition_coloring[i]);
    tasks.join_all();

    // Gather the colors together.
    return internal::gather_colors(partition_coloring);
  }

  unsigned int
  color_sparsity_pattern(const SparsityPattern &    sparsity_pattern,
                         std::vector<unsigned int> &color_indices);

} // namespace GraphColoring

DEAL_II_NAMESPACE_CLOSE


//----------------------------   graph_coloring.h ---------------------------
// end of #ifndef dealii_graph_coloring_h
#endif
//----------------------------   graph_coloring.h ---------------------------