3 #ifndef DUNE_MATRIX_UTILS_HH
4 #define DUNE_MATRIX_UTILS_HH
9 #include <dune/common/typetraits.hh>
10 #include <dune/common/static_assert.hh>
11 #include <dune/common/fmatrix.hh>
12 #include <dune/common/diagonalmatrix.hh>
13 #include <dune/common/unused.hh>
21 template<
typename B,
typename A>
24 template<
typename K,
int n,
int m>
27 template<
class T,
class A>
47 static typename M::size_type
count(
const M& matrix)
49 typedef typename M::ConstRowIterator RowIterator;
51 RowIterator endRow = matrix.end();
52 typename M::size_type nonZeros = 0;
54 for(RowIterator
row = matrix.begin();
row != endRow; ++
row) {
55 typedef typename M::ConstColIterator Entry;
56 Entry endEntry =
row->end();
57 for(Entry entry =
row->begin(); entry != endEntry; ++entry) {
66 struct NonZeroCounter<1>
69 static typename M::size_type
count(
const M& matrix)
71 return matrix.N()*matrix.M();
81 template<
class Matrix, std::
size_t blocklevel, std::
size_t l=blocklevel>
90 DUNE_UNUSED_PARAMETER(mat);
91 #ifdef DUNE_ISTL_WITH_CHECKING
95 Entry diagonal =
row->find(
row.index());
96 if(diagonal==
row->end())
97 DUNE_THROW(
ISTLError,
"Missing diagonal value in row "<<
row.index()
98 <<
" at block recursion level "<<l-blocklevel);
106 template<
class Matrix, std::
size_t l>
114 DUNE_THROW(
ISTLError,
"Missing diagonal value in row "<<
row.index()
115 <<
" at block recursion level "<<l);
120 template<
typename T1,
typename T2,
typename T3,
typename T4,
typename T5,
121 typename T6,
typename T7,
typename T8,
typename T9>
122 class MultiTypeBlockMatrix;
124 template<
typename T1,
typename T2,
typename T3,
typename T4,
typename T5,
125 typename T6,
typename T7,
typename T8,
typename T9, std::size_t blocklevel, std::size_t l>
129 typedef MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>
Matrix;
137 #ifdef DUNE_ISTL_WITH_CHECKING
168 template<
class G,
class M>
169 bool operator()(
const std::pair<G,M>& p1,
const std::pair<G,M>& p2)
171 return p1.first<p2.first;
176 template<
class M,
class C>
179 typedef typename C::ParallelIndexSet::const_iterator IIter;
180 typedef typename C::OwnerSet OwnerSet;
181 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
185 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
187 gmax=std::max(gmax,idx->global());
189 gmax=ooc.communicator().max(gmax);
190 ooc.buildGlobalLookup();
192 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
194 if(OwnerSet::contains(idx->local().attribute()))
196 typedef typename M::block_type Block;
198 std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
201 typedef typename M::ConstColIterator CIter;
202 for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
204 const typename C::ParallelIndexSet::IndexPair* pair
205 =ooc.globalLookup().pair(c.index());
207 entries.insert(std::make_pair(pair->global(), *c));
211 GlobalIndex rowidx = idx->global();
212 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
214 cur=ooc.communicator().min(rowidx);
217 typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
218 for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
219 os<<idx->global()<<
" "<<s->first<<
" "<<s->second<<std::endl;
225 ooc.freeGlobalLookup();
227 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
228 while(cur!=ooc.communicator().min(cur)) ;
236 template<
typename B,
typename TA>
258 if (A.j.get()[k]==c) {
294 std::vector<size_type> coldims(A.
M(),
295 std::numeric_limits<size_type>::max());
300 if (coldims[
col.index()]==std::numeric_limits<size_type>::max())
304 for (
typename std::vector<size_type>::iterator it=coldims.begin();
305 it!=coldims.end(); ++it)
315 template<
typename B,
int n,
int m,
typename TA>
340 template<
typename K,
int n,
int m>
367 template<
typename K,
int n,
int m,
typename TA>
394 template<
typename K,
int n>
421 template<
typename K,
int n>
474 template<
typename T,
typename A>