numaMatrix.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2013-2015 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef NUMAMATRIX_HH
#define NUMAMATRIX_HH
 
#include <algorithm>
#include <cassert>
#include <vector>
 
#include <boost/iterator/iterator_facade.hpp>
 
#include "linalg/dynamicMatrix.hh"
#include "utilities/threading.hh"
 
namespace Kaskade {
  
  // forward declaration
  template <class Entry> class NumaDenseMatrix;
  
  namespace NumaDenseMatrixDetail
  {
    
    template <class Entry>
    class ColumnIterator: public boost::iterator_facade<ColumnIterator<Entry>,Entry,boost::random_access_traversal_tag>
    {
    public:
      ColumnIterator() = default;
      
      ColumnIterator(Entry* base_, size_t stride_, size_t idx_=0)
      : base(base_), stride(stride_), idx(idx_) {}
      
      size_t index() const { return idx; }
      
    private:
      friend class boost::iterator_core_access;
      
      bool equal(ColumnIterator<Entry> const& p) const
      {
        return idx==p.idx;
      }
      
      void increment() { ++idx; }
      void decrement() { --idx; }
      void advance(typename ColumnIterator<Entry>::difference_type n) { idx += n; }
      template <class E>
      typename ColumnIterator<Entry>::difference_type distance_to(ColumnIterator<E> const& i) const { return i.idx - idx; }
      
      
      Entry& dereference() const { return base[idx*stride]; }
      
      Entry* base;
      size_t stride;
      size_t idx;
    };
    
    // ----------------------------------------------------------------------------------------------
    
    template <class Entry>
    class Row 
    {
    public:
      typedef ColumnIterator<Entry> iterator;
      typedef ColumnIterator<Entry const> const_iterator;
      
      Row() = default;
      
      Row(Entry* first_, size_t stride_, std::size_t cols_)
      : first(first_), stride(stride_), cols(cols_) {}
      
      iterator begin() { return iterator(first,stride,0); }
      iterator end()   { return iterator(first,stride,cols); }
      const_iterator begin() const { return const_iterator(first,stride,0); }
      const_iterator end()   const { return const_iterator(first,stride,cols); }
      
      Entry&       operator[](size_t c)       { return first[stride*c]; }
      Entry const& operator[](size_t c) const { return first[stride*c]; }
      
    private:
      Entry* first;
      size_t stride;
      size_t cols;
    };
    
    // ----------------------------------------------------------------------------------------------
 
    template <class Entry>
    class RowIterator: public boost::iterator_facade<RowIterator<Entry>,
                                                     typename std::conditional<std::is_const<Entry>::value,
                                                                               Row<typename std::remove_const<Entry>::type> const,
                                                                               Row<Entry>>::type,
                                                     boost::random_access_traversal_tag>
    {
      typedef typename std::remove_const<Entry>::type NonConstEntry;
      
      typedef typename std::conditional<std::is_const<Entry>::value,
                                        NumaDenseMatrix<NonConstEntry> const,
                                        NumaDenseMatrix<Entry>>::type Matrix;
      
      typedef boost::iterator_facade<RowIterator<Entry>,
                                    typename std::conditional<std::is_const<Entry>::value,
                                                              Row<NonConstEntry> const,
                                                              Row<Entry>>::type,
                                    boost::random_access_traversal_tag> Base;
      
      
    public:
      RowIterator() = default;
      
      RowIterator(Matrix& mat_, size_t idx_)
      : mat(&mat_) 
      {
        update(idx_);
      }
      
      size_t index() const { return idx; }
      
    private:
      friend class boost::iterator_core_access;
      using DifferenceType = Base::difference_type;
      
      bool equal(RowIterator<Entry> const& other) const { return idx==other.idx;}
      void increment() { update(idx+1); }
      void decrement() { update(idx-1); }
      void advance(typename RowIterator<Entry>::difference_type n) { update(idx+n); }
      template <class E> DifferenceType distance_to(RowIterator<E> const& i) const { return static_cast<DifferenceType>(i.idx) - idx; }
      typename Base::reference dereference() const { return row; } 
      
      // moves row as to represent row with index newidx.
      void update(size_t newidx)
      {
        idx = newidx;
        if (idx<mat->N()) // to not trip over the edge when arriving at end()
        {
          size_t chunk = uniformWeightRange(idx,mat->chunks.size(),mat->N());
          size_t rowStart = uniformWeightRangeStart(chunk,mat->chunks.size(),mat->N());
          size_t stride = uniformWeightRangeStart(chunk+1,mat->chunks.size(),mat->N()) - rowStart; // that's "LDA"
          row = typename Base::value_type(const_cast<NonConstEntry*>(&mat->chunks[chunk][idx-rowStart]),stride,mat->M());
        }
      }
      
      Matrix* mat;
      size_t idx; // row index
      mutable typename Base::value_type row; // dereference is const method, but can "point to" mutable row value type..
    };
 
