3 #ifndef DUNE_ISTL_MATRIXUTILS_HH
4 #define DUNE_ISTL_MATRIXUTILS_HH
9 #include <dune/common/typetraits.hh>
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/dynmatrix.hh>
12 #include <dune/common/diagonalmatrix.hh>
13 #include <dune/common/scalarmatrixview.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())
59 DUNE_THROW(
ISTLError,
"Missing diagonal value in row "<<row.index()
60 <<
" at block recursion level "<<l-blocklevel);
62 auto m = Impl::asMatrix(*diagonal);
70 template<
class Matrix, std::
size_t l>
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>
99 #ifdef DUNE_ISTL_WITH_CHECKING
118 [[maybe_unused]]
typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae =
nullptr)
125 [[maybe_unused]]
typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae =
nullptr)
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>
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();
162 gmax=std::max(gmax,idx->global());
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();
187 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
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();
202 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
203 while(cur!=ooc.communicator().min(cur)) ;
210 static_assert(IsNumber<M>::value,
"MatrixDimension is not implemented for this type!");
224 template<
typename B,
typename TA>
258 template<
typename B,
typename TA>
267 const B* row = A.r[i].getptr();
280 if (A.j_.get()[k] == c) {
290 B* a = A.r[i].getptr();
291 for (
size_type k=0; k<A.r[i].getsize(); k++)
316 std::vector<size_type> coldims(A.M(),
317 std::numeric_limits<size_type>::max());
319 for (ConstRowIterator row=A.begin(); row!=A.end(); ++row)
320 for (ConstColIterator
col=row->begin();
col!=row->end(); ++
col)
322 if (coldims[
col.index()]==std::numeric_limits<size_type>::max())
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>
362 template<
typename K,
int n,
int m>
416 template<
typename K,
int n,
int m,
typename TA>
443 template<
typename K,
int n>
470 template<
typename K,
int n>
523 template<
typename T,
typename A>
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
auto countNonZeros(const M &, [[maybe_unused]] typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:117
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:9
void printGlobalSparseMatrix(const M &mat, C &ooc, std::ostream &os)
Definition: matrixutils.hh:152
Definition: matrixutils.hh:209
static auto coldim(const M &A)
Definition: matrixutils.hh:217
static auto rowdim(const M &A)
Definition: matrixutils.hh:212
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:498
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:739
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:489
Iterator access to matrix rows
Definition: bcrsmatrix.hh:577
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
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
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:618
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:612
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
Definition: matrixutils.hh:25
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:46
static void check([[maybe_unused]] const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:51
static void check(const Matrix &mat)
Definition: matrixutils.hh:73
static void check(const Matrix &)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:97
MultiTypeBlockMatrix< T1, Args... > Matrix
Definition: matrixutils.hh:91
static size_type rowdim(const Matrix< B, TA > &A, size_type i)
Definition: matrixutils.hh:230
static size_type coldim(const Matrix< B, TA > &A)
Definition: matrixutils.hh:248
static size_type rowdim(const Matrix< B, TA > &A)
Definition: matrixutils.hh:240
typename Matrix< B, TA >::size_type size_type
Definition: matrixutils.hh:228
static size_type coldim(const Matrix< B, TA > &A, size_type c)
Definition: matrixutils.hh:235
typename Matrix< B, TA >::block_type block_type
Definition: matrixutils.hh:227
BCRSMatrix< B, TA > Matrix
Definition: matrixutils.hh:261
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:309
Matrix::block_type block_type
Definition: matrixutils.hh:262
static size_type coldim(const Matrix &A, size_type c)
Definition: matrixutils.hh:274
Matrix::size_type size_type
Definition: matrixutils.hh:263
static size_type rowdim(const Matrix &A, size_type i)
Definition: matrixutils.hh:265
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:302
static size_type coldim(const Matrix &A)
Definition: matrixutils.hh:357
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:343
static size_type rowdim(const Matrix &A)
Definition: matrixutils.hh:353
Matrix::size_type size_type
Definition: matrixutils.hh:341
BCRSMatrix< FieldMatrix< B, n, m >,TA > Matrix
Definition: matrixutils.hh:340
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:348
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:368
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:383
Matrix::size_type size_type
Definition: matrixutils.hh:366
FieldMatrix< K, n, m > Matrix
Definition: matrixutils.hh:365
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:373
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:378
static size_type coldim(const MatrixType &A)
Definition: matrixutils.hh:410
static size_type rowdim(const MatrixType &A)
Definition: matrixutils.hh:405
static size_type rowdim(const MatrixType &, size_type)
Definition: matrixutils.hh:395
MatrixType::size_type size_type
Definition: matrixutils.hh:393
static size_type coldim(const MatrixType &, size_type)
Definition: matrixutils.hh:400
Dune::DynamicMatrix< T > MatrixType
Definition: matrixutils.hh:392
static size_type coldim(const ThisMatrix &A)
Definition: matrixutils.hh:437
static size_type rowdim(const ThisMatrix &A)
Definition: matrixutils.hh:432
Matrix< FieldMatrix< K, n, m >, TA > ThisMatrix
Definition: matrixutils.hh:419
static size_type coldim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:427
static size_type rowdim(const ThisMatrix &, size_type)
Definition: matrixutils.hh:422
ThisMatrix::size_type size_type
Definition: matrixutils.hh:420
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:454
Matrix::size_type size_type
Definition: matrixutils.hh:447
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:464
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:459
DiagonalMatrix< K, n > Matrix
Definition: matrixutils.hh:446
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:449
static size_type coldim(const Matrix &)
Definition: matrixutils.hh:491
static size_type rowdim(const Matrix &, size_type)
Definition: matrixutils.hh:476
static size_type coldim(const Matrix &, size_type)
Definition: matrixutils.hh:481
Matrix::size_type size_type
Definition: matrixutils.hh:474
ScaledIdentityMatrix< K, n > Matrix
Definition: matrixutils.hh:473
static size_type rowdim(const Matrix &)
Definition: matrixutils.hh:486
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
@ value
True if T is an ISTL matrix.
Definition: matrixutils.hh:507
Definition: matrixutils.hh:536
bool operator()(const T *l, const T *r)
Definition: matrixutils.hh:537
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:28
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41