graph.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/adjacency/base.hpp>
#include <kernel/adjacency/adjactor.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/exception.hpp>
#include <kernel/util/omp_util.hpp>
 
// includes, system
#include <vector>
 
namespace FEAT
{
  namespace Adjacency
  {
 
    // forward declaration
    class Permutation;
 
    class Graph
    {
    public:
      using IndexVector = std::vector<Index>;
 
      typedef IndexVector::const_iterator ImageIterator;
 
      static constexpr std::uint64_t magic = 0x5052474A44413346ull; // "F3ADJGRP"
 
    protected:
      Index _num_nodes_image;
 
      IndexVector _domain_ptr;
 
      IndexVector _image_idx;
 
    public:
      Graph();
 
      explicit Graph(
        Index num_nodes_domain,
        Index num_nodes_image,
        Index num_indices_image);
 
      explicit Graph(
        Index num_nodes_domain,
        Index num_nodes_image,
        Index num_indices_image,
        const Index* domain_ptr,
        const Index* image_idx);
 
      explicit Graph(
        Index num_nodes_image,
        const IndexVector& domain_ptr,
        const IndexVector& image_idx);
 
      explicit Graph(
        Index num_nodes_image,
        IndexVector&& domain_ptr,
        IndexVector&& image_idx);
 
      template<typename Adjactor_>
      explicit Graph(RenderType render_type, const Adjactor_& adjactor) :
        _num_nodes_image(0),
        _domain_ptr(),
        _image_idx()
      {
        switch(render_type)
        {
        case RenderType::as_is:
        case RenderType::as_is_sorted:
          _render_as_is(adjactor);
          if(render_type == RenderType::as_is_sorted)
            this->sort_indices();
          break;
 
        case RenderType::injectify:
        case RenderType::injectify_sorted:
          _render_injectify(adjactor);
          if(render_type == RenderType::injectify_sorted)
            this->sort_indices();
          break;
 
        case RenderType::transpose:
        case RenderType::transpose_sorted:
          _render_transpose(adjactor);
          // transpose is automatically sorted
          break;
 
        case RenderType::injectify_transpose:
        case RenderType::injectify_transpose_sorted:
          _render_injectify_transpose(adjactor);
          // transpose is automatically sorted
          break;
 
        default:
          XABORTM("Invalid render_type parameter!");
        }
      }
 
      template<typename Adjactor1_, typename Adjactor2_>
      explicit Graph(RenderType render_type, const Adjactor1_& adjactor1, const Adjactor2_& adjactor2) :
        _num_nodes_image(0),
        _domain_ptr(),
        _image_idx()
      {
        switch(render_type)
        {
        case RenderType::as_is:
        case RenderType::as_is_sorted:
          _render_as_is(adjactor1, adjactor2);
          if(render_type == RenderType::as_is_sorted)
            this->sort_indices();
          break;
 
        case RenderType::injectify:
        case RenderType::injectify_sorted:
          _render_injectify(adjactor1, adjactor2);
          if(render_type == RenderType::injectify_sorted)
            this->sort_indices();
          break;
 
        case RenderType::transpose:
        case RenderType::transpose_sorted:
          _render_transpose(adjactor1, adjactor2);
          // transpose is automatically sorted
          break;
 
        case RenderType::injectify_transpose:
        case RenderType::injectify_transpose_sorted:
          _render_injectify_transpose(adjactor1, adjactor2);
          // transpose is automatically sorted
          break;
 
        default:
          XABORTM("Invalid render_type parameter!");
        }
      }
 
      Graph(Graph&& other);
 
      Graph& operator=(Graph&& other);
 
      explicit Graph(const Graph& other, const Permutation& domain_perm, const Permutation& image_perm);
 
      explicit Graph(const std::vector<char>& buffer);
 
      virtual ~Graph();
 
      void clear();
 