    // ----------------------------------------------------------------------------------------------
 
    template <class T, int n>
    typename Dune::FieldTraits<T>::field_type frobenius_product(Dune::FieldVector<T,n> const& x, Dune::FieldVector<T,n> const& y) 
    { 
      typename Dune::FieldTraits<T>::field_type result = 0;
      for (int i=0; i<n; ++i)
        result += frobenius_product(x[i],y[i]);
      return result;
    }
    
    template <class T, int n, int m>
    typename Dune::FieldTraits<T>::field_type frobenius_product(Dune::FieldMatrix<T,n,m> const& x, Dune::FieldMatrix<T,n,m> const& y) 
    { 
      typename Dune::FieldTraits<T>::field_type result = 0;
      for (int i=0; i<n; ++i)
        for (int j=0; j<m; ++j)
          result += frobenius_product(x[i][j],y[i][j]);
      return result;
    }
    
    double frobenius_product(double x, double y) { return x*y; }
    float  frobenius_product(float x, float y)  { return x*y; }
    
    
    
    template <class T, int n>
    typename Dune::FieldTraits<typename Dune::FieldTraits<T>::field_type>::real_type frobenius_norm2(Dune::FieldVector<T,n> const& x) { return x.two_norm2(); }
    
    template <class T, int n, int m>
    typename Dune::FieldTraits<typename Dune::FieldTraits<T>::field_type>::real_type frobenius_norm2(Dune::FieldMatrix<T,n,m> const& x) { return x.frobenius_norm2(); }
    
    double frobenius_norm2(double x) { return x*x; }
    float  frobenius_norm2(float x)  { return x*x; }
    
    
    template <class T, int n>
    typename Dune::FieldTraits<typename Dune::FieldTraits<T>::field_type>::real_type one_norm(Dune::FieldVector<T,n> const& x) { return x.one_norm(); }
    
    template <class T, int n, int m>
    typename Dune::FieldTraits<typename Dune::FieldTraits<T>::field_type>::real_type one_norm(Dune::FieldMatrix<T,n,m> const& x) { return x.one_norm(); }
    
    double one_norm(double x) { return std::abs(x); }
    float  one_norm(float x)  { return std::abs(x); }
    
    template <class T, int n>
    void axpy(Dune::FieldVector<T,n>& y, typename Dune::FieldTraits<T>::field_type const& a, Dune::FieldVector<T,n> const& x) { y.axpy(a,x); }
    
    template <class T, int n, int m>
    void axpy(Dune::FieldMatrix<T,n,m>& y, typename Dune::FieldTraits<T>::field_type const& a, Dune::FieldMatrix<T,n,m> const& x) { y += a*x; }
    
    void axpy(double& y, double a, double x) { y += a*x; }
    void axpy(float& y,  float a,  float x) { y += a*x; }
    
    // ----------------------------------------------------------------------------------------------
    
    // forward declaration
    template <class Entry> class VectorIterator;
 
    template <class Entry>
    class NumaDenseBase
    {
    protected:
      typedef std::vector<Entry,NumaAllocator<Entry>> Chunk;  // the actual data storage type
      typedef NumaDenseBase<Entry> Self;                      
      
    public:
      typedef size_t      size_type;
      typedef Entry       value_type;
      
      typedef typename Dune::FieldTraits<Entry>::field_type     field_type;
      typedef typename Dune::FieldTraits<field_type>::real_type real_type;
      
