5#ifndef DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
6#define DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
14#include <dune/istl/bccsmatrix.hh>
18 template<
class I,
class S,
class D>
19 class OverlappingSchwarzInitializer;
22namespace Dune::ISTL::Impl
31 template<
class M,
class S>
38 typedef S RowIndexSet;
45 MatrixRowSubset(
const Matrix& m,
const RowIndexSet& s)
49 const Matrix& matrix()
const
54 const RowIndexSet& rowIndexSet()
const
61 :
public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
64 const_iterator(
typename Matrix::const_iterator firstRow,
65 typename RowIndexSet::const_iterator pos)
66 : firstRow_(firstRow), pos_(pos)
72 return *(firstRow_+ *pos_);
74 bool equals(
const const_iterator& o)
const
82 typename RowIndexSet::value_type index()
const
89 typename Matrix::const_iterator firstRow_;
91 typename RowIndexSet::const_iterator pos_;
95 const_iterator begin()
const
97 return const_iterator(m_.begin(), s_.begin());
100 const_iterator end()
const
102 return const_iterator(m_.begin(), s_.end());
109 const RowIndexSet& s_;
118 template<
class M,
class I =
typename M::
size_type>
119 class BCCSMatrixInitializer
121 template<
class IList,
class S,
class D>
126 typedef Dune::ISTL::Impl::BCCSMatrix<typename Matrix::field_type, I> OutputMatrix;
131 BCCSMatrixInitializer(OutputMatrix& mat_)
132 : mat(&mat_), cols(mat_.M())
141 n = M::block_type::rows;
142 m = M::block_type::cols;
148 BCCSMatrixInitializer()
149 : mat(0), cols(0), n(0), m(0)
152 virtual ~BCCSMatrixInitializer()
155 template<
typename Iter>
156 void addRowNnz(
const Iter& row)
const
158 mat->Nnz_+=row->getsize();
161 template<
typename Iter,
typename FullMatrixIndex>
162 void addRowNnz(
const Iter& row,
const std::set<FullMatrixIndex>& indices)
const
164 auto siter =indices.begin();
165 for (
auto entry=row->begin(); entry!=row->end(); ++entry)
167 for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
168 if(siter==indices.end())
170 if(*siter==entry.index())
176 template<
typename Iter,
typename SubMatrixIndex>
177 void addRowNnz(
const Iter& row,
const std::vector<SubMatrixIndex>& indices)
const
179 for (
auto entry=row->begin(); entry!=row->end(); ++entry)
186 allocateMatrixStorage();
190 template<
typename Iter,
typename CIter>
191 void countEntries([[maybe_unused]]
const Iter& row,
const CIter& col)
const
193 countEntries(col.index());
196 void countEntries(size_type colindex)
const
198 for(size_type i=0; i < m; ++i)
200 assert(colindex*m+i<cols);
201 marker[colindex*m+i]+=n;
205 void calcColstart()
const
208 for(size_type i=0; i < cols; ++i) {
210 mat->colstart[i+1]=mat->colstart[i]+marker[i];
211 marker[i]=mat->colstart[i];
215 template<
typename Iter,
typename CIter>
216 void copyValue(
const Iter& row,
const CIter& col)
const
218 copyValue(col, row.index(), col.index());
221 template<
typename CIter>
222 void copyValue(
const CIter& col, size_type rowindex, size_type colindex)
const
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];
235 virtual void createMatrix()
const
242 void allocateMatrixStorage()
const
246 mat->values=
new typename M::field_type[mat->Nnz_];
247 mat->rowindex=
new I[mat->Nnz_];
248 mat->colstart=
new I[cols+1];
251 void allocateMarker()
254 std::fill(marker.begin(), marker.end(), 0);
264 mutable std::vector<size_type> marker;
267 template<
class F,
class Matrix>
268 void copyToBCCSMatrix(F& initializer,
const Matrix& matrix)
270 for (
auto row=matrix.begin(); row!= matrix.end(); ++row)
271 initializer.addRowNnz(row);
273 initializer.allocate();
275 for (
auto row=matrix.begin(); row!= matrix.end(); ++row) {
277 for (
auto col=row->begin(); col != row->end(); ++col)
278 initializer.countEntries(row, col);
281 initializer.calcColstart();
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);
289 initializer.createMatrix();
292 template<
class F,
class M,
class S>
293 void copyToBCCSMatrix(F& initializer,
const MatrixRowSubset<M,S>& mrs)
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;
302 typedef typename MRS::Matrix::size_type size_type;
308 std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
311 for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
312 subMatrixIndex[*index]=s++;
315 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
316 initializer.addRowNnz(row, subMatrixIndex);
318 initializer.allocate();
320 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
321 for(CIter col=row->begin(); col != row->end(); ++col) {
324 initializer.countEntries(subMatrixIndex[col.index()]);
327 initializer.calcColstart();
329 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
330 for(CIter col=row->begin(); col != row->end(); ++col) {
333 initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
335 initializer.createMatrix();
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