dune-pdelab  2.4-dev
bcrspattern.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_BACKEND_ISTL_BCRSPATTERN_HH
3 #define DUNE_PDELAB_BACKEND_ISTL_BCRSPATTERN_HH
4 
5 #include <utility>
6 #include <vector>
7 #include <algorithm>
8 #include <set>
9 
10 #include <dune/common/iteratorfacades.hh>
11 
16 
17 namespace Dune {
18  namespace PDELab {
19  namespace istl {
20 
22 
42  template<typename RowOrdering, typename ColOrdering>
44  {
45 
46  public:
47 
49  typedef typename RowOrdering::Traits::size_type size_type;
50 
52  typedef void SubPattern;
53 
54  private:
55 
57  static const size_type empty = ~static_cast<size_type>(0);
58 
60 
69  struct PaddedColumnCriterion
70  {
71 
72  bool operator()(size_type k) const
73  {
74  return k == _j || k == empty;
75  }
76 
77  PaddedColumnCriterion(size_type j)
78  : _j(j)
79  {}
80 
81  const size_type _j;
82 
83  };
84 
85 
86  typedef typename std::vector<size_type>::iterator IndicesIterator;
87  typedef typename std::set<std::pair<size_type,size_type> >::iterator OverflowIterator;
88 
89  typedef typename std::vector<size_type>::const_iterator ConstIndicesIterator;
90  typedef typename std::set<std::pair<size_type,size_type> >::const_iterator ConstOverflowIterator;
91 
92  public:
93 
95  template<typename RI, typename CI>
96  void add_link(const RI& ri, const CI& ci)
97  {
98  // extract block indices for current level
99  size_type i = ri.back();
100  size_type j = ci.back();
101 
102  IndicesIterator start = _indices.begin();
103  IndicesIterator begin = start + _entries_per_row*i;
104  IndicesIterator end = start + _entries_per_row*(i+1);
105 
106  // Does the entry (i,j) already exist?
107  IndicesIterator it = std::find_if(begin,end,PaddedColumnCriterion(j));
108  if (it != end)
109  {
110  // Yes, just write out j. This does the right thing regardless of whether
111  // it points at the correct column value or at an empty entry.
112  *it = j;
113  }
114  else
115  {
116  // The row is already full -> spill into map
117  _overflow.insert(std::make_pair(i,j));
118  }
119  }
120 
121 #ifndef DOXYGEN
122 
123  // implementation detail for nested patterns
124  template<typename RI, typename CI>
125  void recursive_add_entry(const RI& ri, const CI& ci)
126  {
127  add_link(ri,ci);
128  }
129 
130 #endif //DOXYGEN
131 
133  template<typename I>
134  void sizes(I rit) const
135  {
136  ConstIndicesIterator it = _indices.begin();
137  ConstIndicesIterator end = _indices.begin() + _entries_per_row;
138  ConstOverflowIterator oit = _overflow.begin();
139  ConstOverflowIterator oend = _overflow.end();
140  for (size_type i = 0; i < _row_ordering.blockCount(); ++i, ++rit, end+=_entries_per_row)
141  {
142  size_type s = 0;
143  // count non-empty column entries, break when first empty one is found.
144  for (; it != end; ++it)
145  {
146  if (*it == empty)
147  break;
148  ++s;
149  }
150  it = end;
151  // add overflow entries
152  for (; oit != oend && oit->first == i; ++oit)
153  ++s;
154  *rit = s;
155  }
156  }
157 
159  std::vector<size_type> sizes() const
160  {
161  std::vector<size_type> r(_row_ordering.blockCount());
162  sizes(r.begin());
163  return std::move(r);
164  }
165 
167  struct iterator
168  : public ForwardIteratorFacade<iterator, const size_type>
169  {
170 
171 #ifndef DOXYGEN
172 
173  const size_type& dereference() const
174  {
175  if (_in_overflow)
176  return _oit->second;
177  else
178  return *_it;
179  }
180 
181  void increment()
182  {
183  if (_in_overflow)
184  {
185  if (_oit != _oend)
186  ++_oit;
187  if (_oit->first == _row)
188  return;
189  // we have exhausted the row, invalidate iterator
190  _at_end = true;
191  }
192  else
193  {
194  if (_it != _end)
195  ++_it;
196  if (_it == _end || *_it == empty)
197  {
198  _in_overflow = true;
199  // we have exhausted the row, invalidate iterator
200  if (_oit == _oend || _oit->first > _row)
201  _at_end = true;
202  }
203  }
204  }
205 
206  bool equals(const iterator& other) const
207  {
208  if (_row != other._row)
209  return false;
210  if (_at_end || other._at_end)
211  return _at_end && other._at_end;
212  if (_in_overflow)
213  return _oit == other._oit;
214  else
215  return _it == other._it;
216  }
217 
218  iterator(const BCRSPattern& p, size_type row, bool at_end)
219  : _row(row)
220  , _in_overflow(false)
221  , _at_end(at_end)
222  , _it(p._indices.begin() + row * p._entries_per_row)
223  , _end(p._indices.begin() + (row+1) * p._entries_per_row)
224  , _oit(p._overflow.lower_bound(std::make_pair(row,0)))
225  , _oend(p._overflow.end())
226  {
227  // catch corner case with completely empty row
228  if ((!_at_end) && (_it == _end || *_it == empty))
229  {
230  _in_overflow = true;
231  _at_end = _oit == _oend || _oit->first != _row;
232  }
233  }
234 
235  size_type _row;
236  bool _in_overflow;
237  bool _at_end;
238  typename std::vector<size_type>::const_iterator _it;
239  typename std::vector<size_type>::const_iterator _end;
240  typename std::set<std::pair<size_type,size_type> >::const_iterator _oit;
241  const typename std::set<std::pair<size_type,size_type> >::const_iterator _oend;
242 
243 #endif // DOXYGEN
244 
245  };
246 
248  iterator begin(size_type i) const
249  {
250  return iterator(*this,i,false);
251  }
252 
254  iterator end(size_type i) const
255  {
256  return iterator(*this,i,true);
257  }
258 
260 
265  BCRSPattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering, size_type entries_per_row)
266  : _row_ordering(row_ordering)
267  , _col_ordering(col_ordering)
268  , _entries_per_row(entries_per_row)
269  , _indices(row_ordering.blockCount()*entries_per_row,size_type(empty))
270  {}
271 
272  const RowOrdering& rowOrdering() const
273  {
274  return _row_ordering;
275  }
276 
277  const ColOrdering& colOrdering() const
278  {
279  return _row_ordering;
280  }
281 
283 
289  void clear()
290  {
291  _indices = std::vector<size_type>();
292  _overflow = std::set<std::pair<size_type,size_type> >();
293  }
294 
295  size_type entriesPerRow() const
296  {
297  return _entries_per_row;
298  }
299 
300  size_type overflowCount() const
301  {
302  return _overflow.size();
303  }
304 
305  private:
306 
307  const RowOrdering& _row_ordering;
308  const ColOrdering& _col_ordering;
309  const size_type _entries_per_row;
310 
311  std::vector<size_type> _indices;
312  std::set<std::pair<size_type,size_type> > _overflow;
313 
314  };
315 
316 
318 
323  template<typename RowOrdering, typename ColOrdering, typename SubPattern_>
325  {
326 
327  public:
328 
330  typedef SubPattern_ SubPattern;
331 
333  typedef typename SubPattern::size_type size_type;
334 
336 
340  template<typename RI, typename CI>
341  void add_link(const RI& ri, const CI& ci)
342  {
343  recursive_add_entry(ri.view(),ci.view());
344  }
345 
346 #ifndef DOXYGEN
347 
348  template<typename RI, typename CI>
349  void recursive_add_entry(const RI& ri, const CI& ci)
350  {
351  _sub_patterns[ri.back() * _col_ordering.blockCount() + ci.back()].recursive_add_entry(ri.back_popped(),ci.back_popped());
352  }
353 
354 #endif // DOXYGEN
355 
356  template<typename EntriesPerRow>
357  NestedPattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering, const EntriesPerRow& entries_per_row)
358  : _row_ordering(row_ordering)
359  , _col_ordering(col_ordering)
360  {
361  size_type rows = row_ordering.blockCount();
362  size_type cols = col_ordering.blockCount();
363  for (size_type i = 0; i < rows; ++i)
364  for (size_type j = 0; j < cols; ++j)
365  _sub_patterns.push_back(
366  SubPattern(
367  _row_ordering.childOrdering(i),
368  _col_ordering.childOrdering(j),
369  entries_per_row[i][j]
370  )
371  );
372  }
373 
374  NestedPattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering, size_type entries_per_row)
375  : _row_ordering(row_ordering)
376  , _col_ordering(col_ordering)
377  {
378  size_type rows = row_ordering.blockCount();
379  size_type cols = col_ordering.blockCount();
380  for (size_type i = 0; i < rows; ++i)
381  for (size_type j = 0; j < cols; ++j)
382  _sub_patterns.push_back(
383  SubPattern(
384  _row_ordering.childOrdering(i),
385  _col_ordering.childOrdering(j),
386  entries_per_row
387  )
388  );
389  }
390 
392  SubPattern& subPattern(size_type i, size_type j)
393  {
394  return _sub_patterns[i * _col_ordering.blockCount() + j];
395  }
396 
397  private:
398 
399  const RowOrdering& _row_ordering;
400  const ColOrdering& _col_ordering;
401  std::vector<SubPattern> _sub_patterns;
402 
403  };
404 
405 
406  } // namespace istl
407  } // namespace PDELab
408 } // namespace Dune
409 
410 #endif // DUNE_PDELAB_BACKEND_ISTL_BCRSPATTERN_HH
SubPattern_ SubPattern
The pattern type used for each block.
Definition: bcrspattern.hh:330
void add_link(const RI &ri, const CI &ci)
Add a link between the row indicated by ri and the column indicated by ci.
Definition: bcrspattern.hh:96
size_type overflowCount() const
Definition: bcrspattern.hh:300
SubPattern::size_type size_type
size type used by NestedPattern.
Definition: bcrspattern.hh:333
void SubPattern
BCRSPattern cannot contain nested subpatterns. This entry is only required for TMP purposes...
Definition: bcrspattern.hh:52
Various tags for influencing backend behavior.
size_type entriesPerRow() const
Definition: bcrspattern.hh:295
Pattern builder for generic BCRS-like sparse matrices.
Definition: bcrspattern.hh:43
iterator end(size_type i) const
Returns an iterator past the last column index of row i.
Definition: bcrspattern.hh:254
BCRSPattern(const RowOrdering &row_ordering, const ColOrdering &col_ordering, size_type entries_per_row)
Constructs a BCRSPattern for the given pair of orderings and reserves space for the provided average ...
Definition: bcrspattern.hh:265
const RowOrdering & rowOrdering() const
Definition: bcrspattern.hh:272
Pattern builder for nested hierarchies of generic BCRS-like sparse matrices.
Definition: bcrspattern.hh:324
Definition: adaptivity.hh:27
void clear()
Discard all internal data.
Definition: bcrspattern.hh:289
OffsetIterator _oit
Definition: datahandleprovider.hh:76
std::vector< size_type > sizes() const
Returns a vector with the size of all rows in the pattern.
Definition: bcrspattern.hh:159
SubPattern & subPattern(size_type i, size_type j)
Returns the subpattern associated with block (i,j).
Definition: bcrspattern.hh:392
Iterator over all column indices for a given row, unique but in arbitrary order.
Definition: bcrspattern.hh:167
void add_link(const RI &ri, const CI &ci)
Add a link between the row indicated by ri and the column indicated by ci.
Definition: bcrspattern.hh:341
const std::string s
Definition: function.hh:1102
RowOrdering::Traits::size_type size_type
size type used by BCRSPattern.
Definition: bcrspattern.hh:49
const P & p
Definition: constraints.hh:146
NestedPattern(const RowOrdering &row_ordering, const ColOrdering &col_ordering, size_type entries_per_row)
Definition: bcrspattern.hh:374
NestedPattern(const RowOrdering &row_ordering, const ColOrdering &col_ordering, const EntriesPerRow &entries_per_row)
Definition: bcrspattern.hh:357
iterator begin(size_type i) const
Returns an iterator to the first column index of row i.
Definition: bcrspattern.hh:248
const ColOrdering & colOrdering() const
Definition: bcrspattern.hh:277
void sizes(I rit) const
Stream the sizes of all rows into the output iterator rit.
Definition: bcrspattern.hh:134