localMatrices.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) 2016-2022 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef LOCALMATRICES_HH
#define LOCALMATRICES_HH
 
#include <ostream>
#include <vector>
 
#include "linalg/threadedMatrix.hh"
 
namespace Kaskade
{
  
  template <class Entry, bool diagonal, class SortedRowIdx, class SortedColIdx>
  class LocalMatrix 
  {
  public:
    static bool const lumped = diagonal; // backwards compatibility, remove once NumaBCRSMatrix::scatter has been updated
    
    LocalMatrix(SortedRowIdx const& ridx, SortedColIdx const& cidx, Entry* data_)
    : ridx_(ridx), cidx_(cidx), data(data_)
    {
    }
    
    SortedRowIdx ridx() const { return ridx_; }
    
    SortedColIdx cidx() const { return cidx_; }
    
    size_t size() const
    {
      return diagonal? ridx().size(): ridx().size()*cidx().size();
    }
    
    void relocate(Entry* newData)
    {
      data = newData;
    }
 
    Entry& operator()(int row, int col)
    {
      // For diagonal matrix blocks, only the diagonal of local matrices are computed
      // and accessed. Store these not as a diagonal of a full matrix, but in a 
      // contiguous vector.
      if (diagonal)
      {
        assert(row==col);
        return data[row];
      }
      else
        // LAPACK dense storage scheme.
        return data[row*cidx().size()+col];
    }
    
    Entry const& operator()(int row, int col) const
    {
      // For diagonal matrix blocks, only the diagonal of local matrices are computed
      // and accessed. Store these not as a diagonal of a full matrix, but in a 
      // contiguous vector.
      if (diagonal)
      {
        assert(row==col);
        return data[row];
      }
      else
        return data[row*cidx().size()+col];
    }
 
    LocalMatrix& operator=(DynamicMatrix<Entry> const& Aloc)
    {
      assert(Aloc.N()==ridx().size() && Aloc.M()==cidx().size());
      if (diagonal)
        for (int r=0; r<Aloc.N(); ++r)
          data[r] = Aloc[r][r];
      else
        for (int c=0; c<Aloc.M(); ++c)      // TODO: consider implementing this as memcpy
          for (int r=0; r<Aloc.N(); ++r)    //       for performance reasons
            (*this)(r,c) = Aloc[r][c];
 
      return *this;
    }
    
 
  private:
    SortedRowIdx ridx_;
    SortedColIdx cidx_;
    
    // a pointer to the raw storage for the matrix entries
    Entry* data;
  };
 
  template <class Entry, bool diagonal, class SortedRowIdx, class SortedColIdx>
  std::ostream& operator <<(std::ostream& out, LocalMatrix<Entry,diagonal,SortedRowIdx,SortedColIdx> const& A)
  {
    int const n = A.ridx().size();
    int const m = A.cidx().size();
 
    for (int i=0; i<n; ++i)
    {
      for (int j=0; j<m; ++j)
        out << (diagonal && i!=j? Entry(): A(i,j)) << " ";
      out << "\n";
    }
    return out;
  }
  
  // ----------------------------------------------------------------------------------------------
  
  template <class Entry, bool diagonal, class SortedRowIdx, class SortedColIdx, 
            class IT=void, class IndexType=std::size_t>
  class LocalMatrices
  {
  public:
    using value_type = LocalMatrix<Entry,diagonal,SortedRowIdx,SortedColIdx>;
    using InfoType = IT;
    
    LocalMatrices(NumaBCRSMatrix<Entry,IndexType>& globalMatrix_, size_t maxStorage_=256*1024)
    : globalMatrix(&globalMatrix_)
    , maxStorage(maxStorage_)
    {
      // Allocate memory as desired. The vector will only grow if there is a local matrix with more entries than defined here.
      localData.reserve(maxStorage/sizeof(Entry));
    }
    
    LocalMatrices()
    : globalMatrix(nullptr), maxStorage(0)
    {}
    
    ~LocalMatrices()
    {
      scatter();
    }
    
    // Storage for local stiffness matrices.
    std::vector<value_type> localMatrices;
        
    void push_back(SortedRowIdx const& ridx, SortedColIdx const& cidx)
    {
      // for diagonal matrices, only the diagonal is stored in a contiguous vector
      size_t size = diagonal? ridx.size(): ridx.size()*cidx.size();
      assert(ridx.size()==ridx.size() || !diagonal);
      
      // Appending the new matrix exceeds the storage capacity. Scatter the existing matrices to their target and 
      // remove them.
      if (localData.capacity() < localData.size()+size)
        scatter();
      
      // Append the new local matrix. Note that either the capacity of localData is sufficient, and 
      // the localMatrices data pointers remain valid, or the local matrices have been scattered,
      // such that localMatrices is empty and no pointer exists that could be invalidated.
      Entry* pos = &*localData.insert(localData.end(),size,Entry(0));
      localMatrices.emplace_back(ridx,cidx,pos);
    }
 
    value_type      & operator[](int n)       { return localMatrices[n]; }
    value_type const& operator[](int n) const { return localMatrices[n]; }
    
    value_type      & back()       { return localMatrices.back(); }
    value_type const& back() const { return localMatrices.back(); }
    
    size_t size() const { return localMatrices.size(); }
    
    void scatter()
    {
      if (globalMatrix)
        globalMatrix->scatter(begin(localMatrices),end(localMatrices));
      clear();
    }
    
    void clear() 
    {
      localData.clear();
      localMatrices.clear();
    }
    
    size_t storageSize() const
    {
      return localData.size() * sizeof(Entry);
    }
    
    size_t storageSizeLimit() const 
    {
      return maxStorage;
    }
    
  private:
    // This allows to hold a single array into which the local matrices are scattered, and
    // thus improves memory locality (for better cache hit rates and less false sharing).
    // Contains the local matrices one after the other, for diagonal matrices, only storage
    // for the diagonal is provided.
    std::vector<Entry>               localData;      // entries of elemental matrices
    NumaBCRSMatrix<Entry,IndexType>* globalMatrix;
    size_t                           maxStorage;     // maximal desired storage space in bytes
  };
  
}
 
#endif