threadedMatrixChunks.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) 2025-2025 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef THREADED_MATRIX_CHUNKS
#define THREADED_MATRIX_CHUNKS
 
#include <cmath>
#include <iostream>
#include <boost/timer/timer.hpp>
 
#include "dune/common/config.h"
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
 
#include "fem/firstless.hh"
#include "fem/fixdune.hh"
#include "linalg/dynamicMatrix.hh"
#include "utilities/duneInterface.hh"
#include "utilities/threading.hh"
#include "utilities/detailed_exception.hh"
 
 
namespace Kaskade
{
  // forward declarations - tell doxygen to ignore those and refer to the original definition
  template <class Entry, class Index>
  class NumaBCRSMatrix;
 
  template <class Scalar, class Index>
  class MatrixAsTriplet;
 
  template <class Target, class Source, class RowIndices, class ColIndices>
  Target submatrix(Source const& A, RowIndices const& ri, ColIndices const& ci);
 
 
  namespace ThreadedMatrixDetail
  {
    // Computes for each row the number of entries in that row. The count is stored in the
    // vector rowCount, which has to be big enough.
    template <class Entry, class Allocator, class Index>
    void getRowCount(Dune::BCRSMatrix<Entry,Allocator> const& matrix, bool isTransposed, std::vector<Index>& rowCount)
    {
      for (auto ri=matrix.begin(); ri!=matrix.end(); ++ri)
        if (isTransposed)
          for (auto ci=ri->begin(); ci!=ri->end(); ++ci)
            ++rowCount[ci.index()];
        else
          rowCount[ri.index()] = ri->size();
    }
 
    template <class Entry, class Index2, class Index>
    void getRowCount(NumaBCRSMatrix<Entry,Index2> const& matrix, bool isTransposed, std::vector<Index>& rowCount)
    {
      for (auto ri=matrix.begin(); ri!=matrix.end(); ++ri)
        if (isTransposed)
          for (auto ci=ri->begin(); ci!=ri->end(); ++ci)
            ++rowCount[ci.index()];
        else
          rowCount[ri.index()] = ri->size();
    }
 
    template <class Entry, class Index, class Index2>
    void getRowCount(MatrixAsTriplet<Entry,Index2> const& matrix, bool isTransposed, std::vector<Index>& rowCount)
    {
      auto const& indices = isTransposed? matrix.cidx: matrix.ridx;
      for (auto r: indices)
        ++rowCount[r];
    }
 
 
    // Class with partial specialization for copying from a Dune ISTL sparse matrix into ThreadedMatrix Chunks.
    // Supported modes of copying are (i) one-to-one direct copy (ii) transposed copy (iii) source is symmetric
    // (superdiagonal part is never accessed).
    template <class Entry, class Matrix, bool symmetric, bool transposed, class Index>
    struct CopyMatrixToChunk
    {
      // The simple case: not symmetric, not transposed
      static void init(Index first, Index last, Matrix const& matrix, std::vector<size_t,NumaAllocator<size_t>>& colStart,
                       std::vector<Index,NumaAllocator<Index>>& cols, std::vector<Entry,NumaAllocator<Entry>>& values,
                       std::vector<Index> const& /* nRowEntries */)
      {
        for (Index i=0; i<last-first; ++i)
        {
          Index j=0;
          for (auto ci=matrix[first+i].begin(); ci!=matrix[first+i].end(); ++ci, ++j)
          {
            values[colStart[i]+j] = *ci;
            cols[colStart[i]+j] = ci.index();
          }
        }
      }
    };
 