      Graph clone() const
      {
        if(!_domain_ptr.empty())
          return Graph(_num_nodes_image, _domain_ptr, _image_idx);
        else
          return Graph();
      }
 
      Index degree(Index domain_node) const
      {
        ASSERTM(domain_node < get_num_nodes_domain(), "Domain node index out of range");
        return _domain_ptr[domain_node+1] - _domain_ptr[domain_node];
      }
 
      Index degree() const;
 
      Index* get_domain_ptr()
      {
        return _domain_ptr.data();
      }
 
      const Index* get_domain_ptr() const
      {
        return _domain_ptr.data();
      }
 
      Index* get_image_idx()
      {
        return _image_idx.data();
      }
 
      const Index* get_image_idx() const
      {
        return _image_idx.data();
      }
 
      Index get_num_indices() const
      {
        return Index(_image_idx.size());
      }
 
      void sort_indices();
 
      void permute_indices(const Adjacency::Permutation& inv_perm);
 
      std::vector<char> serialize() const;
 
      /* *************************************************** */
      /*  R E N D E R   F U N C T I O N   T E M P L A T E S  */
      /* *************************************************** */
    private:
 
      template<typename Adjactor_>
      void _render_as_is(const Adjactor_& adj)
      {
        typedef typename Adjactor_::ImageIterator AImIt;
 
        // get counts
        //_num_nodes_domain = adj.get_num_nodes_domain();
        _num_nodes_image = adj.get_num_nodes_image();
        Index num_indices_image = 0;
 
        // allocate pointer vector
        _domain_ptr.resize(adj.get_num_nodes_domain() + 1);
        _domain_ptr[0] = 0;
 
        // count number of adjacencies and build pointer vector
        FEAT_PRAGMA_OMP(parallel for schedule(dynamic, 1000))
        for(Index i = 0; i < adj.get_num_nodes_domain(); ++i)
        {
          Index num_indices_here = Index(0);
          AImIt cur(adj.image_begin(i));
          AImIt end(adj.image_end(i));
          for(; cur != end; ++cur)
          {
            ++num_indices_here;
          }
          _domain_ptr[i+1] = num_indices_here;
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj.get_num_nodes_domain()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj.get_num_nodes_domain()];
 
        // allocate and build index vector
        _image_idx.resize(num_indices_image);
        FEAT_PRAGMA_OMP(parallel for schedule(dynamic, 1000))
        for(Index i = 0; i < adj.get_num_nodes_domain(); ++i)
        {
          Index k = _domain_ptr[i];
          AImIt cur(adj.image_begin(i));
          AImIt end(adj.image_end(i));
          for(; cur != end; ++cur, ++k)
          {
            _image_idx[k] = *cur;
          }
        }
      }
 
      template<typename Adjactor_>
      void _render_injectify(const Adjactor_& adj)
      {
        // get counts
        _num_nodes_image = adj.get_num_nodes_image();
        Index num_indices_image = 0;
 
        // allocate pointer vector
        _domain_ptr.resize(adj.get_num_nodes_domain() + 1);
        _domain_ptr[0] = 0;
 
        FEAT_PRAGMA_OMP(parallel)
        {
          // allocate auxiliary mask vector
          std::vector<char> vidx_mask(adj.get_num_nodes_image(), 0);
          char* idx_mask = vidx_mask.data();
 
          // count number of adjacencies and build pointer vector
          FEAT_PRAGMA_OMP(for schedule(dynamic, 1000))
          for(Index i = 0; i < adj.get_num_nodes_domain(); ++i)
          {
            Index num_indices_here = Index(0);
            for(auto it = adj.image_begin(i); it != adj.image_end(i); ++it)
            {
              if(idx_mask[*it] == 0)
              {
                ++num_indices_here;
                idx_mask[*it] = 1;
              }
            }
            _domain_ptr[i+1] = num_indices_here;
            for(auto it = adj.image_begin(i); it != adj.image_end(i); ++it)
              idx_mask[*it] = 0;
          }
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj.get_num_nodes_domain()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj.get_num_nodes_domain()];
 
        // allocate and build index vector
        _image_idx.resize(num_indices_image);
 
        FEAT_PRAGMA_OMP(parallel)
        {
          std::vector<char> vidx_mask(adj.get_num_nodes_image(), 0);
          char* idx_mask = vidx_mask.data();
          FEAT_PRAGMA_OMP(for schedule(dynamic, 1000))
          for(Index i = 0; i < adj.get_num_nodes_domain(); ++i)
          {
            Index k = _domain_ptr[i];
            for(auto it = adj.image_begin(i); it != adj.image_end(i); ++it)
            {
              if(idx_mask[*it] == 0)
              {
                _image_idx[k] = *it;
                ++k;
                idx_mask[*it] = 1;
              }
            }
            for(auto it = adj.image_begin(i); it != adj.image_end(i); ++it)
              idx_mask[*it] = 0;
          }
        }
      }
 
