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>
19 template<
typename B,
typename A>
22 template<
typename K,
int n,
int m>
25 template<
class T,
class A>
45 static typename M::size_type
count(
const M& matrix)
47 typedef typename M::ConstRowIterator RowIterator;
49 RowIterator endRow = matrix.end();
50 typename M::size_type nonZeros = 0;
52 for(RowIterator
row = matrix.begin();
row != endRow; ++
row){
53 typedef typename M::ConstColIterator Entry;
54 Entry endEntry =
row->end();
55 for(Entry entry =
row->begin(); entry != endEntry; ++entry){
64 struct NonZeroCounter<1>
67 static typename M::size_type
count(
const M& matrix)
69 return matrix.N()*matrix.M();
79 template<
class Matrix, std::
size_t blocklevel, std::
size_t l=blocklevel>
88 #ifdef DUNE_ISTL_WITH_CHECKING
92 Entry diagonal =
row->find(
row.index());
93 if(diagonal==
row->end())
94 DUNE_THROW(
ISTLError,
"Missing diagonal value in row "<<
row.index()
95 <<
" at block recursion level "<<l-blocklevel);
103 template<
class Matrix, std::
size_t l>
111 DUNE_THROW(
ISTLError,
"Missing diagonal value in row "<<
row.index()
112 <<
" at block recursion level "<<l);
117 template<
typename T1,
typename T2,
typename T3,
typename T4,
typename T5,
118 typename T6,
typename T7,
typename T8,
typename T9>
119 class MultiTypeBlockMatrix;
121 template<
typename T1,
typename T2,
typename T3,
typename T4,
typename T5,
122 typename T6,
typename T7,
typename T8,
typename T9, std::size_t blocklevel, std::size_t l>
126 typedef MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>
Matrix;
134 #ifdef DUNE_ISTL_WITH_CHECKING
165 template<
class G,
class M>
166 bool operator()(
const std::pair<G,M>& p1,
const std::pair<G,M>& p2)
168 return p1.first<p2.first;
173 template<
class M,
class C>
176 typedef typename C::ParallelIndexSet::const_iterator IIter;
177 typedef typename C::OwnerSet OwnerSet;
178 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
182 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
184 gmax=std::max(gmax,idx->global());
186 gmax=ooc.communicator().max(gmax);
187 ooc.buildGlobalLookup();
189 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end();
191 if(OwnerSet::contains(idx->local().attribute()))
193 typedef typename M::block_type Block;
195 std::set<std::pair<GlobalIndex,Block>,CompPair> entries;
198 typedef typename M::ConstColIterator CIter;
199 for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end();
201 const typename C::ParallelIndexSet::IndexPair* pair
202 =ooc.globalLookup().pair(c.index());
204 entries.insert(std::make_pair(pair->global(), *c));
208 GlobalIndex rowidx = idx->global();
209 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
211 cur=ooc.communicator().min(rowidx);
214 typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter;
215 for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s)
216 os<<idx->global()<<
" "<<s->first<<
" "<<s->second<<std::endl;
222 ooc.freeGlobalLookup();
224 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max();
225 while(cur!=ooc.communicator().min(cur));
234 template<
typename B,
typename TA>
256 if (A.j.get()[k]==c) {
292 std::vector<size_type> coldims(A.
M(),
293 std::numeric_limits<size_type>::max());
298 if (coldims[
col.index()]==std::numeric_limits<size_type>::max())
302 for (
typename std::vector<size_type>::iterator it=coldims.begin();
303 it!=coldims.end(); ++it)
313 template<
typename B,
int n,
int m,
typename TA>
338 template<
typename K,
int n,
int m>
365 template<
typename K,
int n,
int m,
typename TA>
418 template<
typename T,
typename A>