DUNE PDELab (git)

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