      template<typename Adjactor_>
      void _render_transpose(const Adjactor_& adj)
      {
        typedef typename Adjactor_::ImageIterator AImIt;
 
        // get counts
        _num_nodes_image = adj.get_num_nodes_domain();
        Index num_indices_image = 0;
 
        // allocate and format pointer vector
        _domain_ptr.resize(adj.get_num_nodes_image() + 1, Index(0));
 
        // count number of adjacencies
        for(Index j(0); j < adj.get_num_nodes_domain(); ++j)
        {
          AImIt cur(adj.image_begin(j));
          AImIt end(adj.image_end(j));
          for(; cur != end; ++cur)
          {
            ++_domain_ptr[(*cur) + 1];
          }
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj.get_num_nodes_image()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj.get_num_nodes_image()];
 
        // allocate and build index vector
        _image_idx.resize(num_indices_image);
        std::vector<Index*> vimg_ptr(adj.get_num_nodes_image(), nullptr);
        Index** image_ptr = vimg_ptr.data();
        Index* image_idx = _image_idx.data();
        for(Index i(0); i < adj.get_num_nodes_image(); ++i)
        {
          image_ptr[i] = &image_idx[_domain_ptr[i]];
        }
 
        for(Index j(0); j < adj.get_num_nodes_domain(); ++j)
        {
          AImIt cur(adj.image_begin(j));
          AImIt end(adj.image_end(j));
          for(; cur != end; ++cur)
          {
            Index*& idx = image_ptr[*cur];
            *idx = j;
            ++idx;
          }
        }
      }
 
      template<typename Adjactor_>
      void _render_injectify_transpose(const Adjactor_& adj)
      {
        // get counts
        _num_nodes_image = adj.get_num_nodes_domain();
        Index num_indices_image = 0;
 
        // allocate pointer vector
        _domain_ptr.resize(adj.get_num_nodes_image() + 1, Index(0));
        // allocate auxiliary mask vector
        std::vector<char> vidx_mask(adj.get_num_nodes_image(), 0);
        char* idx_mask = vidx_mask.data();
 
        // loop over all image nodes
        for(Index j(0); j < adj.get_num_nodes_domain(); ++j)
        {
          for(auto it = adj.image_begin(j); it != adj.image_end(j); ++it)
          {
            if(idx_mask[*it] == 0)
            {
              ++num_indices_image;
              ++_domain_ptr[(*it)+1];
              idx_mask[*it] = 1;
            }
          }
          for(auto it = adj.image_begin(j); it != adj.image_end(j); ++it)
            idx_mask[*it] = 0;
        }
 
        _image_idx.resize(num_indices_image);
        std::vector<Index*> vimg_ptr(adj.get_num_nodes_image(), nullptr);
        Index** image_ptr = vimg_ptr.data();
        Index* image_idx = _image_idx.data();
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj.get_num_nodes_image()+1u, _domain_ptr.data(), _domain_ptr.data());
        for(Index i(0); i < adj.get_num_nodes_image(); ++i)
        {
          image_ptr[i] = &image_idx[_domain_ptr[i]];
        }
 
