communication_policies.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-2015:  G-CSC, Goethe University Frankfurt
 * Authors: Andreas Vogel, Sebastian Reiter
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __H__LIB_ALGEBRA__PARALLELIZATION__COMMUNICATION_POLICIES__
#define __H__LIB_ALGEBRA__PARALLELIZATION__COMMUNICATION_POLICIES__
 
#include "pcl/pcl.h"
#include "common/serialization.h"
#include "parallel_index_layout.h"
#include "algebra_id.h"
 
 
namespace ug{
 
// predeclaration of block_traits
template <typename t> struct block_traits
{
  enum{
    is_static = 1
  };
};
 
 
 
template <class TVector>
class ComPol_VecCopy : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecCopy() : m_pVecDest(NULL), m_pVecSrc(NULL) {}
 
    ComPol_VecCopy(TVector* pVec): m_pVecDest(pVec), m_pVecSrc(pVec)  {}
 
    ComPol_VecCopy(TVector* pVecDest, const TVector* pVecSrc)
      : m_pVecDest(pVecDest), m_pVecSrc(pVecSrc)  {}
 
    void set_vector(TVector* pVec)  {m_pVecDest = pVec; m_pVecSrc = pVec;}
 
    void set_vector(TVector* pVecDest, const TVector* pVecSrc)
      {m_pVecDest = pVecDest; m_pVecSrc = pVecSrc;}
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecCopy_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVecSrc == NULL) return false;
 
    //  rename for convenience
      const TVector& v = *m_pVecSrc;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  write entry into buffer
        Serialize(buff, v[index]);
      }
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecCopy_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVecDest == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVecDest;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  write entry into vector
        Deserialize(buff, v[index]);
      }
      return true;
    }
 
  private:
  //  pointer to current vector
    TVector* m_pVecDest;
    const TVector* m_pVecSrc;
};
 
 
template <class TVector>
class ComPol_VecScaleCopy : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecScaleCopy() : m_pVec(NULL), m_scale(1.0)  {}
 
    ComPol_VecScaleCopy(TVector* pVec, number scale)
      : m_pVec(pVec), m_scale(scale)  {}
 
    void set_vector(TVector* pVec)  {m_pVec = pVec;}
 
    void set_scale(number scale)  {m_scale = scale;}
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecScaleCopy_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy value into buffer
        Serialize(buff, v[index]);
      }
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecScaleCopy_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy value from buffer
        Deserialize(buff, v[index]);
 
      //  scale value
        v[index] *= m_scale;
      }
      return true;
    }
 
  private:
    TVector* m_pVec;
 
    number m_scale;
};
 
 
template <class TVector>
class ComPol_VecAdd : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecAdd() : m_pVecDest(NULL), m_pVecSrc(NULL) {}
 
    ComPol_VecAdd(TVector* pVec) : m_pVecDest(pVec),  m_pVecSrc(pVec) {}
 
    ComPol_VecAdd(TVector* pVecDest, const TVector* pVecSrc)
      : m_pVecDest(pVecDest),  m_pVecSrc(pVecSrc) {}
 
    void set_vector(TVector* pVec)  {m_pVecDest = pVec; m_pVecSrc = pVec;}
 
    void set_vector(TVector* pVecDest, const TVector* pVecSrc)
      {m_pVecDest = pVecDest; m_pVecSrc = pVecSrc;}
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecAdd_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVecSrc == NULL) return false;
 
    //  rename for convenience
      const TVector& v = *m_pVecSrc;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      // copy entry
        Serialize(buff, v[index]);
      }
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecAdd_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVecDest == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVecDest;
 
    //  entry
      typename TVector::value_type entry;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy entry
        Deserialize(buff, entry);
 
      //  add entry
        v[index] += entry;
      }
      return true;
    }
 
  private:
    TVector* m_pVecDest;
    const TVector* m_pVecSrc;
};
 
 
template <class TVector>
class ComPol_VecScaleAdd : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecScaleAdd() : m_pVec(NULL), m_scale(1.0) {}
 
    ComPol_VecScaleAdd(TVector* pVec) : m_pVec(pVec)  {}
 
    void set_vector(TVector* pVec)  {m_pVec = pVec;}
 
    void set_scale(number scale)  {m_scale = scale;}
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecScaleAdd_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      // copy entry
        Serialize(buff, v[index]);
      }
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecScaleAdd_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  entry
      typename TVector::value_type entry;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy entry
        Deserialize(buff, entry);
 
      //  add entry
        v[index] += entry * m_scale;
      }
      return true;
    }
 
  private:
    TVector* m_pVec;
 
    number m_scale;
};
 
 
 
