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_COLCOMPSPMATRIX_HH 4 #define DUNE_COLCOMPSPMATRIX_HH 5 6 #if HAVE_DUNE_ISTL 7 #include <dune/istl/bccsmatrixinitializer.hh> 8 9 // if the original deprecated header has not been included 10 // it is ok to declare the ColCompMatrix here 11 #ifndef DUNE_ISTL_COLCOMPMATRIX_HH 12 namespace Dune 13 { 14 template <class M, class RowIndex = int> 15 struct ColCompMatrix : public ISTL::Impl::BCCSMatrix< typename M::field_type, RowIndex>{}; 16 } // namespace Dune 17 #endif 18 19 #else // #if HAVE_DUNE_ISTL 20 namespace Dune 21 { 22 template<class M,class RowIndex=int> 23 struct ColCompMatrix {}; 24 } // namespace Dune 25 #endif // #if HAVE_DUNE_ISTL 26 27 #include <dune/fem/operator/matrix/spmatrix.hh> 28 29 #include <vector> 30 31 namespace Dune 32 { 33 /** 34 * @brief Converter for SparseRowMatrix to column-compressed matrix. 35 * Specialization for SparseRowMatrix 36 */ 37 template <class T, class IndexT,class ValuesVector, class IndicesVector,class RowIndex> 38 class ColCompMatrix< Fem::SparseRowMatrix<T,IndexT,ValuesVector,IndicesVector>, RowIndex > 39 { 40 public: 41 /** @brief The type of the matrix converted. */ 42 typedef ColCompMatrix< Fem::SparseRowMatrix<T,IndexT,ValuesVector,IndicesVector>> ThisType; 43 /** @brief The type of the matrix to convert. */ 44 typedef Fem::SparseRowMatrix<T,IndexT,ValuesVector,IndicesVector> Matrix; 45 46 typedef typename Matrix::size_type size_type; 47 48 typedef RowIndex RowIndexType; 49 50 /** 51 * @brief Constructor that initializes the data. 52 * @param mat The matrix to convert. 53 */ ColCompMatrix(const Matrix & mat)54 explicit ColCompMatrix(const Matrix& mat) 55 { 56 setMatrix(mat); 57 } 58 59 /** @brief Empty constructor. */ ColCompMatrix()60 ColCompMatrix() : 61 N_(0), M_(0), Nnz_(0), values_(0), rowindex_(0), colstart_(0) 62 {} 63 64 /** @brief Destructor. */ ~ColCompMatrix()65 virtual ~ColCompMatrix() 66 { 67 if(N_+M_+Nnz_ != 0) 68 free(); 69 } 70 71 /** @brief Get the number of rows. */ N() const72 size_type N() const 73 { 74 return N_; 75 } 76 77 /** @brief Get the number of non zero entries. */ nnz() const78 size_type nnz() const 79 { 80 return Nnz_; 81 } 82 83 /** @brief Get the number of columns. */ M() const84 size_type M() const 85 { 86 return M_; 87 } 88 89 /** @brief Get the non-zero entries of the matrix. */ getValues() const90 T* getValues() const 91 { 92 return values_; 93 } 94 95 /** @brief Get the row indices of the non-zero entries of the matrix. */ getRowIndex() const96 RowIndexType* getRowIndex() const 97 { 98 return rowindex_; 99 } 100 101 /** @brief Get the column start indices. */ getColStart() const102 RowIndexType* getColStart() const 103 { 104 return colstart_; 105 } 106 operator =(const Matrix & mat)107 ThisType& operator=(const Matrix& mat) 108 { 109 if(N_+M_+Nnz_ != 0) 110 free(); 111 setMatrix(mat); 112 return *this; 113 } 114 operator =(const ThisType & mat)115 ThisType& operator=(const ThisType& mat) 116 { 117 if(N_+M_+Nnz_ != 0) 118 free(); 119 N_=mat.N_; 120 M_=mat.M_; 121 Nnz_=mat.Nnz_; 122 if(M_>0) 123 { 124 colstart_=new RowIndexType[M_+1]; 125 for(size_type i=0; i<=M_; ++i) 126 colstart_[i]=mat.colstart[i]; 127 } 128 if(Nnz_>0) 129 { 130 values_ = new T[Nnz_]; 131 rowindex_ = new RowIndexType[Nnz_]; 132 for(size_type i=0; i<Nnz_; ++i) 133 values_[i]=mat.values[i]; 134 for(size_type i=0; i<Nnz_; ++i) 135 rowindex_[i]=mat.rowindex[i]; 136 } 137 return *this; 138 } 139 140 /** @brief Free allocated space. */ free()141 void free() 142 { 143 delete[] values_; 144 delete[] rowindex_; 145 delete[] colstart_; 146 N_=0; 147 M_=0; 148 Nnz_=0; 149 } 150 151 /** @brief Initialize data from given matrix. */ setMatrix(const Matrix & mat)152 virtual void setMatrix(const Matrix& mat) 153 { 154 N_=mat.rows(); 155 M_=mat.cols(); 156 157 // count the number of nonzeros per column 158 colstart_= new RowIndexType[M_+1]; 159 for(size_type i=0;i!=(M_+1);++i) 160 colstart_[i]=0; 161 162 Nnz_ = 0; 163 for(size_type row=0; row < N_; ++row) 164 { 165 const size_type endRow = mat.endRow( row ); 166 for( size_type col = mat.startRow( row ); col < endRow; ++col ) 167 { 168 const auto pairIdx(mat.realValue( col )); 169 if( pairIdx.second!=Matrix::defaultCol ) 170 { 171 ++(colstart_[pairIdx.second+1]); 172 ++Nnz_; 173 } 174 } 175 } 176 177 // compute the starting positions and compute colstart 178 std::vector<int> tempPos(M_,0); 179 for(size_type i=1;i!=(M_+1);++i) 180 { 181 colstart_[i]+=colstart_[i-1]; 182 tempPos[i-1]=colstart_[i-1]; 183 } 184 185 // fill the values and the index arrays 186 values_=new T[Nnz_]; 187 rowindex_=new RowIndexType[Nnz_]; 188 for(size_type row = 0; row < N_ ; ++ row) 189 { 190 const size_type endRow = mat.endRow( row ); 191 for( size_type col = mat.startRow( row ); col < endRow; ++col ) 192 { 193 const auto pairIdx(mat.realValue( col )); 194 if(pairIdx.second!=Matrix::defaultCol) 195 { 196 values_[tempPos[pairIdx.second]] = pairIdx.first; 197 rowindex_[tempPos[pairIdx.second]] = row ; 198 ++(tempPos[pairIdx.second]); 199 } 200 } 201 } 202 203 } 204 205 private: 206 size_type N_; 207 size_type M_; 208 size_type Nnz_; 209 T* values_; 210 RowIndexType* rowindex_; 211 RowIndexType* colstart_; 212 }; 213 214 } 215 #endif // #ifndef DUNE_COLCOMPSPMATRIX_HH 216 217