    template <class Entry, class Matrix, bool transposed, class Index>
    struct CopyMatrixToChunk<Entry, Matrix,true,transposed,Index>
    {
      // The hard case: symmetric with only lower triangular part to be accessed
      static void init(Index first, Index last, Matrix const& matrix, std::vector<size_t,NumaAllocator<size_t>>& colStart,
                       std::vector<Index,NumaAllocator<Index>>& cols, std::vector<Entry,NumaAllocator<Entry>>& values,
                       std::vector<Index> const& nRowEntries)
      {
        // Symmetric case is a mixture of regular and transposed. First we insert the regular entries,
        // subsequently the transposed ones.
 
        // Start with regular entries (including the diagonal)
        std::vector<Index> entriesInRow(last-first,0);
        for (Index i=0; i<last-first; ++i)
        {
          Index j=0;
          for (auto ci=matrix[first+i].begin(); ci!=matrix[first+i].end(); ++ci, ++j)
          {
            values[colStart[i]+j] = *ci;
            cols[colStart[i]+j] = ci.index();
          }
          entriesInRow[i] = matrix[first+i].size();
          assert(entriesInRow[i]<=nRowEntries[first+i]);
        }
 
        // Do the transposed entries (excluding the diagonal). Note that as we only want entries with
        // column index >= first+1, we can start (for the transposed part) at row first+1
        for (Index r=first+1; r!=matrix.N(); ++r)
          for (auto c=matrix[r].begin(); c!=matrix[r].end(); ++c)
          {
            // get transposed indices
            Index ridx = c.index();
            if (ridx>=first && ridx<last)
            {
              Index cidx = r;
 
              if (cidx > ridx) // take only superdiagonal entries
                {
                  assert(ridx-first<colStart.size() && ridx-first < entriesInRow.size());
                  Index pos = colStart[ridx-first] + entriesInRow[ridx-first];
                  assert(pos<colStart[ridx-first+1]);
                  cols[pos] = cidx;
                  values[pos] = transpose(*c);
                  ++entriesInRow[ridx-first];
                  assert(entriesInRow[ridx-first]<=nRowEntries[ridx]);
                }
            }
          }
      }
    };
 
 
    template <class Entry, class Matrix, class Index>
    struct CopyMatrixToChunk<Entry, Matrix,false,true,Index>
    {
      // The case of transposed, nonsymmetric supplied matrices.
      static void init(Index first, Index last, Matrix const& matrix, std::vector<size_t,NumaAllocator<size_t>>& colStart,
                       std::vector<Index,NumaAllocator<Index>>& cols, std::vector<Entry,NumaAllocator<Entry>>& values,
                       std::vector<Index> const& nRowEntries)
      {
        // Essentially we have to exchange row and column indices. If we go through the matrix
        // row by row, the entries of the transposed matrix to be stored in our CRS chunk come
        // with random row indices but increasing column indices. Hence we can push the entries
        // at the back of each row. We only have to remember where the current end of each row is.
        std::vector<Index> entriesInRow(last-first,0);
 
        for (auto r=matrix.begin(); r!=matrix.end(); ++r)
          for (auto c=r->begin(); c!=r->end(); ++c)
          {
            // get transposed indices
            Index ridx = c.index();
            if (ridx>=first && ridx<last)
            {
              Index cidx = r.index();
 
              assert(ridx-first<colStart.size() && ridx-first < entriesInRow.size());
              Index pos = colStart[ridx-first] + entriesInRow[ridx-first];
              assert(pos<colStart[ridx-first+1]);
              cols[pos] = cidx;
              values[pos] = transpose(*c);
              ++entriesInRow[ridx-first];
              assert(entriesInRow[ridx-first]<=nRowEntries[ridx]);
            }
          }
      }
    };
 
    //-------------------------------------------------------------------------
 
    // Class for copying matrix entries. If the entry is to be transposed during
    // copying (e.g. because only the lower triangular part of the source block matrix
    // has been stored), the matrix dimensions may change due to transposition
    // (if they are not equal). This case is covered here - we know the entry is
    // to be transposed.
    template <class To, class From, class scalarsMatch = std::false_type>
    struct MatrixEntry
    {
      static void copy(From const& from, To& to, bool isTransposed)
      {
        assert(isTransposed);
        to = transpose(from);
      }
    };
 
    // This specialization is for the case of different floating point types as scalar entries in the matrix blocks.
    // Currently this is only implemented for blocks of quadratic shape and non transposed blocks.
    template <class To, class From>
    struct MatrixEntry<To, From, typename std::is_same<typename To::field_type, typename From::field_type>::type>
    {
      static void copy(From const& from, To& to, bool isTransposed)
      {
        assert(!isTransposed);
        static_assert(To::rows == From::rows && To::cols == From::cols, "Mismatch of dimensions.");
        for(int i=0;i<To::rows;++i)
          for(int j=0;j<To::cols;++j)
            to[i][j] = static_cast<typename To::field_type>(from[i][j]);
      }
    };
 
    // The case of matching dimensions of source and target is covered here.
    // Both cases can occur, transposition or not. Partial specialization of
    // this handler class is necessary due to the runtime switch between
    // transposed or not - this dynamic switch is only legal in C++ if the
    // dimensions are the same in both cases.
    template <class Entry>
    struct MatrixEntry<Entry,Entry>
    {
      static void copy(Entry const& from, Entry& to, bool isTransposed)
      {
        if (isTransposed && from.rows > 1 && from.cols > 1 && from.rows != from.cols)
        {
          std::cerr << "Transposing rectangular Dune::FieldMatrix not implemented!" << '\n'
                    << "#rows=" << from.rows << ", #cols=" << from.cols << '\n';
          exit(1);
        }
         if (isTransposed && from.rows == from.cols)
//           to = transpose(from);
//       the following transposing code also covers the case that
//       the memory locations of from and to are the same
         {
            for (int i=0;i<from.rows;i++)
            {
              to[i][i]=from[i][i];
              for (int j=i+1;j<from.cols;j++)
              {
                typename Entry::field_type s=from[j][i];
                to[j][i]=from[i][j];
                to[i][j]=s;
              }
            }
         }
         else
           to = from;
      }
    };
 
 
    //-------------------------------------------------------------------------

    template <class Index=size_t>
    class CRSChunkPatternInfo
    {
      template <class IndexP>
      friend class CRSChunkPatternInfo;
    public:
      CRSChunkPatternInfo(Index first_, Index last_, Index cols_, bool symmetric, int node)
      : firstRow(first_), lastRow(last_), cols(cols_), numaNode(node), symm(symmetric)
      {}

      template <class IndexP>
      CRSChunkPatternInfo(CRSChunkPatternInfo<IndexP> const& pi)
      : firstRow(pi.firstRow), lastRow(pi.lastRow), cols(pi.cols)
      , numaNode(pi.numaNode), symm(pi.symm)
      {}

      Index first() const { return firstRow; }

      Index last() const { return lastRow; }

      Index columns() const { return cols; }

      int node() const { return numaNode; }

      bool symmetric() const { return symm; }
 