template <class TVector>
class ComPol_VecAddSetZero : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecAddSetZero() : m_pVec(0)  {}
 
    ComPol_VecAddSetZero(TVector* pVec) : m_pVec(pVec)  {}
 
    void set_vector(TVector* pVec)  {m_pVec = pVec;}
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecAddSetZero_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      // copy values
        Serialize(buff, v[index]);
        v[index] *= 0;
      }
 
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecAddSetZero_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  entry
      typename TVector::value_type entry;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy values
        Deserialize(buff, entry);
 
      //  add entry
        v[index] += entry;
      }
      return true;
    }
 
  private:
    TVector* m_pVec;
};
 
 
template <class TVector>
class ComPol_VecSubtract : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecSubtract() : m_pVec(NULL) {}
 
    ComPol_VecSubtract(TVector* pVec) : m_pVec(pVec)  {}
 
    void set_vector(TVector* pVec)  {m_pVec = pVec;}
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecSubtract_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      // get index
        const size_t index = interface.get_element(iter);
 
      // copy value
        Serialize(buff, v[index]);
      }
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecSubtract_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    // entry
      typename TVector::value_type entry;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy vector
        Deserialize(buff, entry);
 
      //  subtract entry
        v[index] -= entry;
      }
      return true;
    }
 
  private:
    TVector* m_pVec;
};
 
