"Fossies" - the Fresh Open Source Software Archive

Member "dune-istl-2.7.1/dune/istl/colcompmatrix.hh" (26 Nov 2020, 18020 Bytes) of package /linux/misc/dune/dune-istl-2.7.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "colcompmatrix.hh" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2.6.0_vs_2.7.0.

    1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
    2 // vi: set et ts=4 sw=2 sts=2:
    3 #ifndef DUNE_ISTL_COLCOMPMATRIX_HH
    4 #define DUNE_ISTL_COLCOMPMATRIX_HH
    5 #include "bcrsmatrix.hh"
    6 #include "bvector.hh"
    7 #include <dune/common/fmatrix.hh>
    8 #include <dune/common/fvector.hh>
    9 #include <dune/common/typetraits.hh>
   10 #include <dune/common/unused.hh>
   11 #include <dune/common/scalarmatrixview.hh>
   12 #include <limits>
   13 
   14 namespace Dune
   15 {
   16   /**
   17    * @brief Provides access to an iterator over all matrix rows.
   18    *
   19    * @tparam M The type of the matrix.
   20    */
   21   template<class M>
   22   class MatrixRowSet
   23   {
   24   public:
   25     //! @brief The type of the matrix.
   26     typedef M Matrix;
   27     //! @brief The matrix row iterator type.
   28     typedef typename Matrix::ConstRowIterator const_iterator;
   29 
   30     /**
   31      * @brief Construct an row set over all matrix rows.
   32      * @param m The matrix for which we manage the rows.
   33      */
   34     MatrixRowSet(const Matrix& m)
   35       : m_(m)
   36     {}
   37 
   38     //! @brief Get the row iterator at the first row.
   39     const_iterator begin() const
   40     {
   41       return m_.begin();
   42     }
   43     //! @brief Get the row iterator at the end of all rows.
   44     const_iterator end() const
   45     {
   46       return m_.end();
   47     }
   48   private:
   49     const Matrix& m_;
   50   };
   51 
   52   /**
   53    * @brief Provides access to an iterator over an arbitrary subset
   54    * of matrix rows.
   55    *
   56    * @tparam M The type of the matrix.
   57    * @tparam S the type of the set of valid row indices.
   58    */
   59   template<class M, class S>
   60   class MatrixRowSubset
   61   {
   62   public:
   63     /** @brief the type of the matrix class. */
   64     typedef M Matrix;
   65     /** @brief the type of the set of valid row indices. */
   66     typedef S RowIndexSet;
   67 
   68     /**
   69      * @brief Construct an row set over all matrix rows.
   70      * @param m The matrix for which we manage the rows.
   71      * @param s The set of row indices we manage.
   72      */
   73     MatrixRowSubset(const Matrix& m, const RowIndexSet& s)
   74       : m_(m), s_(s)
   75     {}
   76 
   77     const Matrix& matrix() const
   78     {
   79       return m_;
   80     }
   81 
   82     const RowIndexSet& rowIndexSet() const
   83     {
   84       return s_;
   85     }
   86 
   87     //! @brief The matrix row iterator type.
   88     class const_iterator
   89       : public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
   90     {
   91     public:
   92       const_iterator(typename Matrix::const_iterator firstRow,
   93                      typename RowIndexSet::const_iterator pos)
   94         : firstRow_(firstRow), pos_(pos)
   95       {}
   96 
   97 
   98       const typename Matrix::row_type& dereference() const
   99       {
  100         return *(firstRow_+ *pos_);
  101       }
  102       bool equals(const const_iterator& o) const
  103       {
  104         return pos_==o.pos_;
  105       }
  106       void increment()
  107       {
  108         ++pos_;
  109       }
  110       typename RowIndexSet::value_type index() const
  111       {
  112         return *pos_;
  113       }
  114 
  115     private:
  116       //! @brief Iterator pointing to the first row of the matrix.
  117       typename Matrix::const_iterator firstRow_;
  118       //! @brief Iterator pointing to the current row index.
  119       typename RowIndexSet::const_iterator pos_;
  120     };
  121 
  122     //! @brief Get the row iterator at the first row.
  123     const_iterator begin() const
  124     {
  125       return const_iterator(m_.begin(), s_.begin());
  126     }
  127     //! @brief Get the row iterator at the end of all rows.
  128     const_iterator end() const
  129     {
  130       return const_iterator(m_.begin(), s_.end());
  131     }
  132 
  133   private:
  134     //! @brief The matrix for which we manage the row subset.
  135     const Matrix& m_;
  136     //! @brief The set of row indices we manage.
  137     const RowIndexSet& s_;
  138   };
  139 
  140   /**
  141    * @brief Inititializer for the ColCompMatrix
  142    * as needed by OverlappingSchwarz
  143    * @tparam M the matrix type
  144    * @tparam I the internal index type
  145    */
  146   template<class M, class I = int>
  147   class ColCompMatrixInitializer;
  148 
  149   template<class M, class X, class TM, class TD, class T1>
  150   class SeqOverlappingSchwarz;
  151 
  152   template<class T, bool flag>
  153   struct SeqOverlappingSchwarzAssemblerHelper;
  154 
  155   /**
  156    * @brief Utility class for converting an ISTL Matrix into a column-compressed matrix
  157    * @tparam M The matrix type
  158    * @tparam I the internal index type
  159    */
  160   template<class Mat, class I = int>
  161   class ColCompMatrix
  162   {
  163     friend class ColCompMatrixInitializer<Mat, I>;
  164 
  165     using B = typename Mat::field_type;
  166 
  167   public:
  168     /** @brief The type of the matrix to convert. */
  169     using Matrix = Mat;
  170 
  171     friend struct SeqOverlappingSchwarzAssemblerHelper<ColCompMatrix<Matrix>, true>;
  172 
  173     typedef typename Matrix::size_type size_type;
  174 
  175     using Index = I;
  176 
  177     /**
  178      * @brief Constructor that initializes the data.
  179      * @param mat The matrix to convert.
  180      */
  181     explicit ColCompMatrix(const Matrix& mat);
  182 
  183     ColCompMatrix();
  184 
  185     /** @brief Destructor */
  186     virtual ~ColCompMatrix();
  187 
  188     /**
  189      * @brief Get the number of rows.
  190      * @return  The number of rows.
  191      */
  192     size_type N() const
  193     {
  194       return N_;
  195     }
  196 
  197     size_type nnz() const
  198     {
  199       // TODO: The following code assumes that the blocks are dense
  200       // and that they all have the same dimensions.
  201       typename Matrix::block_type dummy;
  202       const auto n = MatrixDimension<typename Matrix::block_type>::rowdim(dummy);
  203       const auto m = MatrixDimension<typename Matrix::block_type>::coldim(dummy);
  204       return Nnz_/n/m;
  205     }
  206 
  207     /**
  208      * @brief Get the number of columns.
  209      * @return  The number of columns.
  210      */
  211     size_type M() const
  212     {
  213       return M_;
  214     }
  215 
  216     B* getValues() const
  217     {
  218       return values;
  219     }
  220 
  221     Index* getRowIndex() const
  222     {
  223       return rowindex;
  224     }
  225 
  226     Index* getColStart() const
  227     {
  228       return colstart;
  229     }
  230 
  231     ColCompMatrix& operator=(const Matrix& mat);
  232     ColCompMatrix& operator=(const ColCompMatrix& mat);
  233 
  234     /**
  235      * @brief Initialize data from a given set of matrix rows and columns
  236      * @param mat the matrix with the values
  237      * @param mrs The set of row (and column) indices to remove
  238      */
  239     virtual void setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs);
  240     /** @brief free allocated space. */
  241     virtual void free();
  242 
  243     /** @brief Initialize data from given matrix. */
  244     virtual void setMatrix(const Matrix& mat);
  245 
  246   public:
  247     size_type N_, M_, Nnz_;
  248     B* values;
  249     Index* rowindex;
  250     Index* colstart;
  251   };
  252 
  253   template<class M, class I>
  254   class ColCompMatrixInitializer
  255   {
  256     template<class IList, class S, class D>
  257     friend class OverlappingSchwarzInitializer;
  258   public:
  259     using Matrix = M;
  260     using Index = I;
  261     typedef Dune::ColCompMatrix<Matrix, Index> ColCompMatrix;
  262     typedef typename Matrix::row_type::const_iterator CIter;
  263     typedef typename Matrix::size_type size_type;
  264 
  265     /** \brief Constructor for scalar-valued matrices
  266      *
  267      * \tparam Block Dummy parameter to make SFINAE work
  268      */
  269     template <class Block=typename M::block_type>
  270     ColCompMatrixInitializer(ColCompMatrix& lum,
  271                              typename std::enable_if_t<Dune::IsNumber<Block>::value>* sfinae = nullptr);
  272 
  273     /** \brief Constructor for dense matrix-valued matrices
  274      *
  275      * \tparam Block Dummy parameter to make SFINAE work
  276      */
  277     template <class Block=typename M::block_type>
  278     ColCompMatrixInitializer(ColCompMatrix& lum,
  279                              typename std::enable_if_t<!Dune::IsNumber<Block>::value>* sfinae = nullptr);
  280 
  281     ColCompMatrixInitializer();
  282 
  283     virtual ~ColCompMatrixInitializer()
  284     {}
  285 
  286     template<typename Iter>
  287     void addRowNnz(const Iter& row) const;
  288 
  289     template<typename Iter, typename FullMatrixIndex>
  290     void addRowNnz(const Iter& row, const std::set<FullMatrixIndex>& indices) const;
  291 
  292     template<typename Iter, typename SubMatrixIndex>
  293     void addRowNnz(const Iter& row, const std::vector<SubMatrixIndex>& indices) const;
  294 
  295     void allocate();
  296 
  297     template<typename Iter>
  298     void countEntries(const Iter& row, const CIter& col) const;
  299 
  300     void countEntries(size_type colidx) const;
  301 
  302     void calcColstart() const;
  303 
  304     template<typename Iter>
  305     void copyValue(const Iter& row, const CIter& col) const;
  306 
  307     void copyValue(const CIter& col, size_type rowindex, size_type colidx) const;
  308 
  309     virtual void createMatrix() const;
  310 
  311   protected:
  312 
  313     void allocateMatrixStorage() const;
  314 
  315     void allocateMarker();
  316 
  317     ColCompMatrix* mat;
  318     size_type cols;
  319 
  320     // Number of rows/columns of the matrix entries
  321     // (assumed to be scalars or dense matrices)
  322     size_type n, m;
  323 
  324     mutable std::vector<size_type> marker;
  325   };
  326 
  327   template<class M, class I>
  328   template <class Block>
  329   ColCompMatrixInitializer<M, I>::ColCompMatrixInitializer(ColCompMatrix& mat_,typename std::enable_if_t<Dune::IsNumber<Block>::value>* sfinae)
  330     : mat(&mat_), cols(mat_.M())
  331   {
  332     n = 1;
  333     m = 1;
  334 
  335     mat->Nnz_=0;
  336   }
  337 
  338   template<class M, class I>
  339   template <class Block>
  340   ColCompMatrixInitializer<M, I>::ColCompMatrixInitializer(ColCompMatrix& mat_,typename std::enable_if_t<!Dune::IsNumber<Block>::value>* sfinae)
  341     : mat(&mat_), cols(mat_.M())
  342   {
  343     // WARNING: This assumes that all blocks are dense and identical
  344     n = M::block_type::rows;
  345     m = M::block_type::cols;
  346 
  347     mat->Nnz_=0;
  348   }
  349 
  350   template<class M, class I>
  351   ColCompMatrixInitializer<M, I>::ColCompMatrixInitializer()
  352     : mat(0), cols(0), n(0), m(0)
  353   {}
  354 
  355   template<class M, class I>
  356   template<typename Iter>
  357   void ColCompMatrixInitializer<M, I>::addRowNnz(const Iter& row) const
  358   {
  359     mat->Nnz_+=row->getsize();
  360   }
  361 
  362   template<class M, class I>
  363   template<typename Iter, typename FullMatrixIndex>
  364   void ColCompMatrixInitializer<M, I>::addRowNnz(const Iter& row,
  365                                                                             const std::set<FullMatrixIndex>& indices) const
  366   {
  367     typedef typename  Iter::value_type::const_iterator RIter;
  368     typedef typename std::set<FullMatrixIndex>::const_iterator MIter;
  369     MIter siter =indices.begin();
  370     for(RIter entry=row->begin(); entry!=row->end(); ++entry)
  371     {
  372       for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
  373       if(siter==indices.end())
  374         break;
  375       if(*siter==entry.index())
  376         // index is in subdomain
  377         ++mat->Nnz_;
  378     }
  379   }
  380 
  381   template<class M, class I>
  382   template<typename Iter, typename SubMatrixIndex>
  383   void ColCompMatrixInitializer<M, I>::addRowNnz(const Iter& row,
  384                                                                             const std::vector<SubMatrixIndex>& indices) const
  385   {
  386     using RIter = typename Iter::value_type::const_iterator;
  387     for(RIter entry=row->begin(); entry!=row->end(); ++entry)
  388       if (indices[entry.index()]!=std::numeric_limits<SubMatrixIndex>::max())
  389           ++mat->Nnz_;
  390   }
  391 
  392   template<class M, class I>
  393   void ColCompMatrixInitializer<M, I>::allocate()
  394   {
  395     allocateMatrixStorage();
  396     allocateMarker();
  397   }
  398 
  399   template<class M, class I>
  400   void ColCompMatrixInitializer<M, I>::allocateMatrixStorage() const
  401   {
  402     mat->Nnz_*=n*m;
  403     // initialize data
  404     mat->values=new typename M::field_type[mat->Nnz_];
  405     mat->rowindex=new I[mat->Nnz_];
  406     mat->colstart=new I[cols+1];
  407   }
  408 
  409   template<class M, class I>
  410   void ColCompMatrixInitializer<M, I>::allocateMarker()
  411   {
  412     marker.resize(cols);
  413     std::fill(marker.begin(), marker.end(), 0);
  414   }
  415 
  416   template<class M, class I>
  417   template<typename Iter>
  418   void ColCompMatrixInitializer<M, I>::countEntries(const Iter& row, const CIter& col) const
  419   {
  420     DUNE_UNUSED_PARAMETER(row);
  421     countEntries(col.index());
  422   }
  423 
  424   template<class M, class I>
  425   void ColCompMatrixInitializer<M, I>::countEntries(size_type colindex) const
  426   {
  427     for(size_type i=0; i < m; ++i)
  428     {
  429       assert(colindex*m+i<cols);
  430       marker[colindex*m+i]+=n;
  431     }
  432   }
  433 
  434   template<class M, class I>
  435   void ColCompMatrixInitializer<M, I>::calcColstart() const
  436   {
  437     mat->colstart[0]=0;
  438     for(size_type i=0; i < cols; ++i) {
  439       assert(i<cols);
  440       mat->colstart[i+1]=mat->colstart[i]+marker[i];
  441       marker[i]=mat->colstart[i];
  442     }
  443   }
  444 
  445   template<class M, class I>
  446   template<typename Iter>
  447   void ColCompMatrixInitializer<M, I>::copyValue(const Iter& row, const CIter& col) const
  448   {
  449     copyValue(col, row.index(), col.index());
  450   }
  451 
  452   template<class M, class I>
  453   void ColCompMatrixInitializer<M, I>::copyValue(const CIter& col, size_type rowindex, size_type colindex) const
  454   {
  455     for(size_type i=0; i<n; i++) {
  456       for(size_type j=0; j<m; j++) {
  457         assert(colindex*m+j<cols-1 || (size_type)marker[colindex*m+j]<(size_type)mat->colstart[colindex*m+j+1]);
  458         assert((size_type)marker[colindex*m+j]<mat->Nnz_);
  459         mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
  460         mat->values[marker[colindex*m+j]]=Impl::asMatrix(*col)[i][j];
  461         ++marker[colindex*m+j]; // index for next entry in column
  462       }
  463     }
  464   }
  465 
  466   template<class M, class I>
  467   void ColCompMatrixInitializer<M, I>::createMatrix() const
  468   {
  469     marker.clear();
  470   }
  471 
  472   template<class F, class MRS>
  473   void copyToColCompMatrix(F& initializer, const MRS& mrs)
  474   {
  475     typedef typename MRS::const_iterator Iter;
  476     typedef typename  std::iterator_traits<Iter>::value_type::const_iterator CIter;
  477     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
  478       initializer.addRowNnz(row);
  479 
  480     initializer.allocate();
  481 
  482     for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
  483 
  484       for(CIter col=row->begin(); col != row->end(); ++col)
  485         initializer.countEntries(row, col);
  486     }
  487 
  488     initializer.calcColstart();
  489 
  490     for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
  491       for(CIter col=row->begin(); col != row->end(); ++col) {
  492         initializer.copyValue(row, col);
  493       }
  494 
  495     }
  496     initializer.createMatrix();
  497   }
  498 
  499   template<class F, class M,class S>
  500   void copyToColCompMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
  501   {
  502     typedef MatrixRowSubset<M,S> MRS;
  503     typedef typename MRS::RowIndexSet SIS;
  504     typedef typename SIS::const_iterator SIter;
  505     typedef typename MRS::const_iterator Iter;
  506     typedef typename std::iterator_traits<Iter>::value_type row_type;
  507     typedef typename row_type::const_iterator CIter;
  508 
  509     typedef typename MRS::Matrix::size_type size_type;
  510 
  511     // A vector containing the corresponding indices in
  512     // the to create submatrix.
  513     // If an entry is the maximum of size_type then this index will not appear in
  514     // the submatrix.
  515     std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
  516                                           std::numeric_limits<size_type>::max());
  517     size_type s=0;
  518     for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
  519       subMatrixIndex[*index]=s++;
  520 
  521     // Calculate upper Bound for nonzeros
  522     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
  523       initializer.addRowNnz(row, subMatrixIndex);
  524 
  525     initializer.allocate();
  526 
  527     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
  528       for(CIter col=row->begin(); col != row->end(); ++col) {
  529         if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
  530           // This column is in our subset (use submatrix column index)
  531           initializer.countEntries(subMatrixIndex[col.index()]);
  532       }
  533 
  534     initializer.calcColstart();
  535 
  536     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
  537       for(CIter col=row->begin(); col != row->end(); ++col) {
  538         if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
  539           // This value is in our submatrix -> copy (use submatrix indices
  540           initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
  541       }
  542     initializer.createMatrix();
  543   }
  544 
  545 #ifndef DOXYGEN
  546 
  547   template<class Mat, class I>
  548   ColCompMatrix<Mat, I>::ColCompMatrix()
  549     : N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
  550   {}
  551 
  552   template<class Mat, class I>
  553   ColCompMatrix<Mat, I>
  554   ::ColCompMatrix(const Matrix& mat)
  555   {
  556     // WARNING: This assumes that all blocks are dense and identical
  557     typename Matrix::block_type dummy;
  558     const auto n = MatrixDimension<typename Matrix::block_type>::rowdim(dummy);
  559     const auto m = MatrixDimension<typename Matrix::block_type>::coldim(dummy);
  560     N_ = n*mat.N();
  561     M_ = m*mat.N();
  562     Nnz_ = n*m*mat.nonzeroes();
  563   }
  564 
  565   template<class Mat, class I>
  566   ColCompMatrix<Mat, I>&
  567   ColCompMatrix<Mat, I>::operator=(const Matrix& mat)
  568   {
  569     if(N_+M_+Nnz_!=0)
  570       free();
  571     setMatrix(mat);
  572     return *this;
  573   }
  574 
  575   template<class Mat, class I>
  576   ColCompMatrix<Mat, I>&
  577   ColCompMatrix<Mat, I>::operator=(const ColCompMatrix& mat)
  578   {
  579     if(N_+M_+Nnz_!=0)
  580       free();
  581     N_=mat.N_;
  582     M_=mat.M_;
  583     Nnz_= mat.Nnz_;
  584     if(M_>0) {
  585       colstart=new size_type[M_+1];
  586       for(size_type i=0; i<=M_; ++i)
  587         colstart[i]=mat.colstart[i];
  588     }
  589 
  590     if(Nnz_>0) {
  591       values = new B[Nnz_];
  592       rowindex = new size_type[Nnz_];
  593 
  594       for(size_type i=0; i<Nnz_; ++i)
  595         values[i]=mat.values[i];
  596 
  597       for(size_type i=0; i<Nnz_; ++i)
  598         rowindex[i]=mat.rowindex[i];
  599     }
  600     return *this;
  601   }
  602 
  603   template<class Mat, class I>
  604   void ColCompMatrix<Mat, I>
  605   ::setMatrix(const Matrix& mat)
  606   {
  607     N_=MatrixDimension<Mat>::rowdim(mat);
  608     M_=MatrixDimension<Mat>::coldim(mat);
  609     ColCompMatrixInitializer<Mat, I> initializer(*this);
  610 
  611     copyToColCompMatrix(initializer, MatrixRowSet<Matrix>(mat));
  612   }
  613 
  614   template<class Mat, class I>
  615   void ColCompMatrix<Mat, I>
  616   ::setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs)
  617   {
  618     if(N_+M_+Nnz_!=0)
  619       free();
  620 
  621     N_=mrs.size()*MatrixDimension<Mat>::rowdim(mat) / mat.N();
  622     M_=mrs.size()*MatrixDimension<Mat>::coldim(mat) / mat.M();
  623     ColCompMatrixInitializer<Mat, I> initializer(*this);
  624 
  625     copyToColCompMatrix(initializer, MatrixRowSubset<Mat,std::set<std::size_t> >(mat,mrs));
  626   }
  627 
  628   template<class Mat, class I>
  629   ColCompMatrix<Mat, I>::~ColCompMatrix()
  630   {
  631     if(N_+M_+Nnz_!=0)
  632       free();
  633   }
  634 
  635   template<class Mat, class I>
  636   void ColCompMatrix<Mat, I>::free()
  637   {
  638     delete[] values;
  639     delete[] rowindex;
  640     delete[] colstart;
  641     N_ = 0;
  642     M_ = 0;
  643     Nnz_ = 0;
  644   }
  645 
  646 #endif // DOXYGEN
  647 
  648 }
  649 #endif