    protected:
      Index firstRow, lastRow; // first and one behind last row
      Index cols;              // number of columns
 
    private:
      int numaNode;
      bool symm;
    };
 
    //-------------------------------------------------------------------------

    template <class Index=size_t>
    class CRSChunkPatternCreator: public CRSChunkPatternInfo<Index>
    {
    public:
      typedef std::vector<Index> IndexArray; // column indices are indexed by user-provided Index type
      // We could use a Numa allocator here, but testing as of 2014-03-12 revealed that the performance is
      // *much* worse with Numa allocators, even though the Numa allocator is actually faster in
      // synthetic tests of a structure similar to here. This reason for this disparity is up to now unknown.
 

      CRSChunkPatternCreator(Index firstRow, Index lastRow, Index ncols, bool symmetric, int node)
      : CRSChunkPatternInfo<Index>(firstRow,lastRow,ncols,symmetric,node),
        cols(lastRow-firstRow)
      {
      }

      void reserve(Index nnzPerRow=8)
      {
        // perform allocation
        for (auto& r: cols)
          r.reserve(nnzPerRow);
      }

      void clear()
      {
        for (auto& r: cols)
        {
          r.clear();
          r.shrink_to_fit();
        }
      }

      template <class IterRow, class IterCol>
      void addElements(IterRow fromRow, IterRow const toRow,
                       IterCol const fromCol, IterCol const toCol, bool colIsSorted=false)
      {
        // sort column indices if needed, and remove duplicates (set_union below does NOT remove duplicates...)
        sortedCols.clear();
        sortedCols.insert(end(sortedCols),fromCol,toCol);
        if (!colIsSorted)
          std::sort(begin(sortedCols),end(sortedCols));
        sortedCols.erase(std::unique(begin(sortedCols),end(sortedCols)),end(sortedCols));
 
        // For each row, merge the specified column indices into the column indices that are already present.
        for ( ; fromRow != toRow; ++fromRow)                         // step through all affected rows
          if (this->first() <= *fromRow && *fromRow < this->last())  // but treat only those that actually lie in our chunk
          {
            IndexArray& c = cols[*fromRow-this->first()];
            auto top = this->symmetric()?
                              std::upper_bound(begin(sortedCols),end(sortedCols),*fromRow):
                              end(sortedCols); // omit superdiagonal elements if symmetric
            tmp.resize(c.size() + std::distance(begin(sortedCols),top));         // reserve space for the union
            tmp.erase(std::set_union(begin(c),end(c),begin(sortedCols),top,      // compute the union
                                     begin(tmp)),end(tmp));                      // and cut the container at the proper end
            std::swap(tmp,c);                                                    // write the sorted range into its target
          }
        // Note that instead of a temporary buffer tmp we could append the new column indices and use
        // std::inplace_merge. This, however, uses a temporary buffer of its own behind the scenes, one that is held
        // by the STL library implementation somewhere. It is highly probable that this memory is *not* located
        // on our NUMA node, but frequently accessed from all nodes. Hence, false sharing is bound to occur with
        // inplace_merge even if it performs proper locking.
      }

      void addAllElements(Index columns)
      {
        for (Index i=0; i<cols.size(); ++i)
        {
          cols[i].resize(this->symmetric()? i+this->first()+1: columns);
          std::iota(begin(cols[i]),end(cols[i]),static_cast<Index>(0));
        }
      }

      IndexArray const& row(Index i) const { return cols[i-this->first()]; }

      size_t nonzeroes() const
      {
        size_t nnz = 0;
        for (auto const& r: cols)
          nnz += r.size();
        return nnz;
      }

      size_t balanceForward(size_t const covered, size_t const nnz, int chunks, std::vector<IndexArray>& moveRows);

      size_t balanceBackward(size_t covered, size_t const nnz, int chunks, std::vector<IndexArray>& moveRows);
 
    private:
      std::vector<IndexArray> cols;      // For each row a std::vector of column indices
      IndexArray tmp;                    // some temporary buffer
      IndexArray sortedCols;             // some temporary buffer for sorted column indices
    };
 
    //-------------------------------------------------------------------------

    template <class Index=size_t>
    class CRSChunkPattern: public CRSChunkPatternInfo<Index>
    {
    public:


      CRSChunkPattern(CRSChunkPatternCreator<Index> const& creator);

      template <class IndexP>
      CRSChunkPattern(CRSChunkPattern<IndexP> const& p)
      : CRSChunkPatternInfo<Index>(p)
      , colStarts(p.colStarts)
      , cols(p.cols.begin(),p.cols.end(),NumaAllocator<Index>(p.node()))
      , nnzPerRow(p.nnzPerRow)
      {}

      template <class Expanded, class Condensed, class Matrix>
      CRSChunkPattern(Index first, Index last, Expanded const& eIndices, Condensed const& cIndices,
                      Matrix const& mat, int node)
      : CRSChunkPatternInfo<Index>(first,last,eIndices.size(),false,node)
      , colStarts(last-first+1,0,NumaAllocator<size_t>(node)), cols(NumaAllocator<Index>(node))
      {
        cols.reserve(last-first);                  // prevent frequent reallocation (heuristic)
 
        for (Index i=first; i<last; ++i)           // scan all our rows
        {
          colStarts[i-first] = cols.size();        // note where this row starts
 
          auto const& row = mat[eIndices[i]];      // get a handle to the column indices in that row
 
          auto cend = row.end();
          for (auto c = row.begin(); c!=cend; ++c) // step through all column entries in that row
          {
            Index ci = cIndices[c.index()];         // extract the condensed index
            if (ci<eIndices.size())                 // if that is in the condensed matrix range...
              cols.push_back(ci);                   // ...include it
          }
        }
 
        colStarts.back() = cols.size();            // sentinel
 
        // investigate sparsity
        nnzPerRow = this->last()==this->first()? 0 : nonzeroes() / (this->last()-this->first());
      }