template <class TVector>
class ComPol_CheckConsistency : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_CheckConsistency() : m_pVec(NULL)  {}
 
    ComPol_CheckConsistency(const TVector* pVec) : m_pVec(pVec) {}
 
    void set_vector(const TVector* pVec)  {m_pVec = pVec;}
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_CheckConsistency_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      const TVector& v = *m_pVec;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      // get index
        const size_t index = interface.get_element(iter);
 
      // copy value
        Serialize(buff, v[index]);
      }
      return true;
    }
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_CheckConsistency_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      const TVector& v = *m_pVec;
 
    // entry
      typename TVector::value_type entry, diff;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy vector
        Deserialize(buff, entry);
 
        diff = v[index];
        diff -= entry;
 
      //  check
        if(VecNormSquared(diff) > 1e-20 * VecNormSquared(entry)){
          UG_LOG_ALL_PROCS("Slave and master value at index "<<index<<
                           " differ by: "<<VecNormSquared(diff)<<"\n");
        }
      }
 
      return true;
    }
 
  private:
    const TVector* m_pVec;
};
 
 
template <class TVector>
class ComPol_VecSubtractOnlyOneSlave : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_VecSubtractOnlyOneSlave() : m_pVec(NULL) {}
 
    ComPol_VecSubtractOnlyOneSlave(TVector* pVec) 
      {
      set_vector(pVec);
    }
 
    void set_vector(TVector* pVec) 
    {
      m_pVec = pVec;
      m_vProcessed.resize(m_pVec->size(), false);
    }
 
    void clear()
    {
      m_vProcessed.clear();
      m_vProcessed.resize(m_pVec->size(), false);
    }
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
    {
      if(block_traits<typename TVector::value_type>::is_static)
        return interface.size() * sizeof(typename TVector::value_type);
      else
        return -1;
    }
 
 
    virtual bool
    collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
 
      PROFILE_BEGIN_GROUP(ComPol_VecSubtractOnlyOneSlave_collect, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      // get index
        const size_t index = interface.get_element(iter);
 
      // copy value
        Serialize(buff, v[index]);
      }
      return true;
    }
 
 
    virtual bool
    extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_VecSubtractOnlyOneSlave_extract, "algebra parallelization");
    //  check that vector has been set
      if(m_pVec == NULL) return false;
 
    //  rename for convenience
      TVector& v = *m_pVec;
 
    // entry
      typename TVector::value_type entry;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  copy vector
        Deserialize(buff, entry);
 
      //  subtract entry
        if(m_vProcessed[index] == false)
        {
          v[index] -= entry;
          m_vProcessed[index] = true;
        }
      }
      return true;
    }
 
  private:
    TVector* m_pVec;
    std::vector<bool> m_vProcessed;
};
 
 
 
 
template <class TAlgebra>
class ComPol_MatDistributeDiag
  : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
 
  typedef typename TAlgebra::matrix_type TMatrix;
  typedef typename TAlgebra::vector_type TVector;
 
 
  ComPol_MatDistributeDiag(TMatrix& rMat, TVector &rWeight, double theta=1.0)
      : m_rMat(rMat), m_rWeight(rWeight), m_dTheta(theta)
    {
      UG_ASSERT( m_rMat.num_rows()==m_rWeight.size(), "sizes should match");
    }
 
 
    virtual int
    get_required_buffer_size(const Interface& interface)
      {
        if(block_traits<typename TMatrix::value_type>::is_static)
          return interface.size() * sizeof(typename TMatrix::value_type);
        else
          return -1;
      }
 
 
    virtual bool collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(CMatDistributeDiag_collect, "algebra parallelization");
      typedef typename TMatrix::value_type block_type;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
        block_type& a_ii = m_rMat(index, index);
 
      // send rescaled diagonal entry
        block_type out_ii = a_ii;
        out_ii = out_ii * (1.0/BlockRef(m_rWeight[index], 0));
        //  write diagonal entry to stream
        Serialize(buff, out_ii);
 
      //  std::cerr << "Sending(" << index << ")=" << out_ii <<  " "<< BlockRef(m_rWeight[index], 0) << std::endl;
      // store average
        a_ii = out_ii;
 
      // after the first call, we need to reset the weight (on master)
      // in order to communicate correct entries in subsequent if calls
        m_rWeight[index] = 1.0;
 
      }
 
      return true;
    }
 
    virtual bool
    begin_layout_extraction(const Layout* pLayout)
    { return true; }
 
    virtual bool extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(MatDistributeDiag_extract, "algebra parallelization");
    //  block type of associated matrix
      typedef typename TMatrix::value_type block_type;
 
      block_type in_ii;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
        //  get index
        const size_t index = interface.get_element(iter);
 
        //  read incoming diagonal
        Deserialize(buff, in_ii);
 
        // set matrix entry
      //  std::cerr << "Recv'd(" << index << ")=" << in_ii << std::endl;
        m_rMat(index, index) += in_ii;
      }
 
      return true;
    }
 
  private:
  //  pointer to current vector
    TMatrix& m_rMat;
    TVector& m_rWeight;
    double m_dTheta;
 
 
};
 
 
template <class TMatrix>
class ComPol_MatAddRowsOverlap0
  : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
 
    ComPol_MatAddRowsOverlap0(TMatrix& rMat, AlgebraIDVec& vGlobalID)
      : m_rMat(rMat), m_vGlobalID(vGlobalID)
    {
      UG_ASSERT(vGlobalID.size() >= m_rMat.num_rows(), "too few Global ids");
    }
 
    virtual bool collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_collect, "algebra parallelization");
      typedef typename TMatrix::row_iterator row_iterator;
      typedef typename TMatrix::value_type block_type;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  count number of row entries
        const row_iterator rowEnd = m_rMat.end_row(index);
        size_t numRowEntry = 0;
        for(row_iterator it_k = m_rMat.begin_row(index); it_k != rowEnd; ++it_k)
          numRowEntry++;
 
      //  write number of row entries to stream
        Serialize(buff, numRowEntry);
 
      //  write entries and global id to stream
        for(row_iterator it_k = m_rMat.begin_row(index); it_k != rowEnd; ++it_k)
        {
          const size_t k = it_k.index();
          block_type& a_ik = it_k.value();
 
        //  write global entry to buffer
          Serialize(buff, m_vGlobalID[k]);
 
        //  write entry into buffer
          Serialize(buff, a_ik);
        }
      }
 
      return true;
    }
 
    virtual bool
    begin_layout_extraction(const Layout* pLayout)
    {
    //  fill the map global->local
      GenerateAlgebraIDHashList(m_algIDHash, m_vGlobalID);
      return true;
    }
 
    virtual bool extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_extract, "algebra parallelization");
    //  block type of associated matrix
      typedef typename TMatrix::value_type block_type;
 
    //  we'll read global ids into this variable
      AlgebraID gID;
 
    //  we'll read blocks into this var
      block_type block;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  read the number of connections
        size_t numConnections = 0;
        Deserialize(buff, numConnections);
 
      //  read each connection
        for(size_t i_conn = 0; i_conn < numConnections; ++i_conn){
          Deserialize(buff, gID);
          Deserialize(buff, block);
 
        //  if gID exists on this process, then add the connection to
        //  the matrix.
          size_t conInd;
          if(m_algIDHash.get_entry(conInd, gID)){
          //  add connection between index and conInd to matrix
            m_rMat(index, conInd) += block;
          }
        }
      }
 
      return true;
    }
 
  private:
  //  pointer to current vector
    TMatrix& m_rMat;
 
  //  map localID->globalID
    AlgebraIDVec& m_vGlobalID;
 
  //  map globalID->localID
    AlgebraIDHashList m_algIDHash;
 
};
 
 
 
