matrixOps.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 MATRIXOPS_HH
#define MATRIXOPS_HH
 
namespace Kaskade 
{
  // forward declarations such that we need not include both dense and sparse matrix headers
  template<class K> class DynamicMatrix;
  template <class Entry, class Index> class NumaBCRSMatrix;
  template <class Index> class NumaCRSPatternCreator;
  
  namespace MatrixOpsDetail
  {
    template <class Target>
    class Submatrix;
    
    template <class Entry, class Allocator>
    class Submatrix<Dune::BCRSMatrix<Entry,Allocator>>
    {
    public:
      template <class Source, class RowIndices, class ColIndices>
      static Dune::BCRSMatrix<Entry,Allocator> apply(Source const& A,
                                                     RowIndices const& ri, ColIndices const& ci)
      {
        using Target = Dune::BCRSMatrix<Entry,Allocator>;
        using Index = typename Source::size_type;
        
        // Extract nonzero entries in the requested row and column ranges and store them
        // in sorted triplet format.
        std::vector<std::tuple<Index,Index>> a;
        std::vector<Entry> av;
        for (Index r=0; r<ri.size(); ++r)
        {
          auto row = A[ri[r]];
          for (Index c=0; c<ci.size(); ++c)
          {
            auto coli = row.find(ci[c]);
            if (coli != row.end())
            {
              a.push_back(std::make_tuple(r,c));
              av.push_back(*coli);
            }
          }
        }
        
        // Define the sparsity pattern.
        Target S(ri.size(),ci.size(),a.size(),Target::row_wise);
        using CreateIter = typename Target::CreateIterator;
        auto ai = begin(a);
        for(CreateIter row=S.createbegin(); row!=S.createend(); ++row)
          for ( ; ai!=end(a) && std::get<0>(*ai)==row.index(); ++ai)
            row.insert(std::get<1>(*ai));
        
        // Enter the matrix entries
        auto avi = begin(av);
        for (auto r=S.begin(); r!=S.end(); ++r)
          for (auto c=r->begin(); c!=r->end(); ++c, ++avi)
            *c = *avi;
        
        return S;
      }
    };
    
    
    template <class Entry, class Index>
    class Submatrix<NumaBCRSMatrix<Entry,Index>>
    {
    public:
      template <class Source, class RowIndices, class ColIndices>
      static NumaBCRSMatrix<Entry,Index> apply(Source const& A, 
                                               RowIndices const& ri, ColIndices const& ci)
      {
        NumaCRSPatternCreator<Index> creator(ri.size(), ci.size());
        
        // define sparsity pattern
        // WARNING: This has O(n*m) run time if the target matrix is n x m.
        // TODO: Devise a more sophisticated algorithm with O(nnz+n) run time.
        for(Index r=0; r<ri.size(); ++r)                                // iterate over rows in ri
        {
          auto row = A[ri[r]];
          for (Index c=0; c<ci.size(); ++c)                             // iterate over columns in ci
            if (auto coli = row.find(ci[c]); coli != row.end())
              creator.addElement(r,c);
        }
        
        NumaBCRSMatrix<Entry,Index> S(creator);
        
        // Enter the matrix entries
        for (auto r=S.begin(); r!=S.end(); ++r)
          for (auto c=r->begin(); c!=r->end(); ++c)
            *c = A[ri[r.index()]][ci[c.index()]];
        
        return S;
      }
    };
    
    template <class Entry>
    class Submatrix<DynamicMatrix<Entry>>
    {
    public:
      template <class RowIdx, class ColIdx>
      static DynamicMatrix<Entry> apply(DynamicMatrix<Entry> const& A, 
                                        RowIdx const& rows, ColIdx const& cols)
      {
        assert(*std::min_element(begin(rows),end(rows))>=0);
        assert(*std::min_element(begin(cols),end(cols))>=0);
        assert(*std::max_element(begin(rows),end(rows))<A.N());
        assert(*std::max_element(begin(cols),end(cols))<A.M());
        
        DynamicMatrix<Entry> S(size(rows),size(cols));
        auto ci = begin(cols);
        for (int j=0; ci!=end(cols); ++j, ++ci)
        {
          auto ri = begin(rows);
          for (int i=0; ri!=end(rows); ++i, ++ri)
            S[i][j] = A[*ri][*ci];
        }
        return S;
      }
      
      template <class Source, class RowIndices, class ColIndices>
      static DynamicMatrix<Entry> apply(Source const& A,
                                        RowIndices const& ri, ColIndices const& ci)
      {
        DynamicMatrix<Entry> S(ri.size(),ci.size());
        for (int r=0; r<ri.size(); ++r)
        {
          auto row = A[ri[r]];
          
          for (int c=0; c<ci.size(); ++c)
            if (auto coli = row.find(ci[c]); coli != row.end())
              S[r][c] = *coli;
        }
        
        return S;
      }
    };
  }
  
  template <class Target, class Source, class RowIndices, class ColIndices>
  Target submatrix(Source const& A, RowIndices const& ri, ColIndices const& ci)
  {
    return MatrixOpsDetail::Submatrix<Target>::apply(A,ri,ci);
  }
  
  
  
}
 
#endif