      typedef Entry block_type;
      
      NumaDenseBase (): NumaDenseBase(0,0)
      {
      } 
      
      NumaDenseBase (NumaDenseBase const& other)
      : rows(other.rows), cols(other.cols), chunks(other.chunks.size())
      {
        parallelForNodes([&](int i, int n) // perform data copy on the respective
        {                                  // NUMA node
          chunks[i] = other.chunks[i];
        });
      }
      
      NumaDenseBase(NumaDenseBase&& a) = default;
      
      NumaDenseBase (size_type r, size_type c, value_type v = value_type()) :
      rows(r), cols(c)
      {
        size_type n = NumaThreadPool::instance().nodes();
        chunks.reserve(n);
        for (size_type i=0; i<n; ++i)
          chunks.emplace_back((uniformWeightRangeStart(i+1,n,r)-uniformWeightRangeStart(i,n,r))*c,v,NumaAllocator<Entry>(i));
      }
      
      NumaDenseBase& operator=(NumaDenseBase const& other)
      {
        assert(chunks.size() == other.chunks.size()); // this can only be violated if the number of NUMA nodes changes...
        rows = other.rows;
        cols = other.cols;
        parallelForNodes([&](int i, int n) // perform data copy on the respective
        {                                  // NUMA node
          chunks[i] = other.chunks[i];
        });
      }
      
      Self& operator=(field_type const& a)
      {
        for_each([&](Entry& e) { e = a; }); 
      }
      
 
    protected:
      void resize (size_type r, size_type c)
      {
        // Check whether we need to change anything. Note that due to row blocking,
        // the check rows*cols==n*m is NOT sufficient (the distribution of rows on 
        // NUMA nodes may change nevertheless).
        if (rows==r && cols==c)
          return;
        
        rows = r;
        cols = c;
        for (int i=0; i<chunks.size(); ++i) // TODO: parallelize
          chunks[i].resize((uniformWeightRangeStart(i+1,chunks.size(),r)-uniformWeightRangeStart(i,chunks.size(),r))*c);
      }
 
    public:
      size_type N() const { return rows; }
      
    protected:
      size_type M() const { return cols; }
      
      
      field_type dot(Self const& other) const
      {
        assert(rows==other.rows && cols==other.cols);
        std::vector<field_type> dots(chunks.size());
        parallelForNodes([&](int i, int n) 
        {
          assert(n==chunks.size());
          field_type sum = 0;
          for (size_t k=0; k<chunks[i].size(); ++k)
            sum += NumaDenseMatrixDetail::frobenius_product(chunks[i][k],other.chunks[i][k]);
          dots[i] = sum;
        });
        return std::accumulate(dots.begin(),dots.end(),0);
      }
      
      template <class Result, class F>
      Result accumulate(F f, Result const& init) const
      {
        std::vector<Result> rs(chunks.size());
        parallelForNodes([&](int i, int n) 
        {
          assert(n==this->chunks.size());
          Result r = init;
          for (auto const& e: chunks[i])
            r = f(r,e);
          rs[i] = r;
        });
        return std::accumulate(rs.begin()+1,rs.end(),rs[0],f);
      }
      
      template <class F>
      void for_each(F const& f) 
      {
        parallelForNodes([&](int i, int n) 
        {
          assert(n==this->chunks.size());
          std::for_each(begin(chunks[i]),end(chunks[i]),f);
        });
      }
      
      template <class F>
      void for_each2(F const& f, Self const& other) 
      {
        assert(rows==other.rows && cols==other.cols);
        parallelForNodes([&](int i, int n) 
        {
          assert(n==this->chunks.size());
          for (size_t j=0; j<chunks[i].size(); ++j)
            f(chunks[i][j],other.chunks[i][j]);
        });
      }
      
    private:
      friend class NumaDenseMatrixDetail::RowIterator<Entry>;
      friend class NumaDenseMatrixDetail::RowIterator<Entry const>;
      friend class NumaDenseMatrixDetail::VectorIterator<Entry>;
      friend class NumaDenseMatrixDetail::VectorIterator<Entry const>;
      
