00001
00002
00003 #ifndef DUNE_DENSEMATRIX_HH
00004 #define DUNE_DENSEMATRIX_HH
00005
00006 #include <cmath>
00007 #include <cstddef>
00008 #include <iostream>
00009 #include <type_traits>
00010 #include <vector>
00011
00012 #include <dune/common/boundschecking.hh>
00013 #include <dune/common/exceptions.hh>
00014 #include <dune/common/fvector.hh>
00015 #include <dune/common/precision.hh>
00016 #include <dune/common/classname.hh>
00017 #include <dune/common/math.hh>
00018 #include <dune/common/typetraits.hh>
00019 #include <dune/common/unused.hh>
00020
00021
00022 namespace Dune
00023 {
00024
00025 template<typename M> class DenseMatrix;
00026
00027 template<typename M>
00028 struct FieldTraits< DenseMatrix<M> >
00029 {
00030 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::field_type field_type;
00031 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::real_type real_type;
00032 };
00033
00034
00035
00036
00037
00038 template<class K, int N, int M> class FieldMatrix;
00039 template<class K, int N> class FieldVector;
00040 namespace {
00041 template<class V>
00042 struct VectorSize
00043 {
00044 static typename V::size_type size(const V & v) { return v.size(); }
00045 };
00046
00047 template<class K, int N>
00048 struct VectorSize< const FieldVector<K,N> >
00049 {
00050 typedef FieldVector<K,N> V;
00051 static typename V::size_type size(const V & v) { return N; }
00052 };
00053 }
00054
00073 template< class DenseMatrix, class RHS >
00074 struct DenseMatrixAssigner;
00075
00076 #ifndef DOXYGEN
00077 namespace Impl
00078 {
00079
00080 template< class DenseMatrix, class RHS, class = void >
00081 class DenseMatrixAssigner
00082 {};
00083
00084 template< class DenseMatrix, class RHS >
00085 class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< Dune::IsNumber< RHS >::value > >
00086 {
00087 public:
00088 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
00089 {
00090 typedef typename DenseMatrix::field_type field_type;
00091 std::fill( denseMatrix.begin(), denseMatrix.end(), static_cast< field_type >( rhs ) );
00092 }
00093 };
00094
00095 template< class DenseMatrix, class RHS >
00096 class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< !std::is_same< typename RHS::const_iterator, void >::value > >
00097 {
00098 public:
00099 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
00100 {
00101 DUNE_ASSERT_BOUNDS(rhs.N() == denseMatrix.N());
00102 DUNE_ASSERT_BOUNDS(rhs.M() == denseMatrix.M());
00103 typename DenseMatrix::iterator tIt = std::begin(denseMatrix);
00104 typename RHS::const_iterator sIt = std::begin(rhs);
00105 for(; sIt != std::end(rhs); ++tIt, ++sIt)
00106 std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt));
00107 }
00108 };
00109
00110 }
00111
00112
00113
00114 template< class DenseMatrix, class RHS >
00115 struct DenseMatrixAssigner
00116 : public Impl::DenseMatrixAssigner< DenseMatrix, RHS >
00117 {};
00118
00119
00120 namespace Impl
00121 {
00122
00123 template< class DenseMatrix, class RHS >
00124 std::true_type hasDenseMatrixAssigner ( DenseMatrix &, const RHS &, decltype( Dune::DenseMatrixAssigner< DenseMatrix, RHS >::apply( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) ) * = nullptr );
00125
00126 std::false_type hasDenseMatrixAssigner ( ... );
00127
00128 }
00129
00130 template< class DenseMatrix, class RHS >
00131 struct HasDenseMatrixAssigner
00132 : public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) )
00133 {};
00134
00135 #endif // #ifndef DOXYGEN
00136
00137
00138
00140 class FMatrixError : public MathError {};
00141
00152 template<typename MAT>
00153 class DenseMatrix
00154 {
00155 typedef DenseMatVecTraits<MAT> Traits;
00156
00157
00158 MAT & asImp() { return static_cast<MAT&>(*this); }
00159 const MAT & asImp() const { return static_cast<const MAT&>(*this); }
00160
00161 public:
00162
00163
00165 typedef typename Traits::derived_type derived_type;
00166
00168 typedef typename Traits::value_type value_type;
00169
00171 typedef typename Traits::value_type field_type;
00172
00174 typedef typename Traits::value_type block_type;
00175
00177 typedef typename Traits::size_type size_type;
00178
00180 typedef typename Traits::row_type row_type;
00181
00183 typedef typename Traits::row_reference row_reference;
00184
00186 typedef typename Traits::const_row_reference const_row_reference;
00187
00189 enum {
00191 blocklevel = 1
00192 };
00193
00194
00195
00197 row_reference operator[] ( size_type i )
00198 {
00199 return asImp().mat_access(i);
00200 }
00201
00202 const_row_reference operator[] ( size_type i ) const
00203 {
00204 return asImp().mat_access(i);
00205 }
00206
00208 size_type size() const
00209 {
00210 return rows();
00211 }
00212
00213
00215 typedef DenseIterator<DenseMatrix,row_type,row_reference> Iterator;
00217 typedef Iterator iterator;
00219 typedef Iterator RowIterator;
00221 typedef typename std::remove_reference<row_reference>::type::Iterator ColIterator;
00222
00224 Iterator begin ()
00225 {
00226 return Iterator(*this,0);
00227 }
00228
00230 Iterator end ()
00231 {
00232 return Iterator(*this,rows());
00233 }
00234
00237 Iterator beforeEnd ()
00238 {
00239 return Iterator(*this,rows()-1);
00240 }
00241
00244 Iterator beforeBegin ()
00245 {
00246 return Iterator(*this,-1);
00247 }
00248
00250 typedef DenseIterator<const DenseMatrix,const row_type,const_row_reference> ConstIterator;
00252 typedef ConstIterator const_iterator;
00254 typedef ConstIterator ConstRowIterator;
00256 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator ConstColIterator;
00257
00259 ConstIterator begin () const
00260 {
00261 return ConstIterator(*this,0);
00262 }
00263
00265 ConstIterator end () const
00266 {
00267 return ConstIterator(*this,rows());
00268 }
00269
00272 ConstIterator beforeEnd () const
00273 {
00274 return ConstIterator(*this,rows()-1);
00275 }
00276
00279 ConstIterator beforeBegin () const
00280 {
00281 return ConstIterator(*this,-1);
00282 }
00283
00284
00285
00286 template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > >
00287 DenseMatrix &operator= ( const RHS &rhs )
00288 {
00289 DenseMatrixAssigner< MAT, RHS >::apply( asImp(), rhs );
00290 return *this;
00291 }
00292
00293
00294
00296 template <class Other>
00297 DenseMatrix& operator+= (const DenseMatrix<Other>& y)
00298 {
00299 DUNE_ASSERT_BOUNDS(rows() == y.rows());
00300 for (size_type i=0; i<rows(); i++)
00301 (*this)[i] += y[i];
00302 return *this;
00303 }
00304
00306 template <class Other>
00307 DenseMatrix& operator-= (const DenseMatrix<Other>& y)
00308 {
00309 DUNE_ASSERT_BOUNDS(rows() == y.rows());
00310 for (size_type i=0; i<rows(); i++)
00311 (*this)[i] -= y[i];
00312 return *this;
00313 }
00314
00316 DenseMatrix& operator*= (const field_type& k)
00317 {
00318 for (size_type i=0; i<rows(); i++)
00319 (*this)[i] *= k;
00320 return *this;
00321 }
00322
00324 DenseMatrix& operator/= (const field_type& k)
00325 {
00326 for (size_type i=0; i<rows(); i++)
00327 (*this)[i] /= k;
00328 return *this;
00329 }
00330
00332 template <class Other>
00333 DenseMatrix &axpy (const field_type &k, const DenseMatrix<Other> &y )
00334 {
00335 DUNE_ASSERT_BOUNDS(rows() == y.rows());
00336 for( size_type i = 0; i < rows(); ++i )
00337 (*this)[ i ].axpy( k, y[ i ] );
00338 return *this;
00339 }
00340
00342 template <class Other>
00343 bool operator== (const DenseMatrix<Other>& y) const
00344 {
00345 DUNE_ASSERT_BOUNDS(rows() == y.rows());
00346 for (size_type i=0; i<rows(); i++)
00347 if ((*this)[i]!=y[i])
00348 return false;
00349 return true;
00350 }
00352 template <class Other>
00353 bool operator!= (const DenseMatrix<Other>& y) const
00354 {
00355 return !operator==(y);
00356 }
00357
00358
00359
00360
00362 template<class X, class Y>
00363 void mv (const X& x, Y& y) const
00364 {
00365 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
00366 DUNE_ASSERT_BOUNDS(x.N() == M());
00367 DUNE_ASSERT_BOUNDS(y.N() == N());
00368
00369 using field_type = typename FieldTraits<Y>::field_type;
00370 for (size_type i=0; i<rows(); ++i)
00371 {
00372 y[i] = field_type(0);
00373 for (size_type j=0; j<cols(); j++)
00374 y[i] += (*this)[i][j] * x[j];
00375 }
00376 }
00377
00379 template< class X, class Y >
00380 void mtv ( const X &x, Y &y ) const
00381 {
00382 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
00383 DUNE_ASSERT_BOUNDS(x.N() == N());
00384 DUNE_ASSERT_BOUNDS(y.N() == M());
00385
00386 using field_type = typename FieldTraits<Y>::field_type;
00387 for(size_type i = 0; i < cols(); ++i)
00388 {
00389 y[i] = field_type(0);
00390 for(size_type j = 0; j < rows(); ++j)
00391 y[i] += (*this)[j][i] * x[j];
00392 }
00393 }
00394
00396 template<class X, class Y>
00397 void umv (const X& x, Y& y) const
00398 {
00399 DUNE_ASSERT_BOUNDS(x.N() == M());
00400 DUNE_ASSERT_BOUNDS(y.N() == N());
00401 for (size_type i=0; i<rows(); i++)
00402 for (size_type j=0; j<cols(); j++)
00403 y[i] += (*this)[i][j] * x[j];
00404 }
00405
00407 template<class X, class Y>
00408 void umtv (const X& x, Y& y) const
00409 {
00410 DUNE_ASSERT_BOUNDS(x.N() == N());
00411 DUNE_ASSERT_BOUNDS(y.N() == M());
00412 for (size_type i=0; i<rows(); i++)
00413 for (size_type j=0; j<cols(); j++)
00414 y[j] += (*this)[i][j]*x[i];
00415 }
00416
00418 template<class X, class Y>
00419 void umhv (const X& x, Y& y) const
00420 {
00421 DUNE_ASSERT_BOUNDS(x.N() == N());
00422 DUNE_ASSERT_BOUNDS(y.N() == M());
00423 for (size_type i=0; i<rows(); i++)
00424 for (size_type j=0; j<cols(); j++)
00425 y[j] += conjugateComplex((*this)[i][j])*x[i];
00426 }
00427
00429 template<class X, class Y>
00430 void mmv (const X& x, Y& y) const
00431 {
00432 DUNE_ASSERT_BOUNDS(x.N() == M());
00433 DUNE_ASSERT_BOUNDS(y.N() == N());
00434 for (size_type i=0; i<rows(); i++)
00435 for (size_type j=0; j<cols(); j++)
00436 y[i] -= (*this)[i][j] * x[j];
00437 }
00438
00440 template<class X, class Y>
00441 void mmtv (const X& x, Y& y) const
00442 {
00443 DUNE_ASSERT_BOUNDS(x.N() == N());
00444 DUNE_ASSERT_BOUNDS(y.N() == M());
00445 for (size_type i=0; i<rows(); i++)
00446 for (size_type j=0; j<cols(); j++)
00447 y[j] -= (*this)[i][j]*x[i];
00448 }
00449
00451 template<class X, class Y>
00452 void mmhv (const X& x, Y& y) const
00453 {
00454 DUNE_ASSERT_BOUNDS(x.N() == N());
00455 DUNE_ASSERT_BOUNDS(y.N() == M());
00456 for (size_type i=0; i<rows(); i++)
00457 for (size_type j=0; j<cols(); j++)
00458 y[j] -= conjugateComplex((*this)[i][j])*x[i];
00459 }
00460
00462 template<class X, class Y>
00463 void usmv (const typename FieldTraits<Y>::field_type & alpha,
00464 const X& x, Y& y) const
00465 {
00466 DUNE_ASSERT_BOUNDS(x.N() == M());
00467 DUNE_ASSERT_BOUNDS(y.N() == N());
00468 for (size_type i=0; i<rows(); i++)
00469 for (size_type j=0; j<cols(); j++)
00470 y[i] += alpha * (*this)[i][j] * x[j];
00471 }
00472
00474 template<class X, class Y>
00475 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
00476 const X& x, Y& y) const
00477 {
00478 DUNE_ASSERT_BOUNDS(x.N() == N());
00479 DUNE_ASSERT_BOUNDS(y.N() == M());
00480 for (size_type i=0; i<rows(); i++)
00481 for (size_type j=0; j<cols(); j++)
00482 y[j] += alpha*(*this)[i][j]*x[i];
00483 }
00484
00486 template<class X, class Y>
00487 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
00488 const X& x, Y& y) const
00489 {
00490 DUNE_ASSERT_BOUNDS(x.N() == N());
00491 DUNE_ASSERT_BOUNDS(y.N() == M());
00492 for (size_type i=0; i<rows(); i++)
00493 for (size_type j=0; j<cols(); j++)
00494 y[j] += alpha*conjugateComplex((*this)[i][j])*x[i];
00495 }
00496
00497
00498
00500 typename FieldTraits<value_type>::real_type frobenius_norm () const
00501 {
00502 typename FieldTraits<value_type>::real_type sum=(0.0);
00503 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
00504 return fvmeta::sqrt(sum);
00505 }
00506
00508 typename FieldTraits<value_type>::real_type frobenius_norm2 () const
00509 {
00510 typename FieldTraits<value_type>::real_type sum=(0.0);
00511 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
00512 return sum;
00513 }
00514
00516 template <typename vt = value_type,
00517 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
00518 typename FieldTraits<vt>::real_type infinity_norm() const {
00519 using real_type = typename FieldTraits<vt>::real_type;
00520 using std::max;
00521
00522 real_type norm = 0;
00523 for (auto const &x : *this) {
00524 real_type const a = x.one_norm();
00525 norm = max(a, norm);
00526 }
00527 return norm;
00528 }
00529
00531 template <typename vt = value_type,
00532 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
00533 typename FieldTraits<vt>::real_type infinity_norm_real() const {
00534 using real_type = typename FieldTraits<vt>::real_type;
00535 using std::max;
00536
00537 real_type norm = 0;
00538 for (auto const &x : *this) {
00539 real_type const a = x.one_norm_real();
00540 norm = max(a, norm);
00541 }
00542 return norm;
00543 }
00544
00546 template <typename vt = value_type,
00547 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
00548 typename FieldTraits<vt>::real_type infinity_norm() const {
00549 using real_type = typename FieldTraits<vt>::real_type;
00550 using std::max;
00551
00552 real_type norm = 0;
00553 real_type isNaN = 1;
00554 for (auto const &x : *this) {
00555 real_type const a = x.one_norm();
00556 norm = max(a, norm);
00557 isNaN += a;
00558 }
00559 isNaN /= isNaN;
00560 return norm * isNaN;
00561 }
00562
00564 template <typename vt = value_type,
00565 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
00566 typename FieldTraits<vt>::real_type infinity_norm_real() const {
00567 using real_type = typename FieldTraits<vt>::real_type;
00568 using std::max;
00569
00570 real_type norm = 0;
00571 real_type isNaN = 1;
00572 for (auto const &x : *this) {
00573 real_type const a = x.one_norm_real();
00574 norm = max(a, norm);
00575 isNaN += a;
00576 }
00577 isNaN /= isNaN;
00578 return norm * isNaN;
00579 }
00580
00581
00582
00587 template <class V>
00588 void solve (V& x, const V& b) const;
00589
00594 void invert();
00595
00597 field_type determinant () const;
00598
00600 template<typename M2>
00601 MAT& leftmultiply (const DenseMatrix<M2>& M)
00602 {
00603 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
00604 DUNE_ASSERT_BOUNDS(M.rows() == rows());
00605 MAT C(asImp());
00606
00607 for (size_type i=0; i<rows(); i++)
00608 for (size_type j=0; j<cols(); j++) {
00609 (*this)[i][j] = 0;
00610 for (size_type k=0; k<rows(); k++)
00611 (*this)[i][j] += M[i][k]*C[k][j];
00612 }
00613
00614 return asImp();
00615 }
00616
00618 template<typename M2>
00619 MAT& rightmultiply (const DenseMatrix<M2>& M)
00620 {
00621 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
00622 DUNE_ASSERT_BOUNDS(M.cols() == cols());
00623 MAT C(asImp());
00624
00625 for (size_type i=0; i<rows(); i++)
00626 for (size_type j=0; j<cols(); j++) {
00627 (*this)[i][j] = 0;
00628 for (size_type k=0; k<cols(); k++)
00629 (*this)[i][j] += C[i][k]*M[k][j];
00630 }
00631 return asImp();
00632 }
00633
00634 #if 0
00636 template<int l>
00637 DenseMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
00638 {
00639 FieldMatrix<K,l,cols> C;
00640
00641 for (size_type i=0; i<l; i++) {
00642 for (size_type j=0; j<cols(); j++) {
00643 C[i][j] = 0;
00644 for (size_type k=0; k<rows(); k++)
00645 C[i][j] += M[i][k]*(*this)[k][j];
00646 }
00647 }
00648 return C;
00649 }
00650
00652 template<int l>
00653 FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
00654 {
00655 FieldMatrix<K,rows,l> C;
00656
00657 for (size_type i=0; i<rows(); i++) {
00658 for (size_type j=0; j<l; j++) {
00659 C[i][j] = 0;
00660 for (size_type k=0; k<cols(); k++)
00661 C[i][j] += (*this)[i][k]*M[k][j];
00662 }
00663 }
00664 return C;
00665 }
00666 #endif
00667
00668
00669
00671 size_type N () const
00672 {
00673 return rows();
00674 }
00675
00677 size_type M () const
00678 {
00679 return cols();
00680 }
00681
00683 size_type rows() const
00684 {
00685 return asImp().mat_rows();
00686 }
00687
00689 size_type cols() const
00690 {
00691 return asImp().mat_cols();
00692 }
00693
00694
00695
00697 bool exists (size_type i, size_type j) const
00698 {
00699 DUNE_UNUSED_PARAMETER(i);
00700 DUNE_UNUSED_PARAMETER(j);
00701 DUNE_ASSERT_BOUNDS(i >= 0 && i < rows());
00702 DUNE_ASSERT_BOUNDS(j >= 0 && j < cols());
00703 return true;
00704 }
00705
00706 private:
00707
00708 #ifndef DOXYGEN
00709 struct ElimPivot
00710 {
00711 ElimPivot(std::vector<size_type> & pivot);
00712
00713 void swap(int i, int j);
00714
00715 template<typename T>
00716 void operator()(const T&, int, int)
00717 {}
00718
00719 std::vector<size_type> & pivot_;
00720 };
00721
00722 template<typename V>
00723 struct Elim
00724 {
00725 Elim(V& rhs);
00726
00727 void swap(int i, int j);
00728
00729 void operator()(const typename V::field_type& factor, int k, int i);
00730
00731 V* rhs_;
00732 };
00733
00734 struct ElimDet
00735 {
00736 ElimDet(field_type& sign) : sign_(sign)
00737 { sign_ = 1; }
00738
00739 void swap(int, int)
00740 { sign_ *= -1; }
00741
00742 void operator()(const field_type&, int, int)
00743 {}
00744
00745 field_type& sign_;
00746 };
00747 #endif // DOXYGEN
00748
00749 template<class Func>
00750 void luDecomposition(DenseMatrix<MAT>& A, Func func) const;
00751 };
00752
00753 #ifndef DOXYGEN
00754 template<typename MAT>
00755 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
00756 : pivot_(pivot)
00757 {
00758 typedef typename std::vector<size_type>::size_type size_type;
00759 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
00760 }
00761
00762 template<typename MAT>
00763 void DenseMatrix<MAT>::ElimPivot::swap(int i, int j)
00764 {
00765 pivot_[i]=j;
00766 }
00767
00768 template<typename MAT>
00769 template<typename V>
00770 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
00771 : rhs_(&rhs)
00772 {}
00773
00774 template<typename MAT>
00775 template<typename V>
00776 void DenseMatrix<MAT>::Elim<V>::swap(int i, int j)
00777 {
00778 std::swap((*rhs_)[i], (*rhs_)[j]);
00779 }
00780
00781 template<typename MAT>
00782 template<typename V>
00783 void DenseMatrix<MAT>::
00784 Elim<V>::operator()(const typename V::field_type& factor, int k, int i)
00785 {
00786 (*rhs_)[k] -= factor*(*rhs_)[i];
00787 }
00788 template<typename MAT>
00789 template<typename Func>
00790 inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func) const
00791 {
00792 typedef typename FieldTraits<value_type>::real_type
00793 real_type;
00794 real_type norm = A.infinity_norm_real();
00795 real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() );
00796 real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() );
00797
00798
00799 for (size_type i=0; i<rows(); i++)
00800 {
00801 typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]);
00802
00803
00804 if (pivmax<pivthres)
00805 {
00806
00807 size_type imax=i;
00808 typename FieldTraits<value_type>::real_type abs(0.0);
00809 for (size_type k=i+1; k<rows(); k++)
00810 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
00811 {
00812 pivmax = abs; imax = k;
00813 }
00814
00815 if (imax!=i) {
00816 for (size_type j=0; j<rows(); j++)
00817 std::swap(A[i][j],A[imax][j]);
00818 func.swap(i, imax);
00819 }
00820 }
00821
00822
00823 if (pivmax<singthres)
00824 DUNE_THROW(FMatrixError,"matrix is singular");
00825
00826
00827 for (size_type k=i+1; k<rows(); k++)
00828 {
00829 field_type factor = A[k][i]/A[i][i];
00830 A[k][i] = factor;
00831 for (size_type j=i+1; j<rows(); j++)
00832 A[k][j] -= factor*A[i][j];
00833 func(factor, k, i);
00834 }
00835 }
00836 }
00837
00838 template<typename MAT>
00839 template <class V>
00840 inline void DenseMatrix<MAT>::solve(V& x, const V& b) const
00841 {
00842
00843 if (rows()!=cols())
00844 DUNE_THROW(FMatrixError, "Can't solve for a " << rows() << "x" << cols() << " matrix!");
00845
00846 if (rows()==1) {
00847
00848 #ifdef DUNE_FMatrix_WITH_CHECKING
00849 if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit())
00850 DUNE_THROW(FMatrixError,"matrix is singular");
00851 #endif
00852 x[0] = b[0]/(*this)[0][0];
00853
00854 }
00855 else if (rows()==2) {
00856
00857 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
00858 #ifdef DUNE_FMatrix_WITH_CHECKING
00859 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
00860 DUNE_THROW(FMatrixError,"matrix is singular");
00861 #endif
00862 detinv = 1.0/detinv;
00863
00864 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
00865 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
00866
00867 }
00868 else if (rows()==3) {
00869
00870 field_type d = determinant();
00871 #ifdef DUNE_FMatrix_WITH_CHECKING
00872 if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit())
00873 DUNE_THROW(FMatrixError,"matrix is singular");
00874 #endif
00875
00876 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
00877 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
00878 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
00879
00880 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
00881 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
00882 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
00883
00884 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
00885 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
00886 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
00887
00888 }
00889 else {
00890
00891 V& rhs = x;
00892 rhs = b;
00893 Elim<V> elim(rhs);
00894 MAT A(asImp());
00895
00896 luDecomposition(A, elim);
00897
00898
00899 for(int i=rows()-1; i>=0; i--) {
00900 for (size_type j=i+1; j<rows(); j++)
00901 rhs[i] -= A[i][j]*x[j];
00902 x[i] = rhs[i]/A[i][i];
00903 }
00904 }
00905 }
00906
00907 template<typename MAT>
00908 inline void DenseMatrix<MAT>::invert()
00909 {
00910
00911 if (rows()!=cols())
00912 DUNE_THROW(FMatrixError, "Can't invert a " << rows() << "x" << cols() << " matrix!");
00913
00914 if (rows()==1) {
00915
00916 #ifdef DUNE_FMatrix_WITH_CHECKING
00917 if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit())
00918 DUNE_THROW(FMatrixError,"matrix is singular");
00919 #endif
00920 (*this)[0][0] = field_type( 1 ) / (*this)[0][0];
00921
00922 }
00923 else if (rows()==2) {
00924
00925 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
00926 #ifdef DUNE_FMatrix_WITH_CHECKING
00927 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
00928 DUNE_THROW(FMatrixError,"matrix is singular");
00929 #endif
00930 detinv = field_type( 1 ) / detinv;
00931
00932 field_type temp=(*this)[0][0];
00933 (*this)[0][0] = (*this)[1][1]*detinv;
00934 (*this)[0][1] = -(*this)[0][1]*detinv;
00935 (*this)[1][0] = -(*this)[1][0]*detinv;
00936 (*this)[1][1] = temp*detinv;
00937
00938 }
00939 else if (rows()==3)
00940 {
00941 using K = field_type;
00942
00943 K t4 = (*this)[0][0] * (*this)[1][1];
00944 K t6 = (*this)[0][0] * (*this)[1][2];
00945 K t8 = (*this)[0][1] * (*this)[1][0];
00946 K t10 = (*this)[0][2] * (*this)[1][0];
00947 K t12 = (*this)[0][1] * (*this)[2][0];
00948 K t14 = (*this)[0][2] * (*this)[2][0];
00949
00950 K det = (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
00951 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
00952 K t17 = K(1.0)/det;
00953
00954 K matrix01 = (*this)[0][1];
00955 K matrix00 = (*this)[0][0];
00956 K matrix10 = (*this)[1][0];
00957 K matrix11 = (*this)[1][1];
00958
00959 (*this)[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1])*t17;
00960 (*this)[0][1] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1])*t17;
00961 (*this)[0][2] = (matrix01 * (*this)[1][2] - (*this)[0][2] * (*this)[1][1])*t17;
00962 (*this)[1][0] = -((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0])*t17;
00963 (*this)[1][1] = (matrix00 * (*this)[2][2] - t14) * t17;
00964 (*this)[1][2] = -(t6-t10) * t17;
00965 (*this)[2][0] = (matrix10 * (*this)[2][1] - matrix11 * (*this)[2][0]) * t17;
00966 (*this)[2][1] = -(matrix00 * (*this)[2][1] - t12) * t17;
00967 (*this)[2][2] = (t4-t8) * t17;
00968 }
00969 else {
00970
00971 MAT A(asImp());
00972 std::vector<size_type> pivot(rows());
00973 luDecomposition(A, ElimPivot(pivot));
00974 DenseMatrix<MAT>& L=A;
00975 DenseMatrix<MAT>& U=A;
00976
00977
00978 *this=field_type();
00979
00980 for(size_type i=0; i<rows(); ++i)
00981 (*this)[i][i]=1;
00982
00983
00984 for (size_type i=0; i<rows(); i++)
00985 for (size_type j=0; j<i; j++)
00986 for (size_type k=0; k<rows(); k++)
00987 (*this)[i][k] -= L[i][j]*(*this)[j][k];
00988
00989
00990 for (size_type i=rows(); i>0;) {
00991 --i;
00992 for (size_type k=0; k<rows(); k++) {
00993 for (size_type j=i+1; j<rows(); j++)
00994 (*this)[i][k] -= U[i][j]*(*this)[j][k];
00995 (*this)[i][k] /= U[i][i];
00996 }
00997 }
00998
00999 for(size_type i=rows(); i>0; ) {
01000 --i;
01001 if(i!=pivot[i])
01002 for(size_type j=0; j<rows(); ++j)
01003 std::swap((*this)[j][pivot[i]], (*this)[j][i]);
01004 }
01005 }
01006 }
01007
01008
01009 template<typename MAT>
01010 inline typename DenseMatrix<MAT>::field_type
01011 DenseMatrix<MAT>::determinant() const
01012 {
01013
01014 if (rows()!=cols())
01015 DUNE_THROW(FMatrixError, "There is no determinant for a " << rows() << "x" << cols() << " matrix!");
01016
01017 if (rows()==1)
01018 return (*this)[0][0];
01019
01020 if (rows()==2)
01021 return (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
01022
01023 if (rows()==3) {
01024
01025 field_type t4 = (*this)[0][0] * (*this)[1][1];
01026 field_type t6 = (*this)[0][0] * (*this)[1][2];
01027 field_type t8 = (*this)[0][1] * (*this)[1][0];
01028 field_type t10 = (*this)[0][2] * (*this)[1][0];
01029 field_type t12 = (*this)[0][1] * (*this)[2][0];
01030 field_type t14 = (*this)[0][2] * (*this)[2][0];
01031
01032 return (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
01033 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
01034
01035 }
01036
01037 MAT A(asImp());
01038 field_type det;
01039 try
01040 {
01041 luDecomposition(A, ElimDet(det));
01042 }
01043 catch (FMatrixError&)
01044 {
01045 return 0;
01046 }
01047 for (size_type i = 0; i < rows(); ++i)
01048 det *= A[i][i];
01049 return det;
01050 }
01051
01052 #endif // DOXYGEN
01053
01054 namespace DenseMatrixHelp {
01055
01057 template <typename MAT, typename V1, typename V2>
01058 static inline void multAssign(const DenseMatrix<MAT> &matrix, const DenseVector<V1> & x, DenseVector<V2> & ret)
01059 {
01060 DUNE_ASSERT_BOUNDS(x.size() == matrix.cols());
01061 DUNE_ASSERT_BOUNDS(ret.size() == matrix.rows());
01062 typedef typename DenseMatrix<MAT>::size_type size_type;
01063
01064 for(size_type i=0; i<matrix.rows(); ++i)
01065 {
01066 ret[i] = 0.0;
01067 for(size_type j=0; j<matrix.cols(); ++j)
01068 {
01069 ret[i] += matrix[i][j]*x[j];
01070 }
01071 }
01072 }
01073
01074 #if 0
01076 template <typename K, int rows, int cols>
01077 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
01078 {
01079 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
01080
01081 for(size_type i=0; i<cols(); ++i)
01082 {
01083 ret[i] = 0.0;
01084 for(size_type j=0; j<rows(); ++j)
01085 ret[i] += matrix[j][i]*x[j];
01086 }
01087 }
01088
01090 template <typename K, int rows, int cols>
01091 static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
01092 {
01093 FieldVector<K,rows> ret;
01094 multAssign(matrix,x,ret);
01095 return ret;
01096 }
01097
01099 template <typename K, int rows, int cols>
01100 static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x)
01101 {
01102 FieldVector<K,cols> ret;
01103 multAssignTransposed( matrix, x, ret );
01104 return ret;
01105 }
01106 #endif
01107
01108 }
01109
01111 template<typename MAT>
01112 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
01113 {
01114 for (typename DenseMatrix<MAT>::size_type i=0; i<a.rows(); i++)
01115 s << a[i] << std::endl;
01116 return s;
01117 }
01118
01121 }
01122
01123 #endif