        // build image index vector
        for(Index j(0); j < adj.get_num_nodes_domain(); ++j)
        {
          for(auto it = adj.image_begin(j); it != adj.image_end(j); ++it)
          {
            if(idx_mask[*it] == 0)
            {
              Index*& idx = image_ptr[*it];
              *idx = j;
              ++idx;
              idx_mask[*it] = 1;
            }
          }
          for(auto it = adj.image_begin(j); it != adj.image_end(j); ++it)
            idx_mask[*it] = 0;
        }
      }
 
      template<
        typename Adjactor1_,
        typename Adjactor2_>
      void _render_as_is(
        const Adjactor1_& adj1,
        const Adjactor2_& adj2)
      {
         // validate adjactor dimensions
        XASSERTM(adj1.get_num_nodes_image() == adj2.get_num_nodes_domain(), "Adjactor dimension mismatch!");
 
        typedef typename Adjactor1_::ImageIterator AImIt1;
        typedef typename Adjactor2_::ImageIterator AImIt2;
 
        // get counts
        _num_nodes_image = adj2.get_num_nodes_image();
        Index num_indices_image = 0;
 
        // allocate pointer vector
        _domain_ptr.resize(adj1.get_num_nodes_domain() + 1);
        _domain_ptr[0] = 0;
        // count number of adjacencies and build pointer vector
        FEAT_PRAGMA_OMP(parallel for schedule(dynamic, 1000))
        for(Index i = 0; i < adj1.get_num_nodes_domain(); ++i)
        {
          Index num_indices_here = Index(0);
          AImIt1 cur1(adj1.image_begin(i));
          AImIt1 end1(adj1.image_end(i));
          for(; cur1 != end1; ++cur1)
          {
            AImIt2 cur2(adj2.image_begin(*cur1));
            AImIt2 end2(adj2.image_end(*cur1));
            for(; cur2 != end2; ++cur2)
            {
              ++num_indices_here;
            }
          }
          _domain_ptr[i+1] = num_indices_here;
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj1.get_num_nodes_domain()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj1.get_num_nodes_domain()];
 
        // allocate and build index vector
        _image_idx.resize(num_indices_image);
 
        FEAT_PRAGMA_OMP(parallel for schedule(dynamic, 1000))
        for(Index i = 0; i < adj1.get_num_nodes_domain(); ++i)
        {
          Index k = _domain_ptr[i];
          AImIt1 cur1(adj1.image_begin(i));
          AImIt1 end1(adj1.image_end(i));
          for(; cur1 != end1; ++cur1)
          {
            AImIt2 cur2(adj2.image_begin(*cur1));
            AImIt2 end2(adj2.image_end(*cur1));
            for(; cur2 != end2; ++cur2, ++k)
            {
              _image_idx[k] = *cur2;
            }
          }
        }
      }
 
