5#ifndef DUNE_ISTL_MATRIXUTILS_HH
6#define DUNE_ISTL_MATRIXUTILS_HH
17#include "istlexception.hh"
23 template<
typename B,
typename A>
26 template<
typename K,
int n,
int m>
29 template<
class T,
class A>
46 template<
class Matrix, std::
size_t blocklevel, std::
size_t l=blocklevel>
55#ifdef DUNE_ISTL_WITH_CHECKING
58 for(Row row = mat.begin(); row!=mat.end(); ++row) {
59 Entry diagonal = row->find(row.index());
60 if(diagonal==row->end())
62 <<
" at block recursion level "<<l-blocklevel);
64 auto m = Impl::asMatrix(*diagonal);
72 template<
class Matrix, std::
size_t l>
73 struct CheckIfDiagonalPresent<Matrix,0,l>
75 static void check(
const Matrix& mat)
78 for(Row row = mat.begin(); row!=mat.end(); ++row) {
79 if(row->find(row.index())==row->end())
80 DUNE_THROW(ISTLError,
"Missing diagonal value in row "<<row.index()
81 <<
" at block recursion level "<<l);
86 template<
typename FirstRow,
typename... Args>
87 class MultiTypeBlockMatrix;
89 template<std::size_t blocklevel, std::size_t l,
typename T1,
typename... Args>
90 struct CheckIfDiagonalPresent<MultiTypeBlockMatrix<T1,Args...>,
93 typedef MultiTypeBlockMatrix<T1,Args...> Matrix;
99 static void check(
const Matrix& )
101#ifdef DUNE_ISTL_WITH_CHECKING
129 typename M::size_type nonZeros = 0;
130 for(
auto&& row : matrix)
131 for(
auto&& entry : row)
145 template<
class G,
class M>
146 bool operator()(
const std::pair<G,M>& p1,
const std::pair<G,M>& p2)
const
148 return p1.first<p2.first;
153 template<
class M,
class C>
154 void printGlobalSparseMatrix(
const M& mat, C& ooc, std::ostream& os)
156 typedef typename C::ParallelIndexSet::const_iterator IIter;
157 typedef typename C::OwnerSet OwnerSet;
158 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
162 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
164 gmax=std::max(gmax,idx->global());
166 gmax=ooc.communicator().max(gmax);
167 ooc.buildGlobalLookup();
169 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
171 if(OwnerSet::contains(idx->local().attribute()))
173 typedef typename M::block_type Block;
175 std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
178 typedef typename M::ConstColIterator CIter;
179 for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
181 const typename C::ParallelIndexSet::IndexPair* pair
182 =ooc.globalLookup().pair(c.index());
184 entries.insert(std::make_pair(pair->global(), *c));
188 GlobalIndex rowidx = idx->global();
191 cur=ooc.communicator().min(rowidx);
194 typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
195 for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
196 os<<idx->global()<<
" "<<s->first<<
" "<<s->second<<std::endl;
202 ooc.freeGlobalLookup();
205 while(cur!=ooc.communicator().min(cur)) ;
210 struct MatrixDimension
212 static_assert(IsNumber<M>::value,
"MatrixDimension is not implemented for this type!");
214 static auto rowdim(
const M& A)
219 static auto coldim(
const M& A)
226 template<
typename B,
typename TA>
227 struct MatrixDimension<Matrix<B,TA> >
232 static size_type rowdim (
const Matrix<B,TA>& A, size_type i)
234 return MatrixDimension<block_type>::rowdim(A[i][0]);
237 static size_type coldim (
const Matrix<B,TA>& A, size_type c)
239 return MatrixDimension<block_type>::coldim(A[0][c]);
242 static size_type rowdim (
const Matrix<B,TA>& A)
245 for (size_type i=0; i<A.N(); i++)
250 static size_type coldim (
const Matrix<B,TA>& A)
253 for (size_type i=0; i<A.M(); i++)
260 template<
typename B,
typename TA>
261 struct MatrixDimension<BCRSMatrix<B,TA> >
263 typedef BCRSMatrix<B,TA> Matrix;
267 static size_type rowdim (
const Matrix& A, size_type i)
269 const B* row = A.r[i].getptr();
271 return MatrixDimension<block_type>::rowdim(*row);
276 static size_type coldim (
const Matrix& A, size_type c)
281 for (size_type k=0; k<A.nnz_; k++) {
282 if (A.j_.get()[k] == c) {
283 return MatrixDimension<block_type>::coldim(A.a[k]);
289 for (size_type i=0; i<A.N(); i++)
291 size_type* j = A.r[i].getindexptr();
292 B* a = A.r[i].getptr();
293 for (size_type k=0; k<A.r[i].getsize(); k++)
295 return MatrixDimension<block_type>::coldim(a[k]);
304 static size_type rowdim (
const Matrix& A){
306 for (size_type i=0; i<A.N(); i++)
311 static size_type coldim (
const Matrix& A){
318 std::vector<size_type> coldims(A.M(),
321 for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
322 for (ConstColIterator col=row->begin(); col!=row->end(); ++col)
325 coldims[col.index()] = MatrixDimension<block_type>::coldim(*col);
328 for (
typename std::vector<size_type>::iterator it=coldims.begin();
329 it!=coldims.end(); ++it)
339 template<
typename B,
int n,
int m,
typename TA>
340 struct MatrixDimension<BCRSMatrix<FieldMatrix<B,n,m> ,TA> >
342 typedef BCRSMatrix<FieldMatrix<B,n,m> ,TA> Matrix;
345 static size_type rowdim (
const Matrix& , size_type )
350 static size_type coldim (
const Matrix& , size_type )
355 static size_type rowdim (
const Matrix& A) {
359 static size_type coldim (
const Matrix& A) {
364 template<
typename K,
int n,
int m>
365 struct MatrixDimension<FieldMatrix<K,n,m> >
367 typedef FieldMatrix<K,n,m> Matrix;
370 static size_type rowdim(
const Matrix& , size_type )
375 static size_type coldim(
const Matrix& , size_type )
380 static size_type rowdim(
const Matrix& )
385 static size_type coldim(
const Matrix& )
392 struct MatrixDimension<
Dune::DynamicMatrix<T> >
395 typedef typename MatrixType::size_type size_type;
397 static size_type rowdim(
const MatrixType& , size_type )
402 static size_type coldim(
const MatrixType& , size_type )
407 static size_type rowdim(
const MatrixType& A)
412 static size_type coldim(
const MatrixType& A)
418 template<
typename K,
int n,
int m,
typename TA>
419 struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
421 typedef Matrix<FieldMatrix<K,n,m>, TA> ThisMatrix;
422 typedef typename ThisMatrix::size_type size_type;
424 static size_type rowdim(
const ThisMatrix& , size_type )
429 static size_type coldim(
const ThisMatrix& , size_type )
434 static size_type rowdim(
const ThisMatrix& A)
439 static size_type coldim(
const ThisMatrix& A)
445 template<
typename K,
int n>
446 struct MatrixDimension<DiagonalMatrix<K,n> >
448 typedef DiagonalMatrix<K,n> Matrix;
451 static size_type rowdim(
const Matrix& , size_type )
456 static size_type coldim(
const Matrix& , size_type )
461 static size_type rowdim(
const Matrix& )
466 static size_type coldim(
const Matrix& )
472 template<
typename K,
int n>
473 struct MatrixDimension<ScaledIdentityMatrix<K,n> >
475 typedef ScaledIdentityMatrix<K,n> Matrix;
478 static size_type rowdim(
const Matrix& , size_type )
483 static size_type coldim(
const Matrix& , size_type )
488 static size_type rowdim(
const Matrix& )
493 static size_type coldim(
const Matrix& )
514 struct IsMatrix<DenseMatrix<T> >
525 template<
typename T,
typename A>
526 struct IsMatrix<BCRSMatrix<T,A> >
537 struct PointerCompare
539 bool operator()(
const T* l,
const T* r)
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
derive error class from the base class in common
Definition: istlexception.hh:19
ConstIterator class for sequential access.
Definition: matrix.hh:404
A generic dynamic dense matrix.
Definition: matrix.hh:561
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:577
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:586
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
T block_type
Export the type representing the components.
Definition: matrix.hh:568
This file implements a quadratic diagonal matrix of fixed size.
This file implements a dense matrix with dynamic numbers of rows and columns.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
auto countNonZeros(const M &, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:119
Dune namespace.
Definition: alignedallocator.hh:13
Implements a scalar matrix view wrapper around an existing scalar.
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:48
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:53
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:504
@ value
True if T is an ISTL matrix.
Definition: matrixutils.hh:509
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Traits for type conversions and type information.