template <class TMatrix>
class ComPol_MatCopyRowsOverlap0
  : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
 
    ComPol_MatCopyRowsOverlap0(TMatrix& rMat, AlgebraIDVec& vGlobalID)
      : m_rMat(rMat), m_vGlobalID(vGlobalID)
    {
      UG_ASSERT(vGlobalID.size() >= m_rMat.num_rows(), "too few Global ids");
    }
 
    virtual bool collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_collect, "algebra parallelization");
      typedef typename TMatrix::row_iterator row_iterator;
      typedef typename TMatrix::value_type block_type;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  count number of row entries
        const row_iterator rowEnd = m_rMat.end_row(index);
        size_t numRowEntry = 0;
        for(row_iterator it_k = m_rMat.begin_row(index); it_k != rowEnd; ++it_k)
          numRowEntry++;
 
      //  write number of row entries to stream
        Serialize(buff, numRowEntry);
 
      //  write entries and global id to stream
        for(row_iterator it_k = m_rMat.begin_row(index); it_k != rowEnd; ++it_k)
        {
          const size_t k = it_k.index();
          block_type& a_ik = it_k.value();
 
        //  write global entry to buffer
          Serialize(buff, m_vGlobalID[k]);
 
        //  write entry into buffer
          Serialize(buff, a_ik);
        }
      }
 
      return true;
    }
 
    virtual bool
    begin_layout_extraction(const Layout* pLayout)
    {
    //  fill the map global->local
      GenerateAlgebraIDHashList(m_algIDHash, m_vGlobalID);
      return true;
    }
 
    virtual bool extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_extract, "algebra parallelization");
    //  block type of associated matrix
      typedef typename TMatrix::value_type block_type;
 
    //  we'll read global ids into this variable
      AlgebraID gID;
 
    //  we'll read blocks into this var
      block_type block;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  read the number of connections
        size_t numConnections = 0;
        Deserialize(buff, numConnections);
 
      //  read each connection
        for(size_t i_conn = 0; i_conn < numConnections; ++i_conn){
          Deserialize(buff, gID);
          Deserialize(buff, block);
 
        //  if gID exists on this process, then set the connection to
        //  the new value.
          size_t conInd;
          if(m_algIDHash.get_entry(conInd, gID)){
            m_rMat(index, conInd) = block;
          }
        }
      }
 
      return true;
    }
 
  private:
  //  pointer to current vector
    TMatrix& m_rMat;
 
  //  map localID->globalID
    AlgebraIDVec& m_vGlobalID;
 
  //  map globalID->localID
    AlgebraIDHashList m_algIDHash;
 
};
 
 
 
