dune-common
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 // $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