Dune Core Modules (2.8.0)

bccsmatrixinitializer.hh
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
11
12#include <dune/istl/bccsmatrix.hh>
13
14namespace Dune
15{
16 template<class I, class S, class D>
17 class OverlappingSchwarzInitializer;
18}
19
20namespace 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 {
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(),
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
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
Initializer for SuperLU Matrices representing the subdomains.
Definition: overlappingschwarz.hh:45
Traits for type conversions and type information.
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:400
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:11
Implements a scalar matrix view wrapper around an existing scalar.
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)