      template <class Matrix>
      CRSChunkPattern(Matrix const& matrix, bool isSymmetric, bool isTransposed,
                      Index firstRow, Index lastRow, bool symmetric, int node)
      : CRSChunkPatternInfo<Index>(firstRow,lastRow,matrix.M(),symmetric,node)
      , colStarts(NumaAllocator<size_t>(node)), cols(NumaAllocator<Index>(node))
      {
        if (isSymmetric)         // Make sure that transposed is only flagged if
          isTransposed = false;  // it's not symmetric.
 
        // Compute start of row indices. Step through the supplied matrix and count all elements that
        // fall in one of our rows.
        std::vector<size_t,NumaAllocator<size_t>> rowCount(this->last()-this->first(),0,NumaAllocator<size_t>(node));
 
        // Compute the row range of the supplied matrix we have to scan.
        Index fromRow = this->first(), toRow = this->last();
        if (isSymmetric && !symmetric) // If the supplied matrix is symmetrically stored but we are not, we have to
          toRow = matrix.N();          // scan the later rows for the elements that appear in later columns here.
        if (isTransposed)                   // If the supplied matrix is transposed, we have to scan all the rows
          fromRow = 0, toRow = matrix.N();  // in order to cover all the columns in our range.
 
        for (Index ri=fromRow; ri<toRow; ++ri)
        {
          auto cend = matrix[ri].end();
          for (auto ci=matrix[ri].begin(); ci!=cend; ++ci)
          {
            Index r = ri, c = ci.index();
            if (isTransposed) // supplied matrix is transposed, swap row and column indices
              std::swap(r,c); // to get the row/col indices in our world
            if (this->first()<=r && r<this->last() && (!symmetric || c<=r)) // yep, this element is in our chunk
              ++rowCount[r-this->first()];                                  //
            if (isSymmetric && !symmetric && c<r && this->first()<=c && c<this->last()) // if supplied matrix is stored symmetrically but we are not
              ++rowCount[c-this->first()];                                              // copy truely subdiagonal entries to the upper triangle
          }
        }
 
        // compute the partial sums of number of elements per row.
        colStarts.resize(rowCount.size()+1,0);
        std::partial_sum(rowCount.begin(),rowCount.end(),colStarts.begin()+1);
 
        // Get space for all nonzeros
        cols.resize(colStarts.back());
 
        // Copy the column indices of elements.
        for (Index ri=fromRow; ri<toRow; ++ri)
        {
          auto cend = matrix[ri].end();
          for (auto ci=matrix[ri].begin(); ci!=cend; ++ci)
          {
            Index r = ri, c = ci.index();
            if (isTransposed) // supplied matrix is transposed, swap row and column indices
              std::swap(r,c); // to get the row/col indices in our world
            if (this->first()<=r && r<this->last() && (!symmetric || c<=r))    // yep, this element is in our chunk
            {
              cols[colStarts[r-this->first()]+rowCount[r-this->first()]-1] = c; // enter the column index
              --rowCount[r-this->first()];                                      // one less element to come
            }
            if (isSymmetric && !symmetric && c<r && this->first()<=c && c<this->last()) // if supplied matrix is stored symmetrically but we are not
            {
              cols[colStarts[c-this->first()]+rowCount[c-this->first()]-1] = r;         // copy truely subdiagonal entries to the upper triangle
              --rowCount[c-this->first()];                                              // one less element to come
            }
          }
        }
 
        // sort the column indices in each row
        for (Index ri=0; ri<rowCount.size(); ++ri) {
          assert(rowCount[ri] == 0);
          std::sort(cols.begin()+colStarts[ri],cols.begin()+colStarts[ri+1]);
        }
 
        // investigate sparsity
        nnzPerRow = this->last()==this->first()? 0 : nonzeroes() / (this->last()-this->first());
      }

