dune-istl
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=8 sw=2 sts=2: 00003 #ifndef DUNE_MATRIX_UTILS_HH 00004 #define DUNE_MATRIX_UTILS_HH 00005 00006 #include<set> 00007 #include<vector> 00008 #include<limits> 00009 #include<dune/common/typetraits.hh> 00010 #include<dune/common/static_assert.hh> 00011 #include<dune/common/fmatrix.hh> 00012 #include<dune/istl/bcrsmatrix.hh> 00013 #include"istlexception.hh" 00014 00015 namespace Dune 00016 { 00017 00018 #ifndef DOYXGEN 00019 template<typename B, typename A> 00020 class BCRSMatrix; 00021 00022 template<typename K, int n, int m> 00023 class FieldMatrix; 00024 00025 template<class T, class A> 00026 class Matrix; 00027 #endif 00028 00038 namespace 00039 { 00040 00041 template<int i> 00042 struct NonZeroCounter 00043 { 00044 template<class M> 00045 static typename M::size_type count(const M& matrix) 00046 { 00047 typedef typename M::ConstRowIterator RowIterator; 00048 00049 RowIterator endRow = matrix.end(); 00050 typename M::size_type nonZeros = 0; 00051 00052 for(RowIterator row = matrix.begin(); row != endRow; ++row){ 00053 typedef typename M::ConstColIterator Entry; 00054 Entry endEntry = row->end(); 00055 for(Entry entry = row->begin(); entry != endEntry; ++entry){ 00056 nonZeros += NonZeroCounter<i-1>::count(*entry); 00057 } 00058 } 00059 return nonZeros; 00060 } 00061 }; 00062 00063 template<> 00064 struct NonZeroCounter<1> 00065 { 00066 template<class M> 00067 static typename M::size_type count(const M& matrix) 00068 { 00069 return matrix.N()*matrix.M(); 00070 } 00071 }; 00072 00073 } 00074 00079 template<class Matrix, std::size_t blocklevel, std::size_t l=blocklevel> 00080 struct CheckIfDiagonalPresent 00081 { 00086 static void check(const Matrix& mat) 00087 { 00088 #ifdef DUNE_ISTL_WITH_CHECKING 00089 typedef typename Matrix::ConstRowIterator Row; 00090 typedef typename Matrix::ConstColIterator Entry; 00091 for(Row row = mat.begin(); row!=mat.end(); ++row){ 00092 Entry diagonal = row->find(row.index()); 00093 if(diagonal==row->end()) 00094 DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index() 00095 <<" at block recursion level "<<l-blocklevel); 00096 else 00097 CheckIfDiagonalPresent<typename Matrix::block_type,blocklevel-1,l>::check(*diagonal); 00098 } 00099 #endif 00100 } 00101 }; 00102 00103 template<class Matrix, std::size_t l> 00104 struct CheckIfDiagonalPresent<Matrix,0,l> 00105 { 00106 static void check(const Matrix& mat) 00107 { 00108 typedef typename Matrix::ConstRowIterator Row; 00109 for(Row row = mat.begin(); row!=mat.end(); ++row){ 00110 if(row->find(row.index())==row->end()) 00111 DUNE_THROW(ISTLError, "Missing diagonal value in row "<<row.index() 00112 <<" at block recursion level "<<l); 00113 } 00114 } 00115 }; 00116 00117 template<typename T1, typename T2, typename T3, typename T4, typename T5, 00118 typename T6, typename T7, typename T8, typename T9> 00119 class MultiTypeBlockMatrix; 00120 00121 template<typename T1, typename T2, typename T3, typename T4, typename T5, 00122 typename T6, typename T7, typename T8, typename T9, std::size_t blocklevel, std::size_t l> 00123 struct CheckIfDiagonalPresent<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>, 00124 blocklevel,l> 00125 { 00126 typedef MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> Matrix; 00127 00132 static void check(const Matrix& mat) 00133 { 00134 #ifdef DUNE_ISTL_WITH_CHECKING 00135 // TODO Implement check 00136 #endif 00137 } 00138 }; 00139 00151 template<class M> 00152 inline int countNonZeros(const M& matrix) 00153 { 00154 return NonZeroCounter<M::blocklevel>::count(matrix); 00155 } 00156 /* 00157 template<class M> 00158 struct ProcessOnFieldsOfMatrix 00159 */ 00160 00162 namespace 00163 { 00164 struct CompPair{ 00165 template<class G,class M> 00166 bool operator()(const std::pair<G,M>& p1, const std::pair<G,M>& p2) 00167 { 00168 return p1.first<p2.first; 00169 } 00170 }; 00171 00172 } 00173 template<class M, class C> 00174 void printGlobalSparseMatrix(const M& mat, C& ooc, std::ostream& os) 00175 { 00176 typedef typename C::ParallelIndexSet::const_iterator IIter; 00177 typedef typename C::OwnerSet OwnerSet; 00178 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex; 00179 00180 GlobalIndex gmax=0; 00181 00182 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end(); 00183 idx!=eidx; ++idx) 00184 gmax=std::max(gmax,idx->global()); 00185 00186 gmax=ooc.communicator().max(gmax); 00187 ooc.buildGlobalLookup(); 00188 00189 for(IIter idx=ooc.indexSet().begin(), eidx=ooc.indexSet().end(); 00190 idx!=eidx; ++idx){ 00191 if(OwnerSet::contains(idx->local().attribute())) 00192 { 00193 typedef typename M::block_type Block; 00194 00195 std::set<std::pair<GlobalIndex,Block>,CompPair> entries; 00196 00197 // sort rows 00198 typedef typename M::ConstColIterator CIter; 00199 for(CIter c=mat[idx->local()].begin(), cend=mat[idx->local()].end(); 00200 c!=cend; ++c){ 00201 const typename C::ParallelIndexSet::IndexPair* pair 00202 =ooc.globalLookup().pair(c.index()); 00203 assert(pair); 00204 entries.insert(std::make_pair(pair->global(), *c)); 00205 } 00206 00207 //wait until its the rows turn. 00208 GlobalIndex rowidx = idx->global(); 00209 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max(); 00210 while(cur!=rowidx) 00211 cur=ooc.communicator().min(rowidx); 00212 00213 // print rows 00214 typedef typename std::set<std::pair<GlobalIndex,Block>,CompPair>::iterator SIter; 00215 for(SIter s=entries.begin(), send=entries.end(); s!=send; ++s) 00216 os<<idx->global()<<" "<<s->first<<" "<<s->second<<std::endl; 00217 00218 00219 } 00220 } 00221 00222 ooc.freeGlobalLookup(); 00223 // Wait until everybody is finished 00224 GlobalIndex cur=std::numeric_limits<GlobalIndex>::max(); 00225 while(cur!=ooc.communicator().min(cur)); 00226 } 00227 00228 template<typename M> 00229 struct MatrixDimension 00230 { 00231 }; 00232 00233 00234 template<typename B, typename TA> 00235 struct MatrixDimension<BCRSMatrix<B,TA> > 00236 { 00237 typedef BCRSMatrix<B,TA> Matrix; 00238 typedef typename Matrix::block_type block_type; 00239 typedef typename Matrix::size_type size_type; 00240 00241 static size_type rowdim (const Matrix& A, size_type i) 00242 { 00243 const B* row = A.r[i].getptr(); 00244 if(row) 00245 return MatrixDimension<block_type>::rowdim(*row); 00246 else 00247 return 0; 00248 } 00249 00250 static size_type coldim (const Matrix& A, size_type c) 00251 { 00252 // find an entry in column j 00253 if (A.nnz>0) 00254 { 00255 for (size_type k=0; k<A.nnz; k++) { 00256 if (A.j[k]==c) { 00257 return MatrixDimension<block_type>::coldim(A.a[k]); 00258 } 00259 } 00260 } 00261 else 00262 { 00263 for (size_type i=0; i<A.N(); i++) 00264 { 00265 size_type* j = A.r[i].getindexptr(); 00266 B* a = A.r[i].getptr(); 00267 for (size_type k=0; k<A.r[i].getsize(); k++) 00268 if (j[k]==c) { 00269 return MatrixDimension<block_type>::coldim(a[k]); 00270 } 00271 } 00272 } 00273 00274 // not found 00275 return 0; 00276 } 00277 00278 static size_type rowdim (const Matrix& A){ 00279 size_type nn=0; 00280 for (size_type i=0; i<A.N(); i++) 00281 nn += rowdim(A,i); 00282 return nn; 00283 } 00284 00285 static size_type coldim (const Matrix& A){ 00286 typedef typename Matrix::ConstRowIterator ConstRowIterator; 00287 typedef typename Matrix::ConstColIterator ConstColIterator; 00288 00289 // The following code has a complexity of nnz, and 00290 // typically a very small constant. 00291 // 00292 std::vector<size_type> coldims(A.M(), 00293 std::numeric_limits<size_type>::max()); 00294 00295 for (ConstRowIterator row=A.begin(); row!=A.end(); ++row) 00296 for (ConstColIterator col=row->begin(); col!=row->end(); ++col) 00297 // only compute blocksizes we don't already have 00298 if (coldims[col.index()]==std::numeric_limits<size_type>::max()) 00299 coldims[col.index()] = MatrixDimension<block_type>::coldim(*col); 00300 00301 size_type sum = 0; 00302 for (typename std::vector<size_type>::iterator it=coldims.begin(); 00303 it!=coldims.end(); ++it) 00304 // skip rows for which no coldim could be determined 00305 if ((*it)>=0) 00306 sum += *it; 00307 00308 return sum; 00309 } 00310 }; 00311 00312 00313 template<typename B, int n, int m, typename TA> 00314 struct MatrixDimension<BCRSMatrix<FieldMatrix<B,n,m> ,TA> > 00315 { 00316 typedef BCRSMatrix<FieldMatrix<B,n,m> ,TA> Matrix; 00317 typedef typename Matrix::size_type size_type; 00318 00319 static size_type rowdim (const Matrix& A, size_type i) 00320 { 00321 return n; 00322 } 00323 00324 static size_type coldim (const Matrix& A, size_type c) 00325 { 00326 return m; 00327 } 00328 00329 static size_type rowdim (const Matrix& A) { 00330 return A.N()*n; 00331 } 00332 00333 static size_type coldim (const Matrix& A) { 00334 return A.M()*m; 00335 } 00336 }; 00337 00338 template<typename K, int n, int m> 00339 struct MatrixDimension<FieldMatrix<K,n,m> > 00340 { 00341 typedef FieldMatrix<K,n,m> Matrix; 00342 typedef typename Matrix::size_type size_type; 00343 00344 static size_type rowdim(const Matrix& A, size_type r) 00345 { 00346 return 1; 00347 } 00348 00349 static size_type coldim(const Matrix& A, size_type r) 00350 { 00351 return 1; 00352 } 00353 00354 static size_type rowdim(const Matrix& A) 00355 { 00356 return n; 00357 } 00358 00359 static size_type coldim(const Matrix& A) 00360 { 00361 return m; 00362 } 00363 }; 00364 00365 template<typename K, int n, int m, typename TA> 00366 struct MatrixDimension<Matrix<FieldMatrix<K,n,m>, TA> > 00367 { 00368 typedef Matrix<FieldMatrix<K,n,m>, TA> ThisMatrix; 00369 typedef typename ThisMatrix::size_type size_type; 00370 00371 static size_type rowdim(const ThisMatrix& A, size_type r) 00372 { 00373 return n; 00374 } 00375 00376 static size_type coldim(const ThisMatrix& A, size_type r) 00377 { 00378 return m; 00379 } 00380 00381 static size_type rowdim(const ThisMatrix& A) 00382 { 00383 return A.N()*n; 00384 } 00385 00386 static size_type coldim(const ThisMatrix& A) 00387 { 00388 return A.M()*m; 00389 } 00390 }; 00391 00395 template<typename T> 00396 struct IsMatrix 00397 { 00398 enum{ 00402 value = false 00403 }; 00404 }; 00405 00406 template<typename T> 00407 struct IsMatrix<DenseMatrix<T> > 00408 { 00409 enum{ 00413 value = true 00414 }; 00415 }; 00416 00417 00418 template<typename T, typename A> 00419 struct IsMatrix<BCRSMatrix<T,A> > 00420 { 00421 enum{ 00425 value = true 00426 }; 00427 }; 00428 00429 } 00430 #endif