5#ifndef DUNE_ISTL_MATRIXUTILS_HH
6#define DUNE_ISTL_MATRIXUTILS_HH
11#include <dune/common/typetraits.hh>
12#include <dune/common/fmatrix.hh>
13#include <dune/common/dynmatrix.hh>
14#include <dune/common/diagonalmatrix.hh>
15#include <dune/common/scalarmatrixview.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())
61 DUNE_THROW(
ISTLError,
"Missing diagonal value in row "<<row.index()
62 <<
" at block recursion level "<<l-blocklevel);
64 auto m = Impl::asMatrix(*diagonal);
72 template<
class Matrix, std::
size_t l>
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>
101#ifdef DUNE_ISTL_WITH_CHECKING
120 [[maybe_unused]]
typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae =
nullptr)
127 [[maybe_unused]]
typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae =
nullptr)
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>
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();
189 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
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();
204 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
205 while(cur!=ooc.communicator().min(cur)) ;
212 static_assert(IsNumber<M>::value,
"MatrixDimension is not implemented for this type!");
226 template<
typename B,
typename TA>
260 template<
typename B,
typename TA>
269 const B* row = A.r[i].getptr();
282 if (A.j_.get()[k] == c) {
292 B* a = A.r[i].getptr();
293 for (
size_type k=0; k<A.r[i].getsize(); k++)
318 std::vector<size_type> coldims(A.M(),
319 std::numeric_limits<size_type>::max());
321 for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
322 for (ConstColIterator
col=row->begin();
col!=row->end(); ++
col)
324 if (coldims[
col.index()]==std::numeric_limits<size_type>::max())
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>
364 template<
typename K,
int n,
int m>
418 template<
typename K,
int n,
int m,
typename TA>
445 template<
typename K,
int n>
472 template<
typename K,
int n>
525 template<
typename T,
typename A>
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Col col
Definition: matrixmatrix.hh:351
Matrix & mat
Definition: matrixmatrix.hh:347
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
Definition: allocator.hh:11
void printGlobalSparseMatrix(const M &mat, C &ooc, std::ostream &os)
Definition: matrixutils.hh:154
Definition: matrixutils.hh:211
static auto coldim(const M &A)
Definition: matrixutils.hh:219
static auto rowdim(const M &A)
Definition: matrixutils.hh:214
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:500
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:741
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:491
Iterator access to matrix rows
Definition: bcrsmatrix.hh:579
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:46
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
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:620
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:614
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
Definition: matrixutils.hh:27
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
static void check(const Matrix &mat)
Definition: matrixutils.hh:75
static void check(const Matrix &)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:99
MultiTypeBlockMatrix< T1, Args... > Matrix
Definition: matrixutils.hh:93
static size_type rowdim(const Matrix< B, TA > &A, size_type i)
Definition: matrixutils.hh:232
static size_type coldim(const Matrix< B, TA > &A)
Definition: matrixutils.hh:250
static size_type rowdim(const Matrix< B, TA > &A)
Definition: matrixutils.hh:242
typename Matrix< B, TA >::size_type size_type
Definition: matrixutils.hh:230
static size_type coldim(const Matrix< B, TA > &A, size_type c)
Definition: matrixutils.hh:237
typename Matrix< B, TA >::block_type block_type
Definition: matrixutils.hh:229
BCRSMatrix< B, TA > Matrix
Definition: matrixutils.hh:263
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:311
Matrix::block_type block_type
Definition: matrixutils.hh:264
static size_type coldim(const Matrix &A, size_type c)
Definition: matrixutils.hh:276
Matrix::size_type size_type
Definition: matrixutils.hh:265
static size_type rowdim(const Matrix &A, size_type i)
Definition: matrixutils.hh:267
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:304
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:359
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:345
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:355
Matrix::size_type size_type
Definition: matrixutils.hh:343
BCRSMatrix< FieldMatrix< B, n, m >,TA > Matrix
Definition: matrixutils.hh:342
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:350
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:370
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:385
Matrix::size_type size_type
Definition: matrixutils.hh:368
FieldMatrix< K, n, m > Matrix
Definition: matrixutils.hh:367
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:375
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:380
static size_type coldim(const MatrixType &A)
Definition: matrixutils.hh:412
static size_type rowdim(const MatrixType &A)
Definition: matrixutils.hh:407
static size_type rowdim(const MatrixType &, size_type)
Definition: matrixutils.hh:397
MatrixType::size_type size_type
Definition: matrixutils.hh:395
static size_type coldim(const MatrixType &, size_type)
Definition: matrixutils.hh:402
Dune::DynamicMatrix< T > MatrixType
Definition: matrixutils.hh:394
static size_type coldim(const ThisMatrix &A)
Definition: matrixutils.hh:439
static size_type rowdim(const ThisMatrix &A)
Definition: matrixutils.hh:434
Matrix< FieldMatrix< K, n, m >, TA > ThisMatrix
Definition: matrixutils.hh:421
static size_type coldim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:429
static size_type rowdim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:424
ThisMatrix::size_type size_type
Definition: matrixutils.hh:422
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:456
Matrix::size_type size_type
Definition: matrixutils.hh:449
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:466
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:461
DiagonalMatrix< K, n > Matrix
Definition: matrixutils.hh:448
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:451
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:493
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:478
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:483
Matrix::size_type size_type
Definition: matrixutils.hh:476
ScaledIdentityMatrix< K, n > Matrix
Definition: matrixutils.hh:475
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:488
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:504
@ value
True if T is an ISTL matrix.
Definition: matrixutils.hh:509
Definition: matrixutils.hh:538
bool operator()(const T *l, const T *r)
Definition: matrixutils.hh:539
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:30
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:43