assemblerMatrixExtraction.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) 2024-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef ASSEMBLERMATRIXEXTRACTION_HH
#define ASSEMBLERMATRIXEXTRACTION_HH
 
#include "dune/common/config.h"
 
#include "linalg/triplet.hh"
 
// forward declarations
namespace Dune
{
  template <class Entry, class Allocator>
  class BCRSMatrix;
}
 
namespace Kaskade
{
  template <class Entry, class Index>
  class NumaBCRSMatrix;
}
 
namespace Kaskade::AssemblyDetail
{
  template <class IdxOutIter, class DataOutIter>
  struct BlockToTriplet
  {
    BlockToTriplet(size_t firstRowBlock_, size_t firstColumnBlock_, std::vector<size_t> const& rowOff_, std::vector<size_t> const& colOff_,
        IdxOutIter& ri_, IdxOutIter& ci_, DataOutIter& xi_, bool onlyLowerTriangle_):
          firstRowBlock(firstRowBlock_), firstColumnBlock(firstColumnBlock_), rowOff(rowOff_), colOff(colOff_),
          ri(ri_), ci(ci_), xi(xi_), onlyLowerTriangle(onlyLowerTriangle_)
    {}
 
    template <class MatrixBlock>
    void operator()(MatrixBlock const& mb) const
    {
      // Check if block is in requested range
      if (inRange(mb.rowId,mb.colId))
        Matrix_to_Triplet<typename MatrixBlock::Matrix>::call(mb.globalMatrix(),ri,ci,xi,
                                                        rowOff[MatrixBlock::rowId-firstRowBlock],colOff[MatrixBlock::colId-firstColumnBlock],
                                                        mb.rowId==mb.colId && onlyLowerTriangle,
                                                        mb.symmetric);
      // For mirror blocks, the transposed block needs to be written
      // if the transposed block is in the requested
      // range. Transposition is implicitly achieved by swapping
      // column and row index output iterators.
      if (MatrixBlock::mirror && inRange(mb.colId,mb.rowId) && onlyLowerTriangle==false)
        Matrix_to_Triplet<typename MatrixBlock::Matrix>::call(mb.globalMatrix(),ci,ri,xi,
            colOff[MatrixBlock::rowId-firstColumnBlock],rowOff[MatrixBlock::colId-firstRowBlock],
            false,
            mb.symmetric);
    }
 
    bool inRange(size_t r, size_t c) const
    {
      return r>=firstRowBlock && r<firstRowBlock+rowOff.size() && c>=firstColumnBlock && c<firstColumnBlock+colOff.size();
    }
 
 
    size_t firstRowBlock, firstColumnBlock;
    std::vector<size_t> const& rowOff;
    std::vector<size_t> const& colOff;
    IdxOutIter& ri;
    IdxOutIter& ci;
    DataOutIter& xi;
    bool onlyLowerTriangle;
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  template <class Matrix> struct Fill;
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  template <class Scalar, class SparseIndexInt>
  struct Fill<MatrixAsTriplet<Scalar,SparseIndexInt>>
  {
    typedef MatrixAsTriplet<Scalar,SparseIndexInt> Matrix;
 
    template <class MatrixBlockArray>
    static Matrix apply(MatrixBlockArray const& block,
                        size_t firstRowBlock, size_t firstColumnBlock,
                        std::vector<size_t> const& rowOff, std::vector<size_t> const& colOff,
                        bool extractOnlyLowerTriangle, size_t nnz, size_t, size_t)
    {
      Matrix result(nnz);
 
      typename Matrix::index_iterator rowIterator = result.ridx.begin(),
                                      columnIterator = result.cidx.begin();
      typename Matrix::iterator dataIterator = result.data.begin();
      for_each(block,BlockToTriplet<typename Matrix::index_iterator,
                                    typename Matrix::iterator>(firstRowBlock,firstColumnBlock,
                                                               rowOff,colOff,
                                                               rowIterator,columnIterator,dataIterator, extractOnlyLowerTriangle));
 
      return result;
    }
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  template <class Scalar, class Allocator>
  struct Fill<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>,Allocator> >
  {
    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>,Allocator> Matrix;
 
    template <class MatrixBlockArray>
    static Matrix apply(MatrixBlockArray const& block, size_t firstRowBlock, size_t firstColumnBlock,
                        std::vector<size_t> const& rowOff, std::vector<size_t> const& colOff,
                        bool onlyLowerTriangle, size_t nnz, size_t nrows, size_t ncols)
    {
      // read as triplet
      typedef MatrixAsTriplet<Scalar,size_t> Triplet;
      Triplet triplet( Fill<Triplet>::apply(block, firstRowBlock, firstColumnBlock, rowOff, colOff, onlyLowerTriangle, nnz, nrows, ncols) );
 
      // get column indices for each row
     auto bcrsids = getBCRSIndicesFromTriplet(nrows, triplet.ridx, triplet.cidx);
 
      Matrix result(nrows,ncols,nnz,Matrix::row_wise);
      // create sparsity pattern
      for (typename Matrix::CreateIterator row=result.createbegin(); row!=result.createend(); ++row)
        for(size_t col = 0; col < bcrsids[row.index()].size(); ++col)
          row.insert(bcrsids[row.index()][col]);
 
      // fill matrix
      for (size_t i=0; i<nnz; ++i)
        result[triplet.ridx[i]][triplet.cidx[i]] = triplet.data[i];
      return result;
    }
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  template <class Scalar, class Index>
  struct Fill<NumaBCRSMatrix<Dune::FieldMatrix<Scalar,1,1>,Index>>
  {
    typedef NumaBCRSMatrix<Dune::FieldMatrix<Scalar,1,1>,Index> Matrix;
 
    template <class MatrixBlockArray>
    static Matrix apply(MatrixBlockArray const& block, size_t firstRowBlock, size_t firstColumnBlock,
                        std::vector<size_t> const& rowOff, std::vector<size_t> const& colOff,
                        bool onlyLowerTriangle, size_t nnz, size_t nrows, size_t ncols)
    {
      // read as triplet
      typedef MatrixAsTriplet<Scalar,size_t> Triplet;
      Triplet triplet( Fill<Triplet>::apply(block, firstRowBlock, firstColumnBlock, rowOff, colOff, onlyLowerTriangle, nnz, nrows, ncols) );
 
      // convert to Numa matrix
      return Matrix(triplet,onlyLowerTriangle,false,onlyLowerTriangle);
    }
  };
 
 
} // end of namespace AssemblyDetail
 
#endif