dune-common  2.2.0
densematrix.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 // $Id: fmatrix.hh 6128 2010-09-08 13:50:00Z christi $
00004 #ifndef DUNE_DENSEMATRIX_HH
00005 #define DUNE_DENSEMATRIX_HH
00006 
00007 #include <cmath>
00008 #include <cstddef>
00009 #include <iostream>
00010 #include <vector>
00011 
00012 #include <dune/common/misc.hh>
00013 #include <dune/common/exceptions.hh>
00014 #include <dune/common/fvector.hh>
00015 #include <dune/common/precision.hh>
00016 #include <dune/common/static_assert.hh>
00017 #include <dune/common/classname.hh>
00018 
00019 
00020 namespace Dune
00021 {
00022    
00023   template<typename M> class DenseMatrix;
00024 
00025   template<typename M>
00026   struct FieldTraits< DenseMatrix<M> >
00027   {
00028     typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::field_type field_type;
00029     typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::real_type real_type;
00030   };
00031 
00032   /*
00033     work around a problem of FieldMatrix/FieldVector,
00034     there is no unique way to obtain the size of a class
00035    */
00036   template<class K, int N, int M> class FieldMatrix;
00037   template<class K, int N> class FieldVector;
00038   namespace {
00039     template<class V>
00040     struct VectorSize
00041     {
00042       static typename V::size_type size(const V & v) { return v.size(); }
00043     };
00044 
00045     template<class K, int N>
00046     struct VectorSize< const FieldVector<K,N> >
00047     {
00048       typedef FieldVector<K,N> V;
00049       static typename V::size_type size(const V & v) { return N; }
00050     };
00051   }
00052   
00068   template<typename M, typename T>
00069   void istl_assign_to_fmatrix(DenseMatrix<M>& f, const T& t)
00070   {
00071     DUNE_THROW(NotImplemented, "You need to specialise the method istl_assign_to_fmatrix(DenseMatrix<M>& f, const T& t) "
00072                                << "(with M being " << className<M>() << ") "
00073                                << "for T == " << className<T>() << "!");
00074   }
00075 
00076   namespace
00077   {
00078     template<bool b>
00079     struct DenseMatrixAssigner
00080     {
00081       template<typename M, typename T>
00082       static void assign(DenseMatrix<M>& fm, const T& t)
00083       {
00084         istl_assign_to_fmatrix(fm, t);
00085       }
00086       
00087     };
00088     
00089 
00090     template<>
00091     struct DenseMatrixAssigner<true>
00092     {
00093       template<typename M, typename T>
00094       static void assign(DenseMatrix<M>& fm, const T& t)
00095       {
00096         fm = static_cast<const typename DenseMatVecTraits<M>::value_type>(t);
00097       }
00098     };
00099   }
00100   
00102   class FMatrixError : public Exception {};
00103 
00114   template<typename MAT>
00115   class DenseMatrix
00116   {
00117     typedef DenseMatVecTraits<MAT> Traits;
00118 
00119     // Curiously recurring template pattern
00120     MAT & asImp() { return static_cast<MAT&>(*this); }
00121     const MAT & asImp() const { return static_cast<const MAT&>(*this); }
00122     
00123   public:
00124     //===== type definitions and constants
00125 
00127     typedef typename Traits::derived_type derived_type;
00128 
00130     typedef typename Traits::value_type value_type;
00131 
00133     typedef typename Traits::value_type field_type;
00134 
00136     typedef typename Traits::value_type block_type;
00137 
00139     typedef typename Traits::size_type size_type;
00140     
00142     typedef typename Traits::row_type row_type;
00143 
00145     typedef typename Traits::row_reference row_reference;
00146 
00148     typedef typename Traits::const_row_reference const_row_reference;
00149 
00151     enum {
00153       blocklevel = 1
00154     };
00155 
00156     //===== access to components
00157     
00159     row_reference operator[] ( size_type i )
00160     {
00161       return asImp().mat_access(i);
00162     }
00163     
00164     const_row_reference operator[] ( size_type i ) const
00165     {
00166       return asImp().mat_access(i);
00167     }
00168     
00170     size_type size() const
00171     {
00172       return rows();
00173     }
00174         
00175     //===== iterator interface to rows of the matrix
00177     typedef DenseIterator<DenseMatrix,row_type> Iterator;
00179     typedef Iterator iterator;
00181     typedef Iterator RowIterator;
00183     typedef typename row_type::Iterator ColIterator;
00184 
00186     Iterator begin ()
00187     {
00188       return Iterator(*this,0);
00189     }
00190       
00192     Iterator end ()
00193     {
00194       return Iterator(*this,rows());
00195     }
00196 
00199     Iterator beforeEnd ()
00200     {
00201       return Iterator(*this,rows()-1);
00202     }
00203 
00206     Iterator beforeBegin ()
00207     {
00208       return Iterator(*this,-1);
00209     }
00210 
00212     typedef DenseIterator<const DenseMatrix,const row_type> ConstIterator;
00214     typedef ConstIterator const_iterator;
00216     typedef ConstIterator ConstRowIterator;
00218     typedef typename row_type::ConstIterator ConstColIterator;
00219 
00221     ConstIterator begin () const
00222     {
00223       return ConstIterator(*this,0);
00224     }
00225       
00227     ConstIterator end () const
00228     {
00229       return ConstIterator(*this,rows());
00230     }
00231 
00234     ConstIterator beforeEnd () const
00235     {
00236       return ConstIterator(*this,rows()-1);
00237     }
00238     
00241     ConstIterator beforeBegin () const
00242     {
00243       return ConstIterator(*this,-1);
00244     }
00245 
00246     //===== assignment from scalar
00247     DenseMatrix& operator= (const field_type& f)
00248     {
00249       for (size_type i=0; i<rows(); i++)
00250         (*this)[i] = f;
00251       return *this;   
00252     }
00253 
00254     template<typename T>
00255     DenseMatrix& operator= (const T& t)
00256     {
00257       DenseMatrixAssigner<Conversion<T,field_type>::exists>::assign(*this, t);
00258       return *this;
00259     }
00260     //===== vector space arithmetic
00261 
00263     template <class Other>
00264     DenseMatrix& operator+= (const DenseMatrix<Other>& y)
00265     {
00266       for (size_type i=0; i<rows(); i++)
00267         (*this)[i] += y[i];
00268       return *this;
00269     }
00270 
00272     template <class Other>
00273     DenseMatrix& operator-= (const DenseMatrix<Other>& y)
00274     {
00275       for (size_type i=0; i<rows(); i++)
00276         (*this)[i] -= y[i];
00277       return *this;
00278     }
00279 
00281     DenseMatrix& operator*= (const field_type& k)
00282     {
00283       for (size_type i=0; i<rows(); i++)
00284         (*this)[i] *= k;
00285       return *this;
00286     }
00287 
00289     DenseMatrix& operator/= (const field_type& k)
00290     {
00291       for (size_type i=0; i<rows(); i++)
00292         (*this)[i] /= k;
00293       return *this;
00294     }
00295 
00297     template <class Other>
00298     DenseMatrix &axpy (const field_type &k, const DenseMatrix<Other> &y )
00299     {
00300       for( size_type i = 0; i < rows(); ++i )
00301         (*this)[ i ].axpy( k, y[ i ] );
00302       return *this;
00303     }
00304 
00306     template <class Other>
00307     bool operator== (const DenseMatrix<Other>& y) const
00308     {
00309       for (size_type i=0; i<rows(); i++)
00310         if ((*this)[i]!=y[i])
00311           return false;
00312       return true;
00313     }
00315     template <class Other>
00316     bool operator!= (const DenseMatrix<Other>& y) const
00317     {
00318       return !operator==(y);
00319     }
00320     
00321 
00322     //===== linear maps
00323    
00325     template<class X, class Y>
00326     void mv (const X& x, Y& y) const
00327     {
00328 #ifdef DUNE_FMatrix_WITH_CHECKING
00329       if (x.N()!=M()) DUNE_THROW(FMatrixError,"Index out of range");
00330       if (y.N()!=N()) DUNE_THROW(FMatrixError,"Index out of range");
00331 #endif
00332       for (size_type i=0; i<rows(); ++i)
00333       {
00334         y[i] = 0;
00335         for (size_type j=0; j<cols(); j++)
00336           y[i] += (*this)[i][j] * x[j];
00337       }
00338     }
00339 
00341     template< class X, class Y >
00342     void mtv ( const X &x, Y &y ) const
00343     {
00344 #ifdef DUNE_FMatrix_WITH_CHECKING
00345       //assert( &x != &y ); 
00346       //This assert did not work for me. Compile error:
00347       //  comparison between distinct pointer types ‘const
00348       //  Dune::FieldVector<double, 3>*’ and ‘Dune::FieldVector<double, 2>*’ lacks a cast
00349       if( x.N() != N() )
00350         DUNE_THROW( FMatrixError, "Index out of range." );
00351       if( y.N() != M() )
00352         DUNE_THROW( FMatrixError, "Index out of range." );
00353 #endif
00354       for( size_type i = 0; i < cols(); ++i )
00355       {
00356         y[ i ] = 0;
00357         for( size_type j = 0; j < rows(); ++j )
00358           y[ i ] += (*this)[ j ][ i ] * x[ j ];
00359       }
00360     }
00361 
00363     template<class X, class Y>
00364     void umv (const X& x, Y& y) const
00365     {
00366 #ifdef DUNE_FMatrix_WITH_CHECKING
00367       if (x.N()!=M())
00368         DUNE_THROW(FMatrixError,"y += A x -- index out of range (sizes: x: " << x.N() << ", y: " << y.N() << ", A: " << this->N() << " x " << this->M() << ")" << std::endl);
00369       if (y.N()!=N())
00370         DUNE_THROW(FMatrixError,"y += A x -- index out of range (sizes: x: " << x.N() << ", y: " << y.N() << ", A: " << this->N() << " x " << this->M() << ")" << std::endl);
00371 #endif
00372       for (size_type i=0; i<rows(); i++)
00373         for (size_type j=0; j<cols(); j++)
00374           y[i] += (*this)[i][j] * x[j];
00375     }
00376 
00378     template<class X, class Y>
00379     void umtv (const X& x, Y& y) const
00380     {
00381 #ifdef DUNE_FMatrix_WITH_CHECKING
00382       if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00383       if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00384 #endif
00385       
00386       for (size_type i=0; i<rows(); i++)
00387         for (size_type j=0; j<cols(); j++)
00388           y[j] += (*this)[i][j]*x[i];
00389     }
00390 
00392     template<class X, class Y>
00393     void umhv (const X& x, Y& y) const
00394     {
00395 #ifdef DUNE_FMatrix_WITH_CHECKING
00396       if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00397       if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00398 #endif
00399       
00400       for (size_type i=0; i<rows(); i++)
00401         for (size_type j=0; j<cols(); j++)
00402           y[j] += conjugateComplex((*this)[i][j])*x[i];
00403     }
00404 
00406     template<class X, class Y>
00407     void mmv (const X& x, Y& y) const
00408     {
00409 #ifdef DUNE_FMatrix_WITH_CHECKING
00410       if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00411       if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00412 #endif
00413       for (size_type i=0; i<rows(); i++)
00414         for (size_type j=0; j<cols(); j++)
00415           y[i] -= (*this)[i][j] * x[j];
00416     }
00417 
00419     template<class X, class Y>
00420     void mmtv (const X& x, Y& y) const
00421     {
00422 #ifdef DUNE_FMatrix_WITH_CHECKING
00423       if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00424       if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00425 #endif
00426       
00427       for (size_type i=0; i<rows(); i++)
00428         for (size_type j=0; j<cols(); j++)
00429           y[j] -= (*this)[i][j]*x[i];
00430     }
00431 
00433     template<class X, class Y>
00434     void mmhv (const X& x, Y& y) const
00435     {
00436 #ifdef DUNE_FMatrix_WITH_CHECKING
00437       if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00438       if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00439 #endif
00440       
00441       for (size_type i=0; i<rows(); i++)
00442         for (size_type j=0; j<cols(); j++)
00443           y[j] -= conjugateComplex((*this)[i][j])*x[i];
00444     }
00445 
00447     template<class X, class Y>
00448     void usmv (const field_type& alpha, const X& x, Y& y) const
00449     {
00450 #ifdef DUNE_FMatrix_WITH_CHECKING
00451       if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00452       if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00453 #endif
00454       for (size_type i=0; i<rows(); i++)
00455         for (size_type j=0; j<cols(); j++)
00456           y[i] += alpha * (*this)[i][j] * x[j];
00457     }
00458 
00460     template<class X, class Y>
00461     void usmtv (const field_type& alpha, const X& x, Y& y) const
00462     {
00463 #ifdef DUNE_FMatrix_WITH_CHECKING
00464       if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00465       if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00466 #endif
00467       
00468       for (size_type i=0; i<rows(); i++)
00469         for (size_type j=0; j<cols(); j++)
00470           y[j] += alpha*(*this)[i][j]*x[i];
00471     }
00472 
00474     template<class X, class Y>
00475     void usmhv (const field_type& alpha, const X& x, Y& y) const
00476     {
00477 #ifdef DUNE_FMatrix_WITH_CHECKING
00478       if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00479       if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00480 #endif
00481       
00482       for (size_type i=0; i<rows(); i++)
00483         for (size_type j=0; j<cols(); j++)
00484           y[j] += alpha*conjugateComplex((*this)[i][j])*x[i];
00485     }
00486 
00487     //===== norms
00488 
00490     typename FieldTraits<value_type>::real_type frobenius_norm () const
00491     {
00492       typename FieldTraits<value_type>::real_type sum=(0.0);
00493       for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
00494       return fvmeta::sqrt(sum);
00495     }
00496 
00498     typename FieldTraits<value_type>::real_type frobenius_norm2 () const
00499     {
00500       typename FieldTraits<value_type>::real_type sum=(0.0);
00501       for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
00502       return sum;
00503     }
00504 
00506     typename FieldTraits<value_type>::real_type infinity_norm () const
00507     {
00508       typename remove_const< typename FieldTraits<value_type>::real_type >::type max=(0.0);
00509       for (size_type i=0; i<rows(); ++i) max = std::max(max,(*this)[i].one_norm());
00510       return max;
00511     }
00512 
00514     typename FieldTraits<value_type>::real_type infinity_norm_real () const
00515     {
00516       typename FieldTraits<value_type>::real_type max(0.0);
00517       for (size_type i=0; i<rows(); ++i) max = std::max(max,(*this)[i].one_norm_real());
00518       return max;
00519     }
00520 
00521     //===== solve
00522 
00527     template <class V>
00528     void solve (V& x, const V& b) const;
00529 
00534     void invert();
00535 
00537     field_type determinant () const;
00538 
00540     template<typename M2>
00541     MAT& leftmultiply (const DenseMatrix<M2>& M)
00542     {
00543       assert(M.rows() == M.cols() && M.rows() == rows());
00544       MAT C(asImp());
00545 
00546       for (size_type i=0; i<rows(); i++)
00547         for (size_type j=0; j<cols(); j++) {
00548           (*this)[i][j] = 0;
00549           for (size_type k=0; k<rows(); k++)
00550             (*this)[i][j] += M[i][k]*C[k][j];
00551         }
00552       
00553       return asImp();
00554     }
00555 
00557     template<typename M2>
00558     MAT& rightmultiply (const DenseMatrix<M2>& M)
00559     {
00560       assert(M.rows() == M.cols() && M.cols() == cols());
00561       MAT C(asImp());
00562       
00563       for (size_type i=0; i<rows(); i++)
00564         for (size_type j=0; j<cols(); j++) {
00565           (*this)[i][j] = 0;
00566           for (size_type k=0; k<cols(); k++)
00567             (*this)[i][j] += C[i][k]*M[k][j];
00568         }
00569       return asImp();
00570     }
00571 
00572 #if 0
00573 
00574     template<int l>
00575     DenseMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
00576     {
00577       FieldMatrix<K,l,cols> C;
00578       
00579       for (size_type i=0; i<l; i++) {
00580         for (size_type j=0; j<cols(); j++) {
00581           C[i][j] = 0;
00582           for (size_type k=0; k<rows(); k++)
00583             C[i][j] += M[i][k]*(*this)[k][j];
00584         }
00585       }
00586       return C;
00587     }
00588 
00590     template<int l>
00591     FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
00592     {
00593       FieldMatrix<K,rows,l> C;
00594       
00595       for (size_type i=0; i<rows(); i++) {
00596         for (size_type j=0; j<l; j++) {
00597           C[i][j] = 0;
00598           for (size_type k=0; k<cols(); k++)
00599             C[i][j] += (*this)[i][k]*M[k][j];
00600         }
00601       }
00602       return C;
00603     }
00604 #endif
00605 
00606     //===== sizes
00607 
00609     size_type N () const
00610     {
00611       return rows();
00612     }
00613 
00615     size_type M () const
00616     {
00617       return cols();
00618     }
00619 
00621     size_type rows() const
00622     {
00623       return asImp().mat_rows();
00624     }
00625 
00627     size_type cols() const
00628     {
00629       return asImp().mat_cols();
00630     }
00631 
00632     //===== query
00633     
00635     bool exists (size_type i, size_type j) const
00636     {
00637 #ifdef DUNE_FMatrix_WITH_CHECKING
00638       if (i<0 || i>=rows()) DUNE_THROW(FMatrixError,"row index out of range");
00639       if (j<0 || j>=cols()) DUNE_THROW(FMatrixError,"column index out of range");
00640 #endif
00641       return true;
00642     }
00643 
00644   private:
00645 
00646 #ifndef DOXYGEN
00647     struct ElimPivot
00648     {
00649       ElimPivot(std::vector<size_type> & pivot);
00650       
00651       void swap(int i, int j);
00652       
00653       template<typename T>
00654       void operator()(const T&, int k, int i)
00655       {}
00656       
00657       std::vector<size_type> & pivot_;
00658     };
00659 
00660     template<typename V>
00661     struct Elim
00662     {
00663       Elim(V& rhs);
00664       
00665       void swap(int i, int j);
00666       
00667       void operator()(const typename V::field_type& factor, int k, int i);
00668 
00669       V* rhs_;
00670     };
00671 
00672     struct ElimDet
00673     {
00674       ElimDet(field_type& sign) : sign_(sign)
00675       { sign_ = 1; }
00676 
00677       void swap(int i, int j)
00678       { sign_ *= -1; }
00679 
00680       void operator()(const field_type&, int k, int i)
00681       {}
00682 
00683       field_type& sign_;
00684     };
00685 #endif // DOXYGEN
00686     
00687     template<class Func>
00688     void luDecomposition(DenseMatrix<MAT>& A, Func func) const;
00689   };
00690 
00691 #ifndef DOXYGEN
00692   template<typename MAT>
00693   DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
00694     : pivot_(pivot)
00695   {
00696     typedef typename std::vector<size_type>::size_type size_type;
00697     for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
00698   }
00699 
00700   template<typename MAT>
00701   void DenseMatrix<MAT>::ElimPivot::swap(int i, int j)
00702   {
00703     pivot_[i]=j;
00704   }
00705   
00706   template<typename MAT>
00707   template<typename V>
00708   DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
00709     : rhs_(&rhs)
00710   {}
00711   
00712   template<typename MAT>
00713   template<typename V>
00714   void DenseMatrix<MAT>::Elim<V>::swap(int i, int j)
00715   {
00716     std::swap((*rhs_)[i], (*rhs_)[j]);
00717   }
00718 
00719   template<typename MAT>
00720   template<typename V>
00721   void DenseMatrix<MAT>::
00722   Elim<V>::operator()(const typename V::field_type& factor, int k, int i)
00723   {
00724     (*rhs_)[k] -= factor*(*rhs_)[i];
00725   }
00726   template<typename MAT>
00727   template<typename Func>
00728   inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func) const
00729   {
00730     typedef typename FieldTraits<value_type>::real_type
00731       real_type;
00732     real_type norm = A.infinity_norm_real(); // for relative thresholds
00733     real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() );
00734     real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() );
00735   
00736     // LU decomposition of A in A
00737     for (size_type i=0; i<rows(); i++)  // loop over all rows
00738     {
00739       typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]);
00740       
00741       // pivoting ?
00742       if (pivmax<pivthres)
00743       {
00744         // compute maximum of column
00745         size_type imax=i;
00746         typename FieldTraits<value_type>::real_type abs(0.0);
00747         for (size_type k=i+1; k<rows(); k++)
00748           if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
00749           {
00750             pivmax = abs; imax = k;
00751           }
00752         // swap rows
00753         if (imax!=i){
00754           for (size_type j=0; j<rows(); j++)
00755             std::swap(A[i][j],A[imax][j]);
00756           func.swap(i, imax); // swap the pivot or rhs
00757         }
00758       }
00759     
00760       // singular ?
00761       if (pivmax<singthres)
00762         DUNE_THROW(FMatrixError,"matrix is singular");          
00763     
00764       // eliminate
00765       for (size_type k=i+1; k<rows(); k++)
00766       {
00767         field_type factor = A[k][i]/A[i][i];
00768         A[k][i] = factor;
00769         for (size_type j=i+1; j<rows(); j++)
00770           A[k][j] -= factor*A[i][j];
00771         func(factor, k, i);
00772       }
00773     }
00774   }
00775 
00776   template<typename MAT>
00777   template <class V>
00778   inline void DenseMatrix<MAT>::solve(V& x, const V& b) const
00779   {
00780     // never mind those ifs, because they get optimized away
00781     if (rows()!=cols())
00782       DUNE_THROW(FMatrixError, "Can't solve for a " << rows() << "x" << cols() << " matrix!");
00783 
00784     if (rows()==1) {
00785 
00786 #ifdef DUNE_FMatrix_WITH_CHECKING
00787       if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit())
00788         DUNE_THROW(FMatrixError,"matrix is singular");          
00789 #endif
00790       x[0] = b[0]/(*this)[0][0];
00791 
00792     }
00793     else if (rows()==2) {
00794             
00795       field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
00796 #ifdef DUNE_FMatrix_WITH_CHECKING
00797       if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
00798         DUNE_THROW(FMatrixError,"matrix is singular");
00799 #endif
00800       detinv = 1.0/detinv;
00801             
00802       x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
00803       x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
00804 
00805     }
00806     else if (rows()==3) {
00807 
00808       field_type d = determinant();
00809 #ifdef DUNE_FMatrix_WITH_CHECKING
00810       if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit())
00811         DUNE_THROW(FMatrixError,"matrix is singular");
00812 #endif
00813 
00814       x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
00815         - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
00816         + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
00817 
00818       x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
00819         - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
00820         + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
00821 
00822       x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
00823         - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
00824         + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
00825 
00826     }
00827     else {
00828 
00829       V& rhs = x; // use x to store rhs
00830       rhs = b; // copy data
00831       Elim<V> elim(rhs);
00832       MAT A(asImp());
00833       
00834       luDecomposition(A, elim);
00835       
00836       // backsolve
00837       for(int i=rows()-1; i>=0; i--){
00838         for (size_type j=i+1; j<rows(); j++)
00839           rhs[i] -= A[i][j]*x[j];
00840         x[i] = rhs[i]/A[i][i];
00841       }
00842     }   
00843   }
00844 
00845   template<typename MAT>
00846   inline void DenseMatrix<MAT>::invert()
00847   {
00848     // never mind those ifs, because they get optimized away
00849     if (rows()!=cols())
00850       DUNE_THROW(FMatrixError, "Can't invert a " << rows() << "x" << cols() << " matrix!");
00851 
00852     if (rows()==1) {
00853       
00854 #ifdef DUNE_FMatrix_WITH_CHECKING
00855       if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit())
00856         DUNE_THROW(FMatrixError,"matrix is singular");        
00857 #endif
00858       (*this)[0][0] = 1.0/(*this)[0][0];
00859       
00860     }
00861     else if (rows()==2) {
00862 
00863       field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
00864 #ifdef DUNE_FMatrix_WITH_CHECKING
00865       if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
00866         DUNE_THROW(FMatrixError,"matrix is singular");        
00867 #endif
00868       detinv = 1.0/detinv;
00869 
00870       field_type temp=(*this)[0][0];
00871       (*this)[0][0] =  (*this)[1][1]*detinv;
00872       (*this)[0][1] = -(*this)[0][1]*detinv;
00873       (*this)[1][0] = -(*this)[1][0]*detinv;
00874       (*this)[1][1] =  temp*detinv;
00875 
00876     }
00877     else {
00878       
00879       MAT A(asImp());
00880       std::vector<size_type> pivot(rows());
00881       luDecomposition(A, ElimPivot(pivot));
00882       DenseMatrix<MAT>& L=A;
00883       DenseMatrix<MAT>& U=A;
00884           
00885       // initialize inverse
00886       *this=field_type();
00887         
00888       for(size_type i=0; i<rows(); ++i)
00889         (*this)[i][i]=1;
00890         
00891       // L Y = I; multiple right hand sides
00892       for (size_type i=0; i<rows(); i++)
00893         for (size_type j=0; j<i; j++)
00894           for (size_type k=0; k<rows(); k++)
00895             (*this)[i][k] -= L[i][j]*(*this)[j][k];
00896   
00897       // U A^{-1} = Y
00898       for (size_type i=rows(); i>0;){
00899         --i;
00900         for (size_type k=0; k<rows(); k++){
00901           for (size_type j=i+1; j<rows(); j++)
00902             (*this)[i][k] -= U[i][j]*(*this)[j][k];
00903           (*this)[i][k] /= U[i][i];
00904         }
00905       }
00906 
00907       for(size_type i=rows(); i>0; ){
00908         --i;
00909         if(i!=pivot[i])
00910           for(size_type j=0; j<rows(); ++j)
00911             std::swap((*this)[j][pivot[i]], (*this)[j][i]);
00912       }
00913     }
00914   }
00915 
00916   // implementation of the determinant 
00917   template<typename MAT>
00918   inline typename DenseMatrix<MAT>::field_type
00919   DenseMatrix<MAT>::determinant() const
00920   {
00921     // never mind those ifs, because they get optimized away
00922     if (rows()!=cols())
00923       DUNE_THROW(FMatrixError, "There is no determinant for a " << rows() << "x" << cols() << " matrix!");
00924 
00925     if (rows()==1)
00926       return (*this)[0][0];
00927     
00928     if (rows()==2)
00929       return (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0]; 
00930 
00931     if (rows()==3) {
00932       // code generated by maple 
00933       field_type t4  = (*this)[0][0] * (*this)[1][1];
00934       field_type t6  = (*this)[0][0] * (*this)[1][2];
00935       field_type t8  = (*this)[0][1] * (*this)[1][0];
00936       field_type t10 = (*this)[0][2] * (*this)[1][0];
00937       field_type t12 = (*this)[0][1] * (*this)[2][0];
00938       field_type t14 = (*this)[0][2] * (*this)[2][0];
00939         
00940       return (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
00941         t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
00942 
00943     }
00944 
00945     MAT A(asImp());
00946     field_type det;
00947     try
00948     {
00949       luDecomposition(A, ElimDet(det));
00950     }
00951     catch (FMatrixError&)
00952     {
00953       return 0;
00954     }
00955     for (size_type i = 0; i < rows(); ++i)
00956       det *= A[i][i];
00957     return det;
00958   }
00959 
00960 #endif // DOXYGEN
00961 
00962   namespace DenseMatrixHelp {
00963 #if 0
00964 
00965     template <typename K>
00966     static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
00967     {
00968       inverse[0][0] = 1.0/matrix[0][0];
00969       return matrix[0][0];
00970     }
00971 
00973     template <typename K>
00974     static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
00975     {
00976       return invertMatrix(matrix,inverse); 
00977     }
00978 
00979 
00981     template <typename K>
00982     static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
00983     {
00984       // code generated by maple 
00985       field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
00986       field_type det_1 = 1.0/det;
00987       inverse[0][0] =   matrix[1][1] * det_1;
00988       inverse[0][1] = - matrix[0][1] * det_1;
00989       inverse[1][0] = - matrix[1][0] * det_1;
00990       inverse[1][1] =   matrix[0][0] * det_1;
00991       return det;
00992     }
00993 
00996     template <typename K>
00997     static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
00998     {
00999       // code generated by maple 
01000       field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
01001       field_type det_1 = 1.0/det;
01002       inverse[0][0] =   matrix[1][1] * det_1;
01003       inverse[1][0] = - matrix[0][1] * det_1;
01004       inverse[0][1] = - matrix[1][0] * det_1;
01005       inverse[1][1] =   matrix[0][0] * det_1;
01006       return det;
01007     }
01008 
01010     template <typename K>
01011     static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
01012     {
01013       // code generated by maple 
01014       field_type t4  = matrix[0][0] * matrix[1][1];
01015       field_type t6  = matrix[0][0] * matrix[1][2];
01016       field_type t8  = matrix[0][1] * matrix[1][0];
01017       field_type t10 = matrix[0][2] * matrix[1][0];
01018       field_type t12 = matrix[0][1] * matrix[2][0];
01019       field_type t14 = matrix[0][2] * matrix[2][0];
01020 
01021       field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
01022         t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
01023       field_type t17 = 1.0/det;
01024 
01025       inverse[0][0] =  (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
01026       inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
01027       inverse[0][2] =  (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
01028       inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
01029       inverse[1][1] =  (matrix[0][0] * matrix[2][2] - t14) * t17;
01030       inverse[1][2] = -(t6-t10) * t17;
01031       inverse[2][0] =  (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
01032       inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
01033       inverse[2][2] =  (t4-t8) * t17;
01034 
01035       return det;
01036     }
01037 
01039     template <typename K>
01040     static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
01041     {
01042       // code generated by maple 
01043       field_type t4  = matrix[0][0] * matrix[1][1];
01044       field_type t6  = matrix[0][0] * matrix[1][2];
01045       field_type t8  = matrix[0][1] * matrix[1][0];
01046       field_type t10 = matrix[0][2] * matrix[1][0];
01047       field_type t12 = matrix[0][1] * matrix[2][0];
01048       field_type t14 = matrix[0][2] * matrix[2][0];
01049 
01050       field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
01051         t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
01052       field_type t17 = 1.0/det;
01053 
01054       inverse[0][0] =  (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
01055       inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
01056       inverse[2][0] =  (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
01057       inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
01058       inverse[1][1] =  (matrix[0][0] * matrix[2][2] - t14) * t17;
01059       inverse[2][1] = -(t6-t10) * t17;
01060       inverse[0][2] =  (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
01061       inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
01062       inverse[2][2] =  (t4-t8) * t17;
01063 
01064       return det;
01065     }
01066 
01068     template< class K, int m, int n, int p >
01069     static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
01070       const FieldMatrix< K, n, p > &B,
01071       FieldMatrix< K, m, p > &ret )
01072     {
01073       typedef typename FieldMatrix< K, m, p > :: size_type size_type;
01074 
01075       for( size_type i = 0; i < m; ++i )
01076       {
01077         for( size_type j = 0; j < p; ++j )
01078         {
01079           ret[ i ][ j ] = K( 0 );
01080           for( size_type k = 0; k < n; ++k )
01081             ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
01082         }
01083       }
01084     }
01085 
01087     template <typename K, int rows, int cols>
01088     static inline void multTransposedMatrix(const FieldMatrix<K,rows,cols> &matrix, FieldMatrix<K,cols,cols>& ret)
01089     {
01090       typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
01091   
01092       for(size_type i=0; i<cols(); i++)
01093         for(size_type j=0; j<cols(); j++)
01094         {
01095           ret[i][j]=0.0;
01096           for(size_type k=0; k<rows(); k++)
01097             ret[i][j]+=matrix[k][i]*matrix[k][j];
01098         }
01099     }
01100 #endif
01101 
01103     template <typename MAT, typename V1, typename V2>
01104     static inline void multAssign(const DenseMatrix<MAT> &matrix, const DenseVector<V1> & x, DenseVector<V2> & ret) 
01105     {
01106       assert(x.size() == matrix.cols());
01107       assert(ret.size() == matrix.rows());
01108       typedef typename DenseMatrix<MAT>::size_type size_type;
01109   
01110       for(size_type i=0; i<matrix.rows(); ++i)
01111       {
01112         ret[i] = 0.0;
01113         for(size_type j=0; j<matrix.cols(); ++j)
01114         {
01115           ret[i] += matrix[i][j]*x[j];
01116         }
01117       }
01118     }
01119 
01120 #if 0
01121 
01122     template <typename K, int rows, int cols>
01123     static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret) 
01124     {
01125       typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
01126   
01127       for(size_type i=0; i<cols(); ++i)
01128       {
01129         ret[i] = 0.0;
01130         for(size_type j=0; j<rows(); ++j)
01131           ret[i] += matrix[j][i]*x[j];
01132       }
01133     }
01134 
01136     template <typename K, int rows, int cols>
01137     static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x) 
01138     {
01139       FieldVector<K,rows> ret;
01140       multAssign(matrix,x,ret);
01141       return ret;
01142     }
01143 
01145     template <typename K, int rows, int cols>
01146     static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x) 
01147     {
01148       FieldVector<K,cols> ret;
01149       multAssignTransposed( matrix, x, ret );
01150       return ret; 
01151     }
01152 #endif
01153 
01154   } // end namespace DenseMatrixHelp 
01155 
01157   template<typename MAT>
01158   std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
01159   {
01160     for (typename DenseMatrix<MAT>::size_type i=0; i<a.rows(); i++)
01161       s << a[i] << std::endl;
01162     return s;
01163   }
01164 
01167 } // end namespace Dune
01168 
01169 #endif