      template <class Scalar, class Index2>
      CRSChunkPattern(MatrixAsTriplet<Scalar,Index2> const& matrix, bool isSymmetric, bool isTransposed,  // TODO: fuse with constructor above
                      Index firstRow, Index lastRow, bool symmetric, int node)
      : CRSChunkPatternInfo<Index>(firstRow,lastRow,matrix.M(),symmetric,node)
      , colStarts(NumaAllocator<size_t>(node)), cols(NumaAllocator<Index>(node))
      {
        if (isSymmetric)         // Make sure that transposed is only flagged if
          isTransposed = false;  // it's not symmetric.
 
        // Compute start of row indices. Step through the supplied matrix and count all elements that
        // fall in one of our rows.
        std::vector<size_t,NumaAllocator<size_t>> rowCount(this->last()-this->first(),0,NumaAllocator<size_t>(node));
 
        for (size_t i=0; i<matrix.ridx.size(); ++i)
        {
          Index r = matrix.ridx[i], c = matrix.cidx[i];
          if (isTransposed) // supplied matrix is transposed, swap row and column indices
            std::swap(r,c); // to get the row/col indices in our world
          if (this->first()<=r && r<this->last() && (!symmetric || c<=r)) // yep, this element is in our chunk
            ++rowCount[r-this->first()];                                  //
          if (isSymmetric && !symmetric && c<r && this->first()<=c && c<this->last()) // if supplied matrix is stored symmetrically but we are not
            ++rowCount[c-this->first()];                                              // copy truely subdiagonal entries to the upper triangle
        }
 
        // compute the partial sums of number of elements per row.
        colStarts.resize(rowCount.size()+1,0);
        std::partial_sum(rowCount.begin(),rowCount.end(),colStarts.begin()+1);
 
        // Get space for all nonzeros
        cols.resize(colStarts.back());
 
        // Copy the column indices of elements.
        for (size_t i=0; i<matrix.ridx.size(); ++i)
        {
          Index r = matrix.ridx[i], c = matrix.cidx[i];
          if (isTransposed) // supplied matrix is transposed, swap row and column indices
            std::swap(r,c); // to get the row/col indices in our world
          if (this->first()<=r && r<this->last() && (!symmetric || c<=r))    // yep, this element is in our chunk
          {
            cols[colStarts[r-this->first()]+rowCount[r-this->first()]-1] = c; // enter the column index
            --rowCount[r-this->first()];                                      // one less element to come
          }
          if (isSymmetric && !symmetric && c<r && this->first()<=c && c<this->last()) // if supplied matrix is stored symmetrically but we are not
          {
            cols[colStarts[c-this->first()]+rowCount[c-this->first()]-1] = r;         // copy truely subdiagonal entries to the upper triangle
            --rowCount[c-this->first()];                                              // one less element to come
          }
        }
 
        // sort the column indices in each row
        for (Index ri=0; ri<rowCount.size(); ++ri)    // TODO: do this in parallel on our node?
        {
          assert(rowCount[ri] == 0);
          auto first = cols.begin()+colStarts[ri];
          auto last = cols.begin()+colStarts[ri+1];
          std::sort(first,last);                      // sort in ascending order
          last = std::unique(first,last);             // duplicate entries may come from triplet matrices, remove
          rowCount[ri] = last-first;                  // now we have the correct number of entries.
        }
 
        // Merging duplicate column indices in a row may have reduced the number of entries in that row,
        // hence there may be a gap to the start of the next row. Compactify the column index array now.
        for (Index ri=1; ri<rowCount.size(); ++ri)
        {
          Index newStart = colStarts[ri-1] + rowCount[ri-1];
          if (newStart != colStarts[ri])                // entries need to be shifted
          {
            std::copy(cols.begin()+colStarts[ri],cols.begin()+colStarts[ri]+rowCount[ri],cols.begin()+newStart);
            colStarts[ri] = newStart;
          }
        }
        if (!rowCount.empty())                                  // set the sentinel to point just behind the last entry...
          colStarts.back() = colStarts[rowCount.size()-1]+rowCount.back();
        cols.erase(cols.begin()+colStarts.back(),cols.end());   // ...and drop the rubbish left over in the tail
 
 
        // investigate sparsity
        nnzPerRow = this->last()==this->first()? 0 : nonzeroes() / (this->last()-this->first());
      }

 
      // returns the column index of entry at position idx
      Index col(size_t idx) const { return cols[idx]; }

      size_t colStart(Index row) const { return colStarts[row]; }

      typename std::vector<Index,NumaAllocator<Index>>::const_iterator colStartIterator(Index row) const
      {
        return cols.begin()+colStart(row);
      }

      size_t position(Index r, Index c) const
      {
        auto b = cols.begin();
        r -= this->first();
        assert(0 <= r && r+1 < colStarts.size());
        auto p = std::lower_bound(b+colStarts[r],b+colStarts[r+1],c); // col indices are sorted in each row
        if (p==b+colStarts[r+1] || *p!=c)                             // return sentinel in case the entry is not here
          return std::numeric_limits<size_t>::max();
        return p-b;                                                   // return offset
      }

      bool exists(Index r, Index c) const
      {
        return position(r,c) != std::numeric_limits<size_t>::max();
      }

      size_t storage() const { return cols.size(); }

      size_t nonzeroes() const;

      Index nonzeroesPerRow() const { return nnzPerRow; }
 
    private:
      template <class IndexP>           // Provide internal access even if the index type
      friend class CRSChunkPattern;     // does not coincide.
 
      // Raw data storage. The allocators used are NUMA allocators, guaranteeing local memory access.
      std::vector<size_t,NumaAllocator<size_t>> colStarts; // memory positions are indexed by size_t (since there may be MANY)
      std::vector<Index,NumaAllocator<Index>> cols;        // column indices are indexed by user-provided Index type
      Index nnzPerRow;                                     // average number of nonzeroes per row (rounded)
    };
 
    //----------------------------------------------------------------------------------------------
    //----------------------------------------------------------------------------------------------
 
 
    template <class Arguments, class Operation>
    class NumaBCRSMatrixExpressionChunk
    {
    public:
      class iterator
      {
      };
 
      iterator begin(size_t) const
      {
      }
 