template <class TMatrix>
class ComPol_MatAddSetZeroInnerInterfaceCouplings
  : public pcl::ICommunicationPolicy<IndexLayout>
{
  public:
    ComPol_MatAddSetZeroInnerInterfaceCouplings(TMatrix& rMat)
      : m_rMat(rMat)
    {}
 
    virtual bool collect(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_MatAddInnerInterfaceCouplings_collect, "algebra parallelization");
      typedef typename TMatrix::row_iterator row_iterator;
      typedef typename TMatrix::value_type block_type;
 
    //  build map
      Hash<size_t, size_t> hash((size_t)(interface.size() * 1.1));
      size_t interfaceIndex = 0;
 
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter, ++interfaceIndex){
        hash.insert(interface.get_element(iter), interfaceIndex);
      }
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  count number of row entries
        const row_iterator rowEnd = m_rMat.end_row(index);
        size_t numRowEntry = 0;
        for(row_iterator it_k = m_rMat.begin_row(index); it_k != rowEnd; ++it_k){
          if(!hash.has_entry(it_k.index())) continue;
          numRowEntry++;
        }
 
      //  write number of row entries to stream
        Serialize(buff, numRowEntry);
 
      //  write entries and global id to stream
        for(row_iterator it_k = m_rMat.begin_row(index); it_k != rowEnd; ++it_k)
        {
          const size_t target = it_k.index();
          if(!hash.get_entry(interfaceIndex, target)) continue;
          block_type& a_ik = it_k.value();
 
        //  write global entry to buffer
          Serialize(buff, interfaceIndex);
 
        //  write entry into buffer
          Serialize(buff, a_ik);
 
        //  set entry to zero
          a_ik *= 0;
        }
      }
 
      return true;
    }
 
    virtual bool extract(ug::BinaryBuffer& buff, const Interface& interface)
    {
      PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_extract, "algebra parallelization");
    //  block type of associated matrix
      typedef typename TMatrix::value_type block_type;
 
    //  build map
      std::vector<size_t> vMapInterfaceToGlob(interface.size());
      int k = 0;
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter, ++k){
        vMapInterfaceToGlob[k] = interface.get_element(iter);
      }
 
    //  we'll read blocks into this var
      block_type block;
      size_t locID;
 
    //  loop interface
      for(typename Interface::const_iterator iter = interface.begin();
        iter != interface.end(); ++iter)
      {
      //  get index
        const size_t index = interface.get_element(iter);
 
      //  read the number of connections
        size_t numConnections = 0;
        Deserialize(buff, numConnections);
 
      //  read each connection
        for(size_t i_conn = 0; i_conn < numConnections; ++i_conn){
          Deserialize(buff, locID);
          Deserialize(buff, block);
 
          UG_ASSERT(locID < vMapInterfaceToGlob.size(), "Invalid"
                    " id: "<<locID<<", size: "<<vMapInterfaceToGlob.size());
          m_rMat(index, vMapInterfaceToGlob[locID]) += block;
        }
      }
 
      return true;
    }
 
  private:
  //  pointer to current vector
    TMatrix& m_rMat;
};
 
 
 
 
}// end of namespace
 
 
#endif /* __H__LIB_ALGEBRA__PARALLELIZATION__COMMUNICATION_POLICIES__ */