      template<
        typename Adjactor1_,
        typename Adjactor2_>
      void _render_injectify(
        const Adjactor1_& adj1,
        const Adjactor2_& adj2)
      {
         // validate adjactor dimensions
        XASSERTM(adj1.get_num_nodes_image() == adj2.get_num_nodes_domain(), "Adjactor dimension mismatch!");
 
        // get counts
        _num_nodes_image = adj2.get_num_nodes_image();
        Index num_indices_image = 0;
 
        // allocate pointer vector
        _domain_ptr.resize(adj1.get_num_nodes_domain() + 1);
        _domain_ptr[0] = Index(0);
        FEAT_PRAGMA_OMP(parallel)
        {
          // allocate auxiliary mask vector
          std::vector<char> vidx_mask(adj2.get_num_nodes_image(), 0);
          char* idx_mask = vidx_mask.data();
 
          // count number of adjacencies and build pointer vector
          FEAT_PRAGMA_OMP(for schedule(dynamic, 1000))
          for(Index i=0; i < adj1.get_num_nodes_domain(); ++i)
          {
            Index num_indices_here = Index(0);
            for(auto it = adj1.image_begin(i); it != adj1.image_end(i); ++it)
            {
              for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
              {
                if(idx_mask[*jt] == 0)
                {
                  ++num_indices_here;
                  idx_mask[*jt] = 1;
                }
              }
            }
            _domain_ptr[i+1] = num_indices_here;
            // reset mask
            for(auto it = adj1.image_begin(i); it != adj1.image_end(i); ++it)
              for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
                idx_mask[*jt] = 0;
          }
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj1.get_num_nodes_domain()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj1.get_num_nodes_domain()];
        _image_idx.resize(num_indices_image);
 
        FEAT_PRAGMA_OMP(parallel)
        {
          std::vector<char> vidx_mask(adj2.get_num_nodes_image(), 0);
          char* idx_mask = vidx_mask.data();
          FEAT_PRAGMA_OMP(for schedule(dynamic, 1000))
          for(Index i = 0; i < adj1.get_num_nodes_domain(); ++i)
          {
            Index k = _domain_ptr[i];
            for(auto it = adj1.image_begin(i); it != adj1.image_end(i); ++it)
            {
              for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
              {
                if(idx_mask[*jt] == 0)
                {
                  _image_idx[k] = *jt;
                  ++k;
                  idx_mask[*jt] = 1;
                }
              }
            }
            // reset mask
            for(auto it = adj1.image_begin(i); it != adj1.image_end(i); ++it)
              for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
                idx_mask[*jt] = 0;
          }
        }
      }
 
      template<
        typename Adjactor1_,
        typename Adjactor2_>
      void _render_transpose(
        const Adjactor1_& adj1,
        const Adjactor2_& adj2)
      {
        // validate adjactor dimensions
        XASSERTM(adj1.get_num_nodes_image() == adj2.get_num_nodes_domain(), "Adjactor dimension mismatch!");
 
        typedef typename Adjactor1_::ImageIterator AImIt1;
        typedef typename Adjactor2_::ImageIterator AImIt2;
 
        // get counts
        _num_nodes_image = adj1.get_num_nodes_domain();
        Index num_indices_image = 0;
 
        // allocate and format pointer vector
        _domain_ptr.resize(adj2.get_num_nodes_image() + 1, Index(0));
 
        // count number of adjacencies
        for(Index j(0); j < adj1.get_num_nodes_domain(); ++j)
        {
          AImIt1 cur1(adj1.image_begin(j));
          AImIt1 end1(adj1.image_end(j));
          for(; cur1 != end1; ++cur1)
          {
            AImIt2 cur2(adj2.image_begin(*cur1));
            AImIt2 end2(adj2.image_end(*cur1));
            for(; cur2 != end2; ++cur2)
            {
              ++_domain_ptr[(*cur2) + 1];
            }
          }
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj2.get_num_nodes_image()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj2.get_num_nodes_image()];
 
        // allocate and build index vector
        _image_idx.resize(num_indices_image);
        std::vector<Index*> vimg_ptr(adj2.get_num_nodes_image(), nullptr);
        Index** image_ptr = vimg_ptr.data();
        Index* image_idx = _image_idx.data();
 
        for(Index i(0); i < adj2.get_num_nodes_image(); ++i)
        {
          image_ptr[i] = &image_idx[_domain_ptr[i]];
        }
 
        for(Index j(0); j < adj1.get_num_nodes_domain(); ++j)
        {
          AImIt1 cur1(adj1.image_begin(j));
          AImIt1 end1(adj1.image_end(j));
          for(; cur1 != end1; ++cur1)
          {
            AImIt2 cur2(adj2.image_begin(*cur1));
            AImIt2 end2(adj2.image_end(*cur1));
            for(; cur2 != end2; ++cur2)
            {
              Index*& idx = image_ptr[*cur2];
              *idx = j;
              ++idx;
            }
          }
        }
      }
 
