graphutilities.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_ASSEMBLER_PARALLEL_GRAPHUTILITIES_HH
#define DUNE_ASSEMBLER_PARALLEL_GRAPHUTILITIES_HH
 
#include <algorithm>
#include <cstddef>
#if __has_include(<execution>)
#include <execution>
#else
#include <version>
#endif
#include <iostream>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <dune/common/rangeutilities.hh>
 
namespace Dune::Assembler::Experimental {
 
namespace Impl {
 
 
 
/* A class for building a container that looks like std::vector<std::vector<T>>
 * but uses flat storage. The class should be used as follows:
 *
 * // Create with given outer size
 * auto vec = FlatBlockedVector<T>(size);
 *
 * // Reserve size for all entries.
 * // It's allowed to do this in arbitrary order.
 * vec.capacity(2) = ...;
 * vec.capacity(1) = ...;
 * vec.capacity(2) += ...;
 *
 * // Initialize blocks
 * vec.initialize();
 *
 * // Fill blocks in arbitrary order
 * vec[i].push_back(value);
 * vec[i].push_back(value);
 *
 * // Swap or access entries using operator[] or iterators
 * // or resize within reserved capacity
 * vec[i][0] = ...;
 * auto&& vec_i = vec[i];
 * std::sort(vec_i.begin(), vec_i.end());
 * vec[i].resize(1);
 *
 * // Export a range of ranges. This will steal the internal data.
 * auto r = std::move(vec).toRange();
 *
 * // Use r like a nested range
 * for (auto&& ri : r)
 *   for (auto&& rij : ri)
 *     foo(rij);
 */
template <class T, class st = std::size_t>
class FlatBlockedVector
{
public:
  using size_type = st;
 
private:
 
  struct MutableBlock
  {
    void push_back (const T& t)
    {
      data_[end_] = t;
      ++end_;
    }
 
    T& operator[] (size_type i)
    {
      return data_[begin_+i];
    }
 
    auto begin ()
    {
      return data_.begin() + begin_;
    }
 
    auto end ()
    {
      return data_.begin() + end_;
    }
 
    void resize (size_type size)
    {
      end_ = begin_+size;
    }
 
    std::vector<T>& data_;
    size_type& begin_;
    size_type& end_;
  };
 
public:
 
  FlatBlockedVector (size_type size)
    : blockBegin_(size, 0),
    blockEnd_(size, 0)
  {
    blockEnd_.clear();
    blockEnd_.resize(size, 0);
  }
 
  size_type size () const
  {
    return blockBegin_.size();
  }
 
  size_type& capacity (size_type i)
  {
    return blockEnd_[i];
  }
 
  void initialize ()
  {
    auto size = blockBegin_.size();
    blockBegin_[0] = 0;
    for (auto i : Dune::range(std::size_t(1), size))
    {
      blockBegin_[i] = blockBegin_[i-1] + capacity(i-1);
      blockEnd_[i-1] = blockBegin_[i-1];
    }
    data_.resize(blockBegin_[size-1] + capacity(size-1));
    blockEnd_[size-1] = blockBegin_[size-1];
  }
 
  auto operator[] (size_type i) const
  {
    return Dune::transformedRangeView(Dune::range(blockBegin_[i], blockEnd_[i]), [&](auto j) {
      return data_[j];
    });
  }
 
  auto operator[] (size_type i)
  {
    return MutableBlock{data_, blockBegin_[i], blockEnd_[i]};
  }
 
  auto toRange () &&
  {
    auto size = blockBegin_.size();
    return Dune::transformedRangeView(Dune::range(size),
                                      [setVector=std::move(*this)](auto i) {
      return setVector[i];
    });
  }
 
  std::vector<T> data_;
  std::vector<size_type> blockBegin_;
  std::vector<size_type> blockEnd_;
};
 
 
 
// Sort vector and drop duplicate entries
template <class V>
void sortAndUnifyVector (V&& v)
{
  std::sort(v.begin(), v.end());
  auto lastUnique = std::unique(v.begin(), v.end());
  v.resize(lastUnique-v.begin());
}
 
 
 
} // end namespace Impl within Dune::Assembler::Experimental
 
 
 
template <class CliqueList>
auto cliqueListToAdjacency (const CliqueList& cliqueList, std::size_t size, bool verbose=false)
{
  using size_type = std::decay_t<decltype(*(cliqueList.begin()->begin()))>;
 
  // In verbode mode: Compute size of the largest provided clique.
  // This is a lower bound for the clique number (largest existing clique)
  // which itself is a lower bound for the chromatic number.
  // Notice that even for planar 2d grids the element graph is not planar
  // and thus we may need more than 4 colors.
  if (verbose)
  {
    size_type maxCliqueSize = 0;
    for (const auto& clique : cliqueList)
      maxCliqueSize = std::max(size_type(maxCliqueSize), size_type(clique.size()));
    std::cout << "Largest found element clique (lower bound of required colors) has size " << maxCliqueSize << std::endl;
  }
 
  // Theoretically we could use a std::set for each adjacency[k]
  // and simply insert all entries from the cliques into this set.
  // However, simply pushing back entries and removing duplicates
  // afterwards is faster.
  // Furthermore using windows into a flat std::vector<size_type>
  // is faster than using a std::vector<std::vector<size_type>>.
  // The latter is provided by the FlatBlockedVector helper class.
  auto adjacency = Impl::FlatBlockedVector<size_type, size_type>(size);
 
  // We first must compute an upper bound for the
  // degree of each vertex to reserve the appropriate
  // block in the flat vector.
  for (const auto& clique : cliqueList)
    for (auto k : clique)
      adjacency.capacity(k) += clique.size()-1;
  adjacency.initialize();
 
  // Add edges to adjacency list
  for (const auto& clique : cliqueList)
    for (auto k : clique)
      for (auto l : clique)
        if (k!=l)
          adjacency[k].push_back(l);
 
  // Remove duplicate entries in inner container
  {
    auto indices = Dune::range(adjacency.size());
#if __cpp_lib_execution >= 201603L
    std::for_each(std::execution::par, indices.begin(), indices.end(), [&](auto i) {
#else
    std::for_each(indices.begin(), indices.end(), [&](auto i) {
#endif
      Impl::sortAndUnifyVector(adjacency[i]);
    });
  }
 
  return std::move(adjacency).toRange();
}
 
} // end namespace Dune::Assembler::Experimental
 
#endif // DUNE_ASSEMBLER_PARALLEL_GRAPHUTILITIES_HH