      iterator end(size_t) const
      {
      }
    };
 
    //----------------------------------------------------------------------------------------------

    template <class Entry, class Index=size_t>
    class CRSChunk
    {
      using Self = CRSChunk<Entry,Index>;
 
      template <class EntryA, class IndexA>
      friend class CRSChunk;
 
    public:
      using Scalar = typename EntryTraits<Entry>::field_type;


      CRSChunk(int node)
      : values(NumaAllocator<Entry>(node))
      {
      }

      template <class EntryA, class IndexA,
                typename = std::enable_if_t<std::is_convertible_v<EntryA,Entry>>>
      CRSChunk(std::shared_ptr<CRSChunkPattern<Index>> const& pattern,
              CRSChunk<EntryA,IndexA> const& c)
      : pat(pattern)
      , scatterMutex(4)
      , values(c.values.begin(),c.values.end(),NumaAllocator<Entry>(c.pat->node()))
      {
      }

      CRSChunk(std::shared_ptr<CRSChunkPattern<Index>> const& pattern_, Entry const& init)
      : pat(pattern_)
      , scatterMutex(4)
      , values(pat->nonzeroes(),init,NumaAllocator<Entry>(pat->node()))
      {
      }

      template <class Matrix>
      CRSChunk(std::shared_ptr<CRSChunkPattern<Index>> const& pattern_, Matrix const& matrix,
              bool isSymmetric, bool isTransposed)
      : CRSChunk(pattern_,Entry(0))
      {
        typedef typename Matrix::block_type SuppliedEntry;
 
        Index first = pat->first(), last = pat->last();
 
        if (isSymmetric)         // Make sure that transposed is only flagged if
          isTransposed = false;  // it's not symmetric.
 
        // Compute the row range of the supplied matrix we have to scan.
        Index fromRow = first, toRow = last;
        if (isSymmetric && !pat->symmetric()) // If the supplied matrix is symmetrically stored but we are not, we have to
          toRow = matrix.N();                 // scan the later rows for the elements that appear in later columns here.
        if (isTransposed)                     // If the supplied matrix is transposed, we have to scan all the rows
          fromRow = 0, toRow = matrix.N();    // in order to cover all the columns in our range.
 
        // Copy the column indices of elements.
        for (Index ri=fromRow; ri<toRow; ++ri)
        {
          auto cend = matrix[ri].end();
          for (auto ci=matrix[ri].begin(); ci!=cend; ++ci)
          {
            Index r = ri, c = ci.index();
            if (isTransposed) // supplied matrix is transposed, swap row and column indices
              std::swap(r,c); // to get the row/col indices in our world
            if (first<=r && r<last && (!pat->symmetric() || c<=r))                                   // yep, this element is in our chunk
            {
              assert(pat->position(r,c) < values.size());
              MatrixEntry<Entry,SuppliedEntry>::copy(*ci,values[pat->position(r,c)],isTransposed);   // copy, transpose if needed
            }
            if (isSymmetric && !pat->symmetric() && c<r && first<=c && c<last)                       // if supplied matrix is stored symmetrically but we are not..
              MatrixEntry<Entry,SuppliedEntry>::copy(*ci,values[pat->position(c,r)],true);           // copy, transpose
          }
        }
      }

      template <class Scalar, class Index2>
      CRSChunk(std::shared_ptr<CRSChunkPattern<Index>> const& pattern_, MatrixAsTriplet<Scalar,Index2> const& matrix, bool isSymmetric, bool isTransposed)
      : CRSChunk(pattern_,Entry(0))
      {
        Index first = pat->first(), last = pat->last();
 
        if (isSymmetric)         // Make sure that transposed is only flagged if
          isTransposed = false;  // it's not symmetric.
 
 
        // Copy the column indices of elements.
        for (size_t i=0; i<matrix.ridx.size(); ++i)
        {
          Index r = matrix.ridx[i], c = matrix.cidx[i];
          if (isTransposed) // supplied matrix is transposed, swap row and column indices
            std::swap(r,c); // to get the row/col indices in our world
          if (first<=r && r<last && (!pat->symmetric() || c<=r))                                   // yep, this element is in our chunk
          {
            assert(pat->position(r,c) < values.size());
            values[pat->position(r,c)] = matrix.data[i];                                           // copy (no transpose as triplet entries are scalar)
          }
          if (isSymmetric && !pat->symmetric() && c<r && first<=c && c<last)                       // if supplied matrix is stored symmetrically but we are not..
            values[pat->position(c,r)] = matrix.data[i];                                           // copy (no transpose as triplet entries are scalar)
        }
      }

      template <class Matrix>
      CRSChunk(Matrix const& matrix, bool isSymmetric, bool isTransposed,
              Index firstRow, Index lastRow, bool symmetric, int node)
      : CRSChunk(std::make_shared<CRSChunkPattern<Index>>(matrix,isSymmetric,isTransposed,
                                                          firstRow,lastRow,symmetric,node),
                matrix,isSymmetric,isTransposed)
      {}