      template<
        typename Adjactor1_,
        typename Adjactor2_>
      void _render_injectify_transpose(
        const Adjactor1_& adj1,
        const Adjactor2_& adj2)
      {
        // validate adjactor dimensions
        XASSERTM(adj1.get_num_nodes_image() == adj2.get_num_nodes_domain(), "Adjactor dimension mismatch!");
 
        // get counts
        _num_nodes_image = adj1.get_num_nodes_domain();
        Index num_indices_image = 0;
 
        // allocate pointer vector
        _domain_ptr.resize(adj2.get_num_nodes_image() + 1, Index(0));
 
        // allocate auxiliary mask vector
        std::vector<char> vidx_mask(adj2.get_num_nodes_image(), 0);
        char* idx_mask = vidx_mask.data();
 
        // loop over all image nodes
        for(Index j(0); j < adj1.get_num_nodes_domain(); ++j)
        {
          for(auto it = adj1.image_begin(j); it != adj1.image_end(j); ++it)
          {
            for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
            {
              if(idx_mask[*jt] == 0)
              {
                ++_domain_ptr[(*jt)+1];
                idx_mask[*jt] = 1;
              }
            }
          }
          // reset mask
          for(auto it = adj1.image_begin(j); it != adj1.image_end(j); ++it)
            for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
              idx_mask[*jt] = 0;
        }
 
        // perform inclusive scan to obtain domain pointer
        feat_omp_in_scan(adj2.get_num_nodes_image()+1u, _domain_ptr.data(), _domain_ptr.data());
        num_indices_image = _domain_ptr[adj2.get_num_nodes_image()];
 
        _image_idx.resize(num_indices_image);
        std::vector<Index*> vimg_ptr(adj2.get_num_nodes_image(), nullptr);
        Index** image_ptr = vimg_ptr.data();
        Index* image_idx = _image_idx.data();
 
        // build pointer vector
        for(Index i(0); i < adj2.get_num_nodes_image(); ++i)
        {
          image_ptr[i] = &image_idx[_domain_ptr[i]];
        }
 
        // build image index vector
        for(Index j(0); j < adj1.get_num_nodes_domain(); ++j)
        {
          for(auto it = adj1.image_begin(j); it != adj1.image_end(j); ++it)
          {
            for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
            {
              if(idx_mask[*jt] == 0)
              {
                Index*& idx = image_ptr[*jt];
                *idx = j;
                ++idx;
                idx_mask[*jt] = 1;
              }
            }
          }
          // reset mask
          for(auto it = adj1.image_begin(j); it != adj1.image_end(j); ++it)
            for(auto jt = adj2.image_begin(*it); jt != adj2.image_end(*it); ++jt)
              idx_mask[*jt] = 0;
        }
      }
 
      /* ******************************************************************* */
      /*  A D J A C T O R   I N T E R F A C E   I M P L E M E N T A T I O N  */
      /* ******************************************************************* */
    public:
 
      inline Index get_num_nodes_domain() const
      {
        return (_domain_ptr.empty() ? Index(0) : Index(_domain_ptr.size() - 1));
      }
 
      inline Index get_num_nodes_image() const
      {
        return _num_nodes_image;
      }
 
      inline ImageIterator image_begin(Index domain_node) const
      {
        ASSERTM(domain_node +1 < _domain_ptr.size(), "Domain node index out of range");
 
        return _image_idx.begin() + IndexVector::difference_type(_domain_ptr[domain_node]);
      }
 
      inline ImageIterator image_end(Index domain_node) const
      {
        ASSERTM(domain_node +1< _domain_ptr.size(), "Domain node index out of range");
 
        return _image_idx.begin() + IndexVector::difference_type(_domain_ptr[domain_node + Index(1)]);
      }
    }; // class Graph
  } // namespace Adjacency
} // namespace FEAT