dune-istl  2.2.0
matrixutils.hh
Go to the documentation of this file.
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