      template <class Expanded, class Condensed,  class Matrix>
      CRSChunk(std::shared_ptr<CRSChunkPattern<Index>> const& pattern_,
              Expanded const& eIndices, Condensed const& cIndices, Matrix const& matrix)
      : CRSChunk(pattern_,Entry(0))
      {
        Index first = pat->first(), last = pat->last();
        for (Index r=first; r<last; ++r)
        {
          auto row = matrix[eIndices[r]];
          auto p = begin(values) + pat->colStart(r-first);
          for (auto c=row.begin(); c!=row.end(); ++c)
            if (cIndices[c.index()]<eIndices.size())
            {
              *p = *c;
              assert(!std::isnan(*p));
              ++p;
            }
 
          assert(p==begin(values)+pat->colStart(r-first+1));
        }
      }


      Self& operator=(Self const& c) = default;

      template <class Value>
      Self&  operator=(Value const& a)
      {
        std::fill(begin(values),end(values),a);
        return *this;
      }
 

      template <class Arguments, class Operation>
      Self& operator=(ThreadedMatrixDetail::NumaBCRSMatrixExpressionChunk<Arguments,Operation> const& e)
      {
        Index first = pat->first();
        for (Index r=0; r<pat->last()-first; ++r)
        {
          auto eend = e.end(r);
          auto ci = pat->colStartIterator(r);
          auto cend = pat->colStartIterator(r+1);
          auto vi = begin(values)+pat->colStart(r);
          for (auto ei = e.begin(r); ei != eend; ++ei) // step through the expression row
          {
            auto pos = std::find(ci,cend,ei.index());  // find our entry
            assert(pos != cend);                       // our pattern must be a superset of the expression pattern
            vi += pos-ci;                              // advance to found entry
            ci = pos;
            *vi = *ei;                                 // assign
          }
        }
      }

      template <class Domain, class Range>
      Scalar apply(Scalar a, Domain const& x, Range& y, bool initialize) const
      {
        Scalar dp = 0;
        auto const& p = *pat;
        Index first = p.first();
        Index last = p.last();
 
        // Perform matrix vector multiplication. Step through all rows in our chunk...
        // On Numa nodes with several cores, one may think that using multiple threads
        // on this node could improve the performance. Alas, this appears not to be the
        // case as of 2014-05 (both aged AMD Opteron 32cores/8nodes and a newer 32/4 Intel
        // Xeon). It seems that a single thread saturates the memory channel(s).
        // Hence we stay with the sequential version here.
        for (Index row=first; row<last; ++row)
        {
          // and compute the scalar product of (sparse) row and vector
          typename Range::block_type z(0);
 
          size_t jEnd = p.colStart(row+1-first);
          for (size_t j=p.colStart(row-first); j<jEnd; ++j)
            if (Entry::rows==1 && Entry::cols==1)
              z.axpy(values[j][0][0],x[p.col(j)]);              // this is necessary as otherwise scalar matrix entries and
            else                                                // vecotr-valued rhs entries do not work together
              values[j].umv(x[p.col(j)],z);
 
          // scale result as requested
          z *= a;
          if (initialize)   y[row] = z;
          else              y[row] += z;
 
          // accumulate duality product
          if (Entry::rows==Entry::cols && row<x.size())
            dp += y[row]*x[row];
        }
 
        return dp;
      }

      template <class Domain, class Range>
      void applyTransposed(Scalar a, Domain const& x, Range& y) const
      {
        auto const& p = *pat;
        Index first = p.first(); // first row in our chunk
        Index last = p.last();   // last row in our chunk
 
        assert(x.size()>=last);
 
        // We compute y_i = a sum_j A_ji^T x_j, with j in the outer loop
        for (Index j=first; j<last; ++j)                        // looping over the rows of A, i.e. the columns of A^T
        {
          size_t jEnd = p.colStart(j+1-first);
          for (size_t k=p.colStart(j-first); k<jEnd; ++k)       // looping over the columns in the row
          {
            auto i = p.col(k);                                  // the column index in A, i.e. the row index in A^T
            if (Entry::rows==1 && Entry::cols==1)
              y[i] += a*values[k][0][0] * x[j];                 // this is necessary as otherwise scalar matrix entries and
            else                                                // vecotr-valued rhs entries do not work together
              values[k].usmtv(a,x[j],y[i]);                     // y_i += a * (A^T)_ij x_j = a * (A_ji)^T x_j
          }
        }
      }

      template <class Domain, class Range>
      void gatherMirrored(Scalar a, Domain const& x, Range& y, CRSChunk<Entry,Index> const& block, bool subdiagonal, bool initToZero) const
      {
        auto const& p = *block.pat;
        Index firstRow = p.first();
        Index lastRow = p.last();
        Index firstCol = pat->first();
        Index lastCol = pat->last();
 
        assert(lastRow<=y.size());
 
        if (initToZero)
          for (Index row=firstRow; row<lastRow; ++row)
            y[row] = 0;
 
        for (Index row=firstRow; row<lastRow; ++row)
        {
          auto const& xrow = x[row];
 
          auto cend = p.colStartIterator(row+1-firstRow);
          auto ci = firstCol==0?                                                       // We consider all elements with column index in our
                    p.colStartIterator(row-firstRow):                                  // row chunk. Compute the starting point by binary search
                    std::lower_bound(p.colStartIterator(row-firstRow),cend,firstCol);  // (or used the shortcut if we know we start from the beginning).
          auto vi = block.values.begin() + (ci-p.colStartIterator(0));
          if (subdiagonal)
            for ( ; ci!=cend && *ci<lastCol && *ci<row; ++ci, ++vi)                    // Then step through the row until we hit the upper bound,
              vi->usmtv(a,xrow,y[*ci]);                                                // which is either the upper row range end or the diagonal,
          else                                                                         // depending on the matrix symmetry (encoded in subdiagonal).
            for ( ; ci!=cend && *ci<lastCol; ++ci, ++vi)
              vi->usmtv(a,xrow,y[*ci]);
        }
      }

