dune-istl  2.8.0
bccsmatrixinitializer.hh
Go to the documentation of this file.
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_BCCSMATRIX_INITIALIZER_HH
4 #define DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
5 
6 #include <limits>
7 #include <set>
8 
9 #include <dune/common/typetraits.hh>
10 #include <dune/common/scalarmatrixview.hh>
11 
12 #include <dune/istl/bccsmatrix.hh>
13 
14 namespace Dune
15 {
16  template<class I, class S, class D>
17  class OverlappingSchwarzInitializer;
18 }
19 
20 namespace Dune::ISTL::Impl
21 {
29  template<class M, class S>
30  class MatrixRowSubset
31  {
32  public:
34  typedef M Matrix;
36  typedef S RowIndexSet;
37 
43  MatrixRowSubset(const Matrix& m, const RowIndexSet& s)
44  : m_(m), s_(s)
45  {}
46 
47  const Matrix& matrix() const
48  {
49  return m_;
50  }
51 
52  const RowIndexSet& rowIndexSet() const
53  {
54  return s_;
55  }
56 
58  class const_iterator
59  : public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
60  {
61  public:
62  const_iterator(typename Matrix::const_iterator firstRow,
63  typename RowIndexSet::const_iterator pos)
64  : firstRow_(firstRow), pos_(pos)
65  {}
66 
67 
68  const typename Matrix::row_type& dereference() const
69  {
70  return *(firstRow_+ *pos_);
71  }
72  bool equals(const const_iterator& o) const
73  {
74  return pos_==o.pos_;
75  }
76  void increment()
77  {
78  ++pos_;
79  }
80  typename RowIndexSet::value_type index() const
81  {
82  return *pos_;
83  }
84 
85  private:
87  typename Matrix::const_iterator firstRow_;
89  typename RowIndexSet::const_iterator pos_;
90  };
91 
93  const_iterator begin() const
94  {
95  return const_iterator(m_.begin(), s_.begin());
96  }
98  const_iterator end() const
99  {
100  return const_iterator(m_.begin(), s_.end());
101  }
102 
103  private:
105  const Matrix& m_;
107  const RowIndexSet& s_;
108  };
109 
116  template<class M, class I = typename M::size_type>
117  class BCCSMatrixInitializer
118  {
119  template<class IList, class S, class D>
121  public:
122  using Matrix = M;
123  using Index = I;
124  typedef Dune::ISTL::Impl::BCCSMatrix<typename Matrix::field_type, I> OutputMatrix;
125  typedef typename Matrix::size_type size_type;
126 
129  BCCSMatrixInitializer(OutputMatrix& mat_)
130  : mat(&mat_), cols(mat_.M())
131  {
132  if constexpr (Dune::IsNumber<typename M::block_type>::value)
133  {
134  n = m = 1;
135  }
136  else
137  {
138  // WARNING: This assumes that all blocks are dense and identical
139  n = M::block_type::rows;
140  m = M::block_type::cols;
141  }
142 
143  mat->Nnz_=0;
144  }
145 
146  BCCSMatrixInitializer()
147  : mat(0), cols(0), n(0), m(0)
148  {}
149 
150  virtual ~BCCSMatrixInitializer()
151  {}
152 
153  template<typename Iter>
154  void addRowNnz(const Iter& row) const
155  {
156  mat->Nnz_+=row->getsize();
157  }
158 
159  template<typename Iter, typename FullMatrixIndex>
160  void addRowNnz(const Iter& row, const std::set<FullMatrixIndex>& indices) const
161  {
162  auto siter =indices.begin();
163  for (auto entry=row->begin(); entry!=row->end(); ++entry)
164  {
165  for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
166  if(siter==indices.end())
167  break;
168  if(*siter==entry.index())
169  // index is in subdomain
170  ++mat->Nnz_;
171  }
172  }
173 
174  template<typename Iter, typename SubMatrixIndex>
175  void addRowNnz(const Iter& row, const std::vector<SubMatrixIndex>& indices) const
176  {
177  for (auto entry=row->begin(); entry!=row->end(); ++entry)
178  if (indices[entry.index()]!=std::numeric_limits<SubMatrixIndex>::max())
179  ++mat->Nnz_;
180  }
181 
182  void allocate()
183  {
184  allocateMatrixStorage();
185  allocateMarker();
186  }
187 
188  template<typename Iter, typename CIter>
189  void countEntries([[maybe_unused]] const Iter& row, const CIter& col) const
190  {
191  countEntries(col.index());
192  }
193 
194  void countEntries(size_type colindex) const
195  {
196  for(size_type i=0; i < m; ++i)
197  {
198  assert(colindex*m+i<cols);
199  marker[colindex*m+i]+=n;
200  }
201  }
202 
203  void calcColstart() const
204  {
205  mat->colstart[0]=0;
206  for(size_type i=0; i < cols; ++i) {
207  assert(i<cols);
208  mat->colstart[i+1]=mat->colstart[i]+marker[i];
209  marker[i]=mat->colstart[i];
210  }
211  }
212 
213  template<typename Iter, typename CIter>
214  void copyValue(const Iter& row, const CIter& col) const
215  {
216  copyValue(col, row.index(), col.index());
217  }
218 
219  template<typename CIter>
220  void copyValue(const CIter& col, size_type rowindex, size_type colindex) const
221  {
222  for(size_type i=0; i<n; i++) {
223  for(size_type j=0; j<m; j++) {
224  assert(colindex*m+j<cols-1 || (size_type)marker[colindex*m+j]<(size_type)mat->colstart[colindex*m+j+1]);
225  assert((size_type)marker[colindex*m+j]<mat->Nnz_);
226  mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
227  mat->values[marker[colindex*m+j]] = Dune::Impl::asMatrix(*col)[i][j];
228  ++marker[colindex*m+j]; // index for next entry in column
229  }
230  }
231  }
232 
233  virtual void createMatrix() const
234  {
235  marker.clear();
236  }
237 
238  protected:
239 
240  void allocateMatrixStorage() const
241  {
242  mat->Nnz_*=n*m;
243  // initialize data
244  mat->values=new typename M::field_type[mat->Nnz_];
245  mat->rowindex=new I[mat->Nnz_];
246  mat->colstart=new I[cols+1];
247  }
248 
249  void allocateMarker()
250  {
251  marker.resize(cols);
252  std::fill(marker.begin(), marker.end(), 0);
253  }
254 
255  OutputMatrix* mat;
256  size_type cols;
257 
258  // Number of rows/columns of the matrix entries
259  // (assumed to be scalars or dense matrices)
260  size_type n, m;
261 
262  mutable std::vector<size_type> marker;
263  };
264 
265  template<class F, class Matrix>
266  void copyToBCCSMatrix(F& initializer, const Matrix& matrix)
267  {
268  for (auto row=matrix.begin(); row!= matrix.end(); ++row)
269  initializer.addRowNnz(row);
270 
271  initializer.allocate();
272 
273  for (auto row=matrix.begin(); row!= matrix.end(); ++row) {
274 
275  for (auto col=row->begin(); col != row->end(); ++col)
276  initializer.countEntries(row, col);
277  }
278 
279  initializer.calcColstart();
280 
281  for (auto row=matrix.begin(); row!= matrix.end(); ++row) {
282  for (auto col=row->begin(); col != row->end(); ++col) {
283  initializer.copyValue(row, col);
284  }
285 
286  }
287  initializer.createMatrix();
288  }
289 
290  template<class F, class M,class S>
291  void copyToBCCSMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
292  {
293  typedef MatrixRowSubset<M,S> MRS;
294  typedef typename MRS::RowIndexSet SIS;
295  typedef typename SIS::const_iterator SIter;
296  typedef typename MRS::const_iterator Iter;
297  typedef typename std::iterator_traits<Iter>::value_type row_type;
298  typedef typename row_type::const_iterator CIter;
299 
300  typedef typename MRS::Matrix::size_type size_type;
301 
302  // A vector containing the corresponding indices in
303  // the to create submatrix.
304  // If an entry is the maximum of size_type then this index will not appear in
305  // the submatrix.
306  std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
307  std::numeric_limits<size_type>::max());
308  size_type s=0;
309  for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
310  subMatrixIndex[*index]=s++;
311 
312  // Calculate upper Bound for nonzeros
313  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
314  initializer.addRowNnz(row, subMatrixIndex);
315 
316  initializer.allocate();
317 
318  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
319  for(CIter col=row->begin(); col != row->end(); ++col) {
320  if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
321  // This column is in our subset (use submatrix column index)
322  initializer.countEntries(subMatrixIndex[col.index()]);
323  }
324 
325  initializer.calcColstart();
326 
327  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
328  for(CIter col=row->begin(); col != row->end(); ++col) {
329  if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
330  // This value is in our submatrix -> copy (use submatrix indices
331  initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
332  }
333  initializer.createMatrix();
334  }
335 
336 }
337 #endif
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
void addRowNnz(const Iter &row)
Definition: overlappingschwarz.hh:893
void calcColstart() const
Definition: overlappingschwarz.hh:924
void copyValue(const Iter &row, const CIter &col) const
Definition: overlappingschwarz.hh:931
void createMatrix() const
Definition: overlappingschwarz.hh:945
void allocate()
Definition: overlappingschwarz.hh:903
std::size_t countEntries(const BlockVector< T, A > &vector)
Definition: matrixmarket.hh:1060
Definition: allocator.hh:9
Initializer for SuperLU Matrices representing the subdomains.
Definition: overlappingschwarz.hh:45
Matrix::row_type::const_iterator CIter
Definition: overlappingschwarz.hh:54
Matrix::const_iterator Iter
Definition: overlappingschwarz.hh:53
IndexSet::size_type size_type
Definition: overlappingschwarz.hh:57
AtomInitializer::Matrix Matrix
Definition: overlappingschwarz.hh:52
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:575
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:572