3#ifndef DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
4#define DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
12#include <dune/istl/bccsmatrix.hh>
16 template<
class I,
class S,
class D>
17 class OverlappingSchwarzInitializer;
20namespace Dune::ISTL::Impl
29 template<
class M,
class S>
36 typedef S RowIndexSet;
43 MatrixRowSubset(
const Matrix& m,
const RowIndexSet& s)
47 const Matrix& matrix()
const
52 const RowIndexSet& rowIndexSet()
const
59 :
public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
62 const_iterator(
typename Matrix::const_iterator firstRow,
63 typename RowIndexSet::const_iterator pos)
64 : firstRow_(firstRow), pos_(pos)
70 return *(firstRow_+ *pos_);
72 bool equals(
const const_iterator& o)
const
80 typename RowIndexSet::value_type index()
const
87 typename Matrix::const_iterator firstRow_;
89 typename RowIndexSet::const_iterator pos_;
93 const_iterator begin()
const
95 return const_iterator(m_.begin(), s_.begin());
98 const_iterator end()
const
100 return const_iterator(m_.begin(), s_.end());
107 const RowIndexSet& s_;
116 template<
class M,
class I =
typename M::
size_type>
117 class BCCSMatrixInitializer
119 template<
class IList,
class S,
class D>
124 typedef Dune::ISTL::Impl::BCCSMatrix<typename Matrix::field_type, I> OutputMatrix;
129 BCCSMatrixInitializer(OutputMatrix& mat_)
130 : mat(&mat_), cols(mat_.M())
139 n = M::block_type::rows;
140 m = M::block_type::cols;
146 BCCSMatrixInitializer()
147 : mat(0), cols(0), n(0), m(0)
150 virtual ~BCCSMatrixInitializer()
153 template<
typename Iter>
154 void addRowNnz(
const Iter& row)
const
156 mat->Nnz_+=row->getsize();
159 template<
typename Iter,
typename FullMatrixIndex>
160 void addRowNnz(
const Iter& row,
const std::set<FullMatrixIndex>& indices)
const
162 auto siter =indices.begin();
163 for (
auto entry=row->begin(); entry!=row->end(); ++entry)
165 for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
166 if(siter==indices.end())
168 if(*siter==entry.index())
174 template<
typename Iter,
typename SubMatrixIndex>
175 void addRowNnz(
const Iter& row,
const std::vector<SubMatrixIndex>& indices)
const
177 for (
auto entry=row->begin(); entry!=row->end(); ++entry)
184 allocateMatrixStorage();
188 template<
typename Iter,
typename CIter>
189 void countEntries([[maybe_unused]]
const Iter& row,
const CIter& col)
const
191 countEntries(col.index());
194 void countEntries(size_type colindex)
const
196 for(size_type i=0; i < m; ++i)
198 assert(colindex*m+i<cols);
199 marker[colindex*m+i]+=n;
203 void calcColstart()
const
206 for(size_type i=0; i < cols; ++i) {
208 mat->colstart[i+1]=mat->colstart[i]+marker[i];
209 marker[i]=mat->colstart[i];
213 template<
typename Iter,
typename CIter>
214 void copyValue(
const Iter& row,
const CIter& col)
const
216 copyValue(col, row.index(), col.index());
219 template<
typename CIter>
220 void copyValue(
const CIter& col, size_type rowindex, size_type colindex)
const
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];
233 virtual void createMatrix()
const
240 void allocateMatrixStorage()
const
244 mat->values=
new typename M::field_type[mat->Nnz_];
245 mat->rowindex=
new I[mat->Nnz_];
246 mat->colstart=
new I[cols+1];
249 void allocateMarker()
252 std::fill(marker.begin(), marker.end(), 0);
262 mutable std::vector<size_type> marker;
265 template<
class F,
class Matrix>
266 void copyToBCCSMatrix(F& initializer,
const Matrix& matrix)
268 for (
auto row=matrix.begin(); row!= matrix.end(); ++row)
269 initializer.addRowNnz(row);
271 initializer.allocate();
273 for (
auto row=matrix.begin(); row!= matrix.end(); ++row) {
275 for (
auto col=row->begin(); col != row->end(); ++col)
276 initializer.countEntries(row, col);
279 initializer.calcColstart();
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);
287 initializer.createMatrix();
290 template<
class F,
class M,
class S>
291 void copyToBCCSMatrix(F& initializer,
const MatrixRowSubset<M,S>& mrs)
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;
300 typedef typename MRS::Matrix::size_type size_type;
306 std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
309 for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
310 subMatrixIndex[*index]=s++;
313 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
314 initializer.addRowNnz(row, subMatrixIndex);
316 initializer.allocate();
318 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
319 for(CIter col=row->begin(); col != row->end(); ++col) {
322 initializer.countEntries(subMatrixIndex[col.index()]);
325 initializer.calcColstart();
327 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
328 for(CIter col=row->begin(); col != row->end(); ++col) {
331 initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
333 initializer.createMatrix();
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