      CRSChunkPattern<Index> const& pattern() const { return *pat; }

      auto valStartIterator(Index row)
      {
        return values.begin()+pat->colStart(row);
      }

      template <class LMIterator>
      void scatter(LMIterator first, LMIterator last)
      {
        Index firstRow = pat->first();
        Index lastRow = pat->last();
 
        if (first==last || firstRow==lastRow)    // empty data or empty chunk
          return;                                // -> nothing to do
 
        Index nnzPerRow = pat->nonzeroesPerRow();
 
        // scatter into uniform ranges
        // TODO: does it make sense to parallelize this within each NUMA domain? Right now it
        // is not parallelized, only prepared.
        Index n = scatterMutex.size();
        for (Index i=0; i<n; ++i)
        {
          boost::lock_guard<boost::mutex> lock(scatterMutex[i].get());
 
          Index rowRangeStart = uniformWeightRangeStart(i,  n,lastRow-firstRow)+firstRow;
          Index rowRangeEnd   = uniformWeightRangeStart(i+1,n,lastRow-firstRow)+firstRow;
 
          if (nnzPerRow == pat->columns())                     // completely dense
            scatterRowRange<0>(first,last,firstRow,rowRangeStart,rowRangeEnd);
          else if (nnzPerRow > 10*(first->cidx().size()))      // many more entries than we scatter right now
            scatterRowRange<1>(first,last,firstRow,rowRangeStart,rowRangeEnd);
          else                                                 // rather few entries per row
            scatterRowRange<2>(first,last,firstRow,rowRangeStart,rowRangeEnd);
        }
      }

      template <class F, class EntryB, class IndexB>
      void entrywiseOp(F const& f, NumaBCRSMatrix<EntryB,IndexB> const& B)
      {
        auto const& p = *pat;
        Index first = p.first();  // first row index
        Index last = p.last();    // end row index
 
        for (Index row=first; row<last; ++row)
        {
          auto const& Br = B[row];                              // get the other matrix row
          size_t jEnd = p.colStart(row+1-first);                // step through all entries in this row
          auto bri = Br.begin();
          auto brend = Br.end();
          for (size_t j=p.colStart(row-first); j<jEnd; ++j)
          {
            while (bri!= brend && bri.index()<p.col(j))         // advance other row pointer until we find our entry (but do not overshoot)
              ++bri;
            if (bri!=brend && bri.index()==p.col(j))            // if the column indices coincide, apply the functor
              values[j] = f(values[j],*bri);
          }
        }
      }
 
    private:
      std::shared_ptr<CRSChunkPattern<Index>> pat;
      std::vector<Mutex> scatterMutex;
 
      // Raw data storage. The allocators used are NUMA allocators to guarantee local memory access.
      std::vector<Entry,NumaAllocator<Entry>> values;     // matrix entries
 
      // scatters given local matrices into the provided row range. The expected sparsity is
      // encoded as static template parameter in order to allow an efficient if-statement in
      // the innermost loop
      template <int sparsity, class LMIterator>
      void scatterRowRange(LMIterator first, LMIterator last, Index startRow, Index firstRow, Index lastRow)
      {
        // simple-minded scatter implementation
        for ( ; first!=last; ++first)    // step through all local matrices
        {
          for (auto rgl: first->ridx())  // step through all global rows affected by this matrix
          {
            Index rg = rgl.first;
            if (rg >= lastRow)           // global row indices in rgl.first are sorted ascendingly,
              break;                     // if we pass our last row, we can stop with this local matrix
            if (firstRow<=rg)
              scatterRow<sparsity>(*first,rg-startRow,rgl.second);
          }
        }
      }
 
      template <int sparsity, class LM, class LIndex>
      void scatterRow(LM const& a, Index rg, LIndex rl)
      {
        auto cbegin = pat->colStartIterator(rg);
        auto cend = pat->colStartIterator(rg+1);
        auto vbegin = valStartIterator(rg);
 
        if (LM::lumped) // only diagonal, i.e. column index is row index both in global and local matrix
        {
          auto pos = std::find(cbegin,cend,rg+pat->first());
          assert(pos != cend);
          vbegin[pos-cbegin] += a(rl,rl);
        }
        else
        {
          for (auto cgl: a.cidx())
          {
            auto pos = sparsity==0? std::lower_bound(cbegin,cend,cgl.first) // find entry in row with global column number cgl.first
                    : sparsity==1? std::lower_bound(cbegin,cend,cgl.first) // TODO: specialized find for dense matrices.
                    :              std::find(cbegin,cend,cgl.first);
 
            if (pos==cend)
              return;                                    // reached end of row - skip the rest (there may be more, e.g. if we are symmetric)
 
            vbegin += pos-cbegin;                        // move on to the found entry
            cbegin =  pos;
 
            *vbegin += a(rl,cgl.second);                 // add up appropriate entry from local matrix
          }
        }
      }
    };
 
    // ---------------------------------------------------------------------------------------------
  }
}
 
 
#endif