      size_type rows, cols;      // the matrix dimensions
      std::vector<Chunk> chunks; // for each NUMA node a separate array to store the entries
    };
    
  }
  //-------------------------------------------------------------------------------------------
  
  template <class Entry>
  class NumaDenseMatrix: public NumaDenseMatrixDetail::NumaDenseBase<Entry>
  {
    typedef NumaDenseMatrix<Entry> Self;
    typedef NumaDenseMatrixDetail::NumaDenseBase<Entry> Base;
    
  public:
    typedef size_t      size_type;
    typedef Entry       value_type;
    
    typedef typename Dune::FieldTraits<Entry>::field_type         field_type;
    typedef typename Dune::FieldTraits<field_type>::real_type     real_type;
    
    typedef NumaDenseMatrixDetail::RowIterator<Entry>       RowIterator;
    typedef NumaDenseMatrixDetail::RowIterator<Entry const> ConstRowIterator;
 
    typedef typename RowIterator::value_type      row_type;
    typedef typename ConstRowIterator::value_type const_row_type;
    
    NumaDenseMatrix () = default;
 
    NumaDenseMatrix (NumaDenseMatrix const& a) = default;
 
    NumaDenseMatrix (NumaDenseMatrix&& a) = default;
 
    NumaDenseMatrix (size_type r, size_type c, value_type v = value_type()) 
    : Base(r,c,v)
    {}
    
    NumaDenseMatrix& operator=(NumaDenseMatrix const& a) = default;
    
    Self& operator=(field_type const& a)
    {
      Base::operator=(a);
      return *this;
    }
     
    void resize (size_type r, size_type c)
    {
      Base::resize(r,c);
    }
 
    void setSize (size_type r, size_type c)
    {
      resize(r,c);
    }
    
    
    RowIterator begin() { return RowIterator(*this,0); }
    
    RowIterator end() { return RowIterator(*this,this->N()); }
    
    ConstRowIterator begin() const { return ConstRowIterator(*this,0); }
    
    ConstRowIterator end() const { return ConstRowIterator(*this,this->N()); }
    
    row_type  operator[] (size_type row) { return *RowIterator(*this,row); }
    
    const_row_type const operator[] (size_type row) const { return *ConstRowIterator(*this,row); }
    
    
    size_type M() const { return Base::M(); }
    
    Self& operator*= (field_type const& a)
    {
      this->for_each([&](Entry& e) { e *= a; });
      return *this;
    }
    
    Self& operator/= (field_type const& scalar)
    {
      return (*this) *= (1/scalar);
    }
    
    Self& operator+= (Self const& b)
    {
      this->for_each2([](Entry& y, Entry const& x) { y += x; },b);
      return *this;
    }
   
    Self& operator-= (Self const& b)
    {
      this->for_each2([](Entry& y, Entry const& x) { y -= x; },b);
      return *this;
    }
    
    
    real_type frobenius_norm () const { return std::sqrt(frobenius_norm2()); }
    
    real_type frobenius_norm2 () const
    {
      return this->accumulate([](real_type s, Entry const& e) { return s+NumaDenseMatrixDetail::frobenius_norm2(e); }, real_type(0)); 
    }
  };
  
  template <class Entry>
  std::ostream& operator<<(std::ostream& out, NumaDenseMatrix<Entry> const& mat)
  {
    for (auto ri=mat.begin(); ri!=mat.end(); ++ri)
    {
      std::copy(ri->begin(),ri->end(),std::ostream_iterator<Entry>(out," "));
      out << '\n';
    }
    return out;
  }
  
  // --------------------------------------------------------------------------
  
  // forward declaration
  template <class Entry> class NumaVector;
  
  namespace NumaDenseMatrixDetail
  {
    
    template <class Entry>
    class VectorIterator: public boost::iterator_facade<VectorIterator<Entry>,Entry,boost::random_access_traversal_tag>
    {
      typedef typename std::remove_const<Entry>::type NonConstEntry;
      
      typedef typename std::conditional<std::is_const<Entry>::value,
                                        NumaVector<NonConstEntry> const,
                                        NumaVector<Entry>>::type Vector;
    public:
      VectorIterator() = default;
      
