dune-istl  2.2.0
matrix.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002 // vi: set et ts=8 sw=4 sts=4:
00003 #ifndef DUNE_MATRIX_HH
00004 #define DUNE_MATRIX_HH
00005 
00010 #include <vector>
00011 #include <memory>
00012 
00013 #include <dune/istl/vbvector.hh>
00014 
00015 namespace Dune {
00016 
00022     template<class T, class A=std::allocator<T> >
00023 class Matrix
00024 {
00025 public:
00026 
00028     typedef typename T::field_type field_type;
00029 
00031     typedef T block_type;
00032 
00034     typedef A allocator_type;
00035 
00037     typedef typename VariableBlockVector<T,A>::window_type row_type;
00038 
00040     typedef typename A::size_type size_type;
00041 
00043     typedef typename VariableBlockVector<T,A>::Iterator RowIterator;
00044 
00046     typedef typename row_type::iterator ColIterator;
00047 
00049     typedef typename VariableBlockVector<T,A>::ConstIterator ConstRowIterator;
00050 
00052     typedef typename row_type::const_iterator ConstColIterator;
00053 
00054     enum {
00056         blocklevel = T::blocklevel+1
00057     };
00058 
00060     Matrix() : data_(0), cols_(0) 
00061     {}
00062 
00065     Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
00066     {}
00067 
00072     void setSize(size_type rows, size_type cols) {
00073         data_.resize(rows,cols);
00074         cols_ = cols;
00075     }
00076 
00078     RowIterator begin()
00079     {
00080         return data_.begin();
00081     }
00082 
00084     RowIterator end()
00085     {
00086         return data_.end();
00087     }
00088 
00091     RowIterator beforeEnd ()
00092     {
00093         return data_.beforeEnd();
00094     }
00095 
00098     RowIterator beforeBegin ()
00099     {
00100         return data_.beforeBegin();
00101     }
00102 
00104     ConstRowIterator begin() const
00105     {
00106         return data_.begin();
00107     }
00108 
00110     ConstRowIterator end() const
00111     {
00112         return data_.end();
00113     }
00114 
00117     ConstRowIterator beforeEnd() const
00118     {
00119         return data_.beforeEnd();
00120     }
00121 
00124     ConstRowIterator beforeBegin () const
00125     {
00126         return data_.beforeBegin();
00127     }
00128 
00130     Matrix& operator= (const field_type& t)
00131     {
00132         data_ = t;
00133       return *this;
00134     }
00135 
00137     row_type& operator[](size_type row) {
00138 #ifdef DUNE_ISTL_WITH_CHECKING
00139         if (row<0)
00140             DUNE_THROW(ISTLError, "Can't access negative rows!");
00141         if (row>=N())
00142             DUNE_THROW(ISTLError, "Row index out of range!");
00143 #endif
00144         return data_[row];
00145     }
00146 
00148     const row_type& operator[](size_type row) const {
00149 #ifdef DUNE_ISTL_WITH_CHECKING
00150         if (row<0)
00151             DUNE_THROW(ISTLError, "Can't access negative rows!");
00152         if (row>=N())
00153             DUNE_THROW(ISTLError, "Row index out of range!");
00154 #endif
00155         return data_[row];
00156     }
00157 
00159     size_type N() const {
00160         return data_.N();
00161     }
00162 
00164     size_type M() const {
00165         return cols_;
00166     }
00167 
00169     size_type rowdim() const {
00170 #ifdef DUNE_ISTL_WITH_CHECKING
00171         if (M()==0)
00172             DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
00173 #endif
00174         size_type dim = 0;
00175         for (size_type i=0; i<data_.N(); i++)
00176             dim += data_[i][0].rowdim();
00177 
00178         return dim;
00179     }
00180 
00182     size_type coldim() const {
00183 #ifdef DUNE_ISTL_WITH_CHECKING
00184         if (N()==0)
00185             DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
00186 #endif
00187         size_type dim = 0;
00188         for (size_type i=0; i<data_[0].size(); i++)
00189             dim += data_[0][i].coldim();
00190 
00191         return dim;
00192     }
00193 
00195     size_type rowdim(size_type r) const {
00196 #ifdef DUNE_ISTL_WITH_CHECKING
00197         if (r<0 || r>=N())
00198             DUNE_THROW(ISTLError, "Rowdim for nonexisting row " << r << " requested!");
00199         if (M()==0)
00200             DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
00201 #endif
00202         return data_[r][0].rowdim();
00203     }
00204 
00206     size_type coldim(size_type c) const {
00207 #ifdef DUNE_ISTL_WITH_CHECKING
00208         if (c<0 || c>=M())
00209             DUNE_THROW(ISTLError, "Coldim for nonexisting column " << c << " requested!");
00210         if (N()==0)
00211             DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
00212 #endif
00213         return data_[0][c].coldim();
00214     }
00215 
00217     Matrix<T>& operator*=(const field_type& scalar) {
00218         data_ *= scalar;
00219         return (*this);
00220     }
00221 
00223     Matrix<T>& operator/=(const field_type& scalar) {
00224         data_ /= scalar;
00225         return (*this);
00226     }
00227 
00233     Matrix& operator+= (const Matrix& b) {
00234 #ifdef DUNE_ISTL_WITH_CHECKING
00235         if(N()!=b.N() || M() != b.M())
00236             DUNE_THROW(RangeError, "Matrix sizes do not match!");
00237 #endif
00238         data_ += b.data_;
00239         return (*this);
00240     }
00241 
00247     Matrix& operator-= (const Matrix& b) {
00248 #ifdef DUNE_ISTL_WITH_CHECKING
00249         if(N()!=b.N() || M() != b.M())
00250             DUNE_THROW(RangeError, "Matrix sizes do not match!");
00251 #endif
00252         data_ -= b.data_;
00253         return (*this);
00254     }
00255 
00257     Matrix transpose() const {
00258         Matrix out(N(), M());
00259         for (size_type i=0; i<M(); i++)
00260             for (size_type j=0; j<N(); j++)
00261                 out[j][i] = (*this)[i][j];
00262 
00263         return out;
00264     }
00265 
00269     template <class X, class Y>
00270     Y transposedMult(const X& vec) DUNE_DEPRECATED;
00271     
00273     friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
00274         Matrix<T> out(m1.N(), m2.M());
00275         out = 0;
00276         
00277         for (size_type i=0; i<out.N(); i++ ) {
00278             for ( size_type j=0; j<out.M(); j++ ) 
00279                 for (size_type k=0; k<m1.M(); k++)
00280                     out[i][j] += m1[i][k]*m2[k][j];
00281         }
00282 
00283         return out;
00284     }
00285 
00287     template <class X, class Y>
00288     friend Y operator*(const Matrix<T>& m, const X& vec) {
00289 #ifdef DUNE_ISTL_WITH_CHECKING
00290         if (m.M()!=vec.size())
00291             DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix columns!");
00292 #endif
00293         Y out(m.N());
00294         out = 0;
00295 
00296         for (size_type i=0; i<out.size(); i++ ) {
00297             for ( size_type j=0; j<vec.size(); j++ ) 
00298                 out[i] += m[i][j]*vec[j];
00299         }
00300 
00301         return out;
00302     }
00303 
00305     template <class X, class Y>
00306     void mv(const X& x, Y& y) const
00307     {
00308 #ifdef DUNE_ISTL_WITH_CHECKING
00309         if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00310         if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00311 #endif
00312 
00313         for (size_type i=0; i<data_.N(); i++) {
00314           y[i]=0;
00315             for (size_type j=0; j<cols_; j++)
00316                 (*this)[i][j].umv(x[j], y[i]);
00317 
00318         }
00319 
00320     }
00321 
00323     template<class X, class Y>
00324     void mtv (const X& x, Y& y) const
00325     {
00326 #ifdef DUNE_ISTL_WITH_CHECKING
00327         if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
00328         if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
00329 #endif
00330         
00331         for(size_type i=0; i<y.N(); ++i)
00332             y[i]=0;
00333         umtv(x,y);
00334     }
00335 
00337     template <class X, class Y>
00338     void umv(const X& x, Y& y) const
00339     {
00340 #ifdef DUNE_ISTL_WITH_CHECKING
00341         if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00342         if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00343 #endif
00344 
00345         for (size_type i=0; i<data_.N(); i++) {
00346 
00347             for (size_type j=0; j<cols_; j++)
00348                 (*this)[i][j].umv(x[j], y[i]);
00349 
00350         }
00351 
00352     }
00353 
00355     template<class X, class Y>
00356     void mmv (const X& x, Y& y) const
00357     {
00358 #ifdef DUNE_ISTL_WITH_CHECKING
00359         if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00360         if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00361 #endif
00362         ConstRowIterator endi=end();
00363         for (ConstRowIterator i=begin(); i!=endi; ++i)
00364             {
00365                 ConstColIterator endj = (*i).end();
00366                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00367                     (*j).mmv(x[j.index()],y[i.index()]);
00368             }
00369     }
00370 
00372     template <class X, class Y>
00373     void usmv(const field_type& alpha, const X& x, Y& y) const
00374     {
00375 #ifdef DUNE_ISTL_WITH_CHECKING
00376         if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00377         if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00378 #endif
00379 
00380         for (size_type i=0; i<data_.N(); i++) {
00381 
00382             for (size_type j=0; j<cols_; j++)
00383                 (*this)[i][j].usmv(alpha, x[j], y[i]);
00384             
00385         }
00386         
00387     }
00388     
00390     template<class X, class Y>
00391     void umtv (const X& x, Y& y) const
00392     {
00393 #ifdef DUNE_ISTL_WITH_CHECKING
00394         if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00395         if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00396 #endif
00397         ConstRowIterator endi=end();
00398         for (ConstRowIterator i=begin(); i!=endi; ++i)
00399             {
00400                 ConstColIterator endj = (*i).end();
00401                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00402                     (*j).umtv(x[i.index()],y[j.index()]);
00403             }
00404     }
00405     
00407     template<class X, class Y>
00408     void mmtv (const X& x, Y& y) const
00409     {
00410 #ifdef DUNE_ISTL_WITH_CHECKING
00411         if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00412         if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00413 #endif
00414         ConstRowIterator endi=end();
00415         for (ConstRowIterator i=begin(); i!=endi; ++i)
00416             {
00417                 ConstColIterator endj = (*i).end();
00418                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00419                     (*j).mmtv(x[i.index()],y[j.index()]);
00420             }
00421     }
00422     
00424     template<class X, class Y>
00425     void usmtv (const field_type& alpha, const X& x, Y& y) const
00426     {
00427 #ifdef DUNE_ISTL_WITH_CHECKING
00428         if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00429         if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00430 #endif
00431         ConstRowIterator endi=end();
00432         for (ConstRowIterator i=begin(); i!=endi; ++i)
00433             {
00434                 ConstColIterator endj = (*i).end();
00435                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00436                     (*j).usmtv(alpha,x[i.index()],y[j.index()]);
00437             }
00438     }
00439     
00441     template<class X, class Y>
00442     void umhv (const X& x, Y& y) const
00443     {
00444 #ifdef DUNE_ISTL_WITH_CHECKING
00445         if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00446         if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00447 #endif
00448         ConstRowIterator endi=end();
00449         for (ConstRowIterator i=begin(); i!=endi; ++i)
00450             {
00451                 ConstColIterator endj = (*i).end();
00452                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00453                     (*j).umhv(x[i.index()],y[j.index()]);
00454             }
00455     }
00456     
00458     template<class X, class Y>
00459     void mmhv (const X& x, Y& y) const
00460     {
00461 #ifdef DUNE_ISTL_WITH_CHECKING
00462         if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00463         if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00464 #endif
00465         ConstRowIterator endi=end();
00466         for (ConstRowIterator i=begin(); i!=endi; ++i)
00467             {
00468                 ConstColIterator endj = (*i).end();
00469                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00470                     (*j).mmhv(x[i.index()],y[j.index()]);
00471             }
00472     }
00473     
00475     template<class X, class Y>
00476     void usmhv (const field_type& alpha, const X& x, Y& y) const
00477     {
00478 #ifdef DUNE_ISTL_WITH_CHECKING
00479         if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00480         if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
00481 #endif
00482         ConstRowIterator endi=end();
00483         for (ConstRowIterator i=begin(); i!=endi; ++i)
00484             {
00485                 ConstColIterator endj = (*i).end();
00486                 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
00487                     (*j).usmhv(alpha,x[i.index()],y[j.index()]);
00488             }
00489     }
00490 
00491     //===== norms
00492     
00494     double frobenius_norm () const
00495     {
00496         return std::sqrt(frobenius_norm2());
00497     }
00498     
00500     double frobenius_norm2 () const
00501     {
00502         double sum=0;
00503         for (size_type i=0; i<N(); ++i) 
00504             for (size_type j=0; j<M(); ++j)
00505                 sum += data_[i][j].frobenius_norm2();
00506         return sum;
00507     }
00508     
00510     double infinity_norm () const
00511     {
00512         double max=0;
00513         for (size_type i=0; i<N(); ++i) {
00514             double sum=0;
00515             for (size_type j=0; j<M(); j++)
00516                 sum += data_[i][j].infinity_norm();
00517             max = std::max(max,sum);
00518         }
00519         return max;
00520     }
00521     
00523     double infinity_norm_real () const
00524     {
00525         double max=0;
00526         for (size_type i=0; i<N(); ++i) {
00527             double sum=0;
00528             for (size_type j=0; j<M(); j++)
00529                 sum += data_[i][j].infinity_norm_real();
00530             max = std::max(max,sum);
00531         }
00532         return max;
00533     }
00534     
00535     //===== query
00536     
00538     bool exists (size_type i, size_type j) const
00539     {
00540 #ifdef DUNE_ISTL_WITH_CHECKING
00541         if (i<0 || i>=N()) DUNE_THROW(ISTLError,"row index out of range");
00542         if (j<0 || i>=M()) DUNE_THROW(ISTLError,"column index out of range");
00543 #endif
00544         return true;
00545     }
00546 
00547 protected:
00548 
00554     VariableBlockVector<T,A> data_;
00555 
00561     size_type cols_;
00562 };
00563 
00564 template<class T, class A>
00565 template<class X, class Y>
00566 inline Y Matrix<T,A>::transposedMult(const X& vec)
00567 {
00568 #ifdef DUNE_ISTL_WITH_CHECKING
00569     if (N()!=vec.size())
00570         DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix rows!");
00571 #endif
00572     Y out(M());
00573     out = 0;
00574     
00575     for (size_type i=0; i<out.size(); i++ ) {
00576         for ( size_type j=0; j<vec.size(); j++ ) 
00577             out[i] += (*this)[j][i]*vec[j];
00578     }
00579     
00580     return out;
00581 }
00582 
00584 } // end namespace Dune
00585 
00586 #endif