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