      VectorIterator(Vector& vec_, typename Vector::size_type idx_)
      : vec(&vec_), idx(idx_) {}
      
      size_t index() const { return idx; }
      
    private:
      friend class boost::iterator_core_access;
      
      bool equal(VectorIterator<Entry> const& p) const
      {
        return idx==p.idx;
      }
      
      void increment() { ++idx; }
      void decrement() { --idx; }
      void advance(typename VectorIterator<Entry>::difference_type n) { idx += n; }
      template <class E>
      typename VectorIterator<Entry>::difference_type distance_to(VectorIterator<E> const& i) const { return i.idx - idx; }
      Entry& dereference() const 
      { 
        size_t chunk = uniformWeightRange(idx,vec->chunks.size(),vec->N());
        size_t chunkStart = uniformWeightRangeStart(chunk,vec->chunks.size(),vec->N());
        return vec->chunks[chunk][idx-chunkStart];
      }
      
      Vector* vec;
      typename Vector::size_type idx;
    };
  }
  template <class Entry>
  class NumaVector: public NumaDenseMatrixDetail::NumaDenseBase<Entry>
  {
    typedef NumaVector<Entry> Self;
    typedef NumaDenseMatrixDetail::NumaDenseBase<Entry> Base;
    
   public:
    typedef typename Dune::FieldTraits<Entry>::field_type         field_type;
    typedef typename Dune::FieldTraits<field_type>::real_type     real_type;
    typedef typename Base::size_type                              size_type;
    typedef typename Base::value_type                             value_type;
    typedef NumaDenseMatrixDetail::VectorIterator<Entry>          iterator;
    typedef NumaDenseMatrixDetail::VectorIterator<Entry const>    const_iterator;
    
 
    NumaVector(): Base(0,1) 
    {} 
 
    NumaVector(NumaVector const& a) = default;
 
    NumaVector (NumaVector&& a) = default;
 
    NumaVector (size_type r, value_type v = value_type()) :
      Base(r,1,v)
    {}
    
    NumaVector& operator=(NumaVector const& a) = default;
    
    void resize (size_type r) { Base::resize(r,1); }
    
    
    iterator begin() { return iterator(*this,0); }
    
    iterator end() { return iterator(*this,this->N()); }
    
    const_iterator begin() const { return const_iterator(*this,0); }
    
    const_iterator end() const { return const_iterator(*this,this->N()); }
    
    Entry&  operator[] (size_type row) { return *iterator(*this,row); }
    
    Entry const& operator[] (size_type row) const { return *const_iterator(*this,row); }
    
     
    Self& operator*= (field_type const& a)
    {
      this->for_each([&](Entry& e) { e *= a; });
      return *this;
    }
    
    Self& operator/= (field_type const& a)
    {
      return (*this) *= (1/a);
    }
    
    Self& operator+= (Self const& x)
    {
      this->for_each2([](Entry& yi, Entry const& xi) { yi += xi; },x);
      return *this;
    }
   
    Self& operator-= (Self const& x)
    {
      this->for_each2([](Entry& yi, Entry const& xi) { yi -= xi; },x);
      return *this;
    }
 
    Self& axpy(field_type const& a, Self const& x)
    {
      this->for_each2([&](Entry& yi, Entry const& xi) { NumaDenseMatrixDetail::axpy(yi,a,xi); },x);
      return *this;
    }
    
    
    real_type two_norm () const { return std::sqrt(two_norm2()); }
    
    real_type two_norm2 () const 
    { 
      return this->accumulate([](real_type s, Entry const& e) { return s+NumaDenseMatrixDetail::frobenius_norm2(e); }, real_type(0)); 
    }
    
    real_type one_norm() const 
    { 
      return this->accumulate([](real_type s, Entry const& e) { return s+NumaDenseMatrixDetail::one_norm(e); }, real_type(0)); 
    }
    
    field_type dot(Self const& other) const { return Base::dot(other); }
    
  };
  
  template <class Entry>
  std::ostream& operator<<(std::ostream& out, NumaVector<Entry> const& vec)
  {
    using namespace std;
    copy(begin(vec),end(vec),ostream_iterator<Entry>(out," "));
    return out;
  }
}
 
#endif