3#ifndef DUNE_ISTL_MATRIXUTILS_HH
4#define DUNE_ISTL_MATRIXUTILS_HH
15#include "istlexception.hh"
21 template<
typename B,
typename A>
24 template<
typename K,
int n,
int m>
27 template<
class T,
class A>
44 template<
class Matrix, std::
size_t blocklevel, std::
size_t l=blocklevel>
53#ifdef DUNE_ISTL_WITH_CHECKING
56 for(Row row = mat.begin(); row!=mat.end(); ++row) {
57 Entry diagonal = row->find(row.index());
58 if(diagonal==row->end())
60 <<
" at block recursion level "<<l-blocklevel);
62 auto m = Impl::asMatrix(*diagonal);
70 template<
class Matrix, std::
size_t l>
71 struct CheckIfDiagonalPresent<Matrix,0,l>
73 static void check(
const Matrix& mat)
76 for(Row row = mat.begin(); row!=mat.end(); ++row) {
77 if(row->find(row.index())==row->end())
78 DUNE_THROW(ISTLError,
"Missing diagonal value in row "<<row.index()
79 <<
" at block recursion level "<<l);
84 template<
typename FirstRow,
typename... Args>
85 class MultiTypeBlockMatrix;
87 template<std::size_t blocklevel, std::size_t l,
typename T1,
typename... Args>
88 struct CheckIfDiagonalPresent<MultiTypeBlockMatrix<T1,Args...>,
91 typedef MultiTypeBlockMatrix<T1,Args...> Matrix;
97 static void check(
const Matrix& )
99#ifdef DUNE_ISTL_WITH_CHECKING
127 typename M::size_type nonZeros = 0;
128 for(
auto&& row : matrix)
129 for(
auto&& entry : row)
143 template<
class G,
class M>
144 bool operator()(
const std::pair<G,M>& p1,
const std::pair<G,M>& p2)
const
146 return p1.first<p2.first;
151 template<
class M,
class C>
152 void printGlobalSparseMatrix(
const M& mat, C& ooc, std::ostream& os)
154 typedef typename C::ParallelIndexSet::const_iterator IIter;
155 typedef typename C::OwnerSet OwnerSet;
156 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
160 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
164 gmax=ooc.communicator().max(gmax);
165 ooc.buildGlobalLookup();
167 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
169 if(OwnerSet::contains(idx->local().attribute()))
171 typedef typename M::block_type Block;
173 std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
176 typedef typename M::ConstColIterator CIter;
177 for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
179 const typename C::ParallelIndexSet::IndexPair* pair
180 =ooc.globalLookup().pair(c.index());
182 entries.insert(std::make_pair(pair->global(), *c));
186 GlobalIndex rowidx = idx->global();
189 cur=ooc.communicator().min(rowidx);
192 typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
193 for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
194 os<<idx->global()<<
" "<<s->first<<
" "<<s->second<<std::endl;
200 ooc.freeGlobalLookup();
203 while(cur!=ooc.communicator().min(cur)) ;
208 struct MatrixDimension
210 static_assert(IsNumber<M>::value,
"MatrixDimension is not implemented for this type!");
212 static auto rowdim(
const M& A)
217 static auto coldim(
const M& A)
224 template<
typename B,
typename TA>
225 struct MatrixDimension<Matrix<B,TA> >
230 static size_type rowdim (
const Matrix<B,TA>& A, size_type i)
232 return MatrixDimension<block_type>::rowdim(A[i][0]);
235 static size_type coldim (
const Matrix<B,TA>& A, size_type c)
237 return MatrixDimension<block_type>::coldim(A[0][c]);
240 static size_type rowdim (
const Matrix<B,TA>& A)
243 for (size_type i=0; i<A.N(); i++)
248 static size_type coldim (
const Matrix<B,TA>& A)
251 for (size_type i=0; i<A.M(); i++)
258 template<
typename B,
typename TA>
259 struct MatrixDimension<BCRSMatrix<B,TA> >
261 typedef BCRSMatrix<B,TA> Matrix;
265 static size_type rowdim (
const Matrix& A, size_type i)
267 const B* row = A.r[i].getptr();
269 return MatrixDimension<block_type>::rowdim(*row);
274 static size_type coldim (
const Matrix& A, size_type c)
279 for (size_type k=0; k<A.nnz_; k++) {
280 if (A.j_.get()[k] == c) {
281 return MatrixDimension<block_type>::coldim(A.a[k]);
287 for (size_type i=0; i<A.N(); i++)
289 size_type* j = A.r[i].getindexptr();
290 B* a = A.r[i].getptr();
291 for (size_type k=0; k<A.r[i].getsize(); k++)
293 return MatrixDimension<block_type>::coldim(a[k]);
302 static size_type rowdim (
const Matrix& A){
304 for (size_type i=0; i<A.N(); i++)
309 static size_type coldim (
const Matrix& A){
316 std::vector<size_type> coldims(A.M(),
319 for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
320 for (ConstColIterator col=row->begin(); col!=row->end(); ++col)
323 coldims[col.index()] = MatrixDimension<block_type>::coldim(*col);
326 for (
typename std::vector<size_type>::iterator it=coldims.begin();
327 it!=coldims.end(); ++it)
337 template<
typename B,
int n,
int m,
typename TA>
338 struct MatrixDimension<BCRSMatrix<FieldMatrix<B,n,m> ,TA> >
340 typedef BCRSMatrix<FieldMatrix<B,n,m> ,TA> Matrix;
343 static size_type rowdim (
const Matrix& , size_type )
348 static size_type coldim (
const Matrix& , size_type )
353 static size_type rowdim (
const Matrix& A) {
357 static size_type coldim (
const Matrix& A) {
362 template<
typename K,
int n,
int m>
363 struct MatrixDimension<FieldMatrix<K,n,m> >
365 typedef FieldMatrix<K,n,m> Matrix;
368 static size_type rowdim(
const Matrix& , size_type )
373 static size_type coldim(
const Matrix& , size_type )
378 static size_type rowdim(
const Matrix& )
383 static size_type coldim(
const Matrix& )
390 struct MatrixDimension<
Dune::DynamicMatrix<T> >
393 typedef typename MatrixType::size_type size_type;
395 static size_type rowdim(
const MatrixType& , size_type )
400 static size_type coldim(
const MatrixType& , size_type )
405 static size_type rowdim(
const MatrixType& A)
410 static size_type coldim(
const MatrixType& A)
416 template<
typename K,
int n,
int m,
typename TA>
417 struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> >
419 typedef Matrix<FieldMatrix<K,n,m>, TA> ThisMatrix;
420 typedef typename ThisMatrix::size_type size_type;
422 static size_type rowdim(
const ThisMatrix& , size_type )
427 static size_type coldim(
const ThisMatrix& , size_type )
432 static size_type rowdim(
const ThisMatrix& A)
437 static size_type coldim(
const ThisMatrix& A)
443 template<
typename K,
int n>
444 struct MatrixDimension<DiagonalMatrix<K,n> >
446 typedef DiagonalMatrix<K,n> Matrix;
449 static size_type rowdim(
const Matrix& , size_type )
454 static size_type coldim(
const Matrix& , size_type )
459 static size_type rowdim(
const Matrix& )
464 static size_type coldim(
const Matrix& )
470 template<
typename K,
int n>
471 struct MatrixDimension<ScaledIdentityMatrix<K,n> >
473 typedef ScaledIdentityMatrix<K,n> Matrix;
476 static size_type rowdim(
const Matrix& , size_type )
481 static size_type coldim(
const Matrix& , size_type )
486 static size_type rowdim(
const Matrix& )
491 static size_type coldim(
const Matrix& )
512 struct IsMatrix<DenseMatrix<T> >
523 template<
typename T,
typename A>
524 struct IsMatrix<BCRSMatrix<T,A> >
535 struct PointerCompare
537 bool operator()(
const T* l,
const T* r)
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:59
derive error class from the base class in common
Definition: istlexception.hh:17
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:575
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:584
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
T block_type
Export the type representing the components.
Definition: matrix.hh:566
This file implements a quadratic diagonal matrix of fixed size.
Traits for type conversions and type information.
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:216
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:117
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.
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:46
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:51
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
@ value
True if T is an ISTL matrix.
Definition: matrixutils.hh:507
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194