dune-common 3.0-git
densematrix.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_DENSEMATRIX_HH
4#define DUNE_DENSEMATRIX_HH
5
6#include <cmath>
7#include <cstddef>
8#include <iostream>
9#include <type_traits>
10#include <vector>
11
17#include <dune/common/math.hh>
18#include <dune/common/unused.hh>
19
20
21namespace Dune
22{
23
24 template<typename M> class DenseMatrix;
25
26 template<typename M>
32
33 /*
34 work around a problem of FieldMatrix/FieldVector,
35 there is no unique way to obtain the size of a class
36 */
37 template<class K, int N, int M> class FieldMatrix;
38 template<class K, int N> class FieldVector;
39 namespace {
40 template<class V>
41 struct VectorSize
42 {
43 static typename V::size_type size(const V & v) { return v.size(); }
44 };
45
46 template<class K, int N>
47 struct VectorSize< const FieldVector<K,N> >
48 {
49 typedef FieldVector<K,N> V;
50 static typename V::size_type size(const V & v) { return N; }
51 };
52 }
53
72 template< class DenseMatrix, class RHS >
74
75#ifndef DOXYGEN
76 namespace
77 {
78 template< class DenseMatrix, class RHS,
79 bool primitive = std::is_convertible< RHS, typename DenseMatrix::field_type >::value >
80 class DenseMatrixAssignerImplementation;
81
82 template< class DenseMatrix, class RHS >
83 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, true >
84 {
85 public:
86 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
87 {
88 typedef typename DenseMatrix::field_type field_type;
89 std::fill( denseMatrix.begin(), denseMatrix.end(), static_cast< field_type >( rhs ) );
90 }
91 };
92
93 template< class DenseMatrix, class RHS >
94 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, false >
95 {
96 public:
97 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
98 {
99 static_assert( (std::is_convertible< const RHS, const DenseMatrix >::value),
100 "No template specialization of DenseMatrixAssigner found" );
101 denseMatrix = static_cast< const DenseMatrix & >( rhs );
102 }
103 };
104 }
105
106
107
108 template< class DenseMatrix, class RHS >
109 struct DenseMatrixAssigner
110 {
111 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
112 {
113 DenseMatrixAssignerImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs );
114 }
115 };
116#endif // #ifndef DOXYGEN
117
118
119
121 class FMatrixError : public MathError {};
122
133 template<typename MAT>
135 {
136 typedef DenseMatVecTraits<MAT> Traits;
137
138 // Curiously recurring template pattern
139 MAT & asImp() { return static_cast<MAT&>(*this); }
140 const MAT & asImp() const { return static_cast<const MAT&>(*this); }
141
142 public:
143 //===== type definitions and constants
144
146 typedef typename Traits::derived_type derived_type;
147
149 typedef typename Traits::value_type value_type;
150
152 typedef typename Traits::value_type field_type;
153
155 typedef typename Traits::value_type block_type;
156
158 typedef typename Traits::size_type size_type;
159
161 typedef typename Traits::row_type row_type;
162
164 typedef typename Traits::row_reference row_reference;
165
167 typedef typename Traits::const_row_reference const_row_reference;
168
170 enum {
172 blocklevel = 1
173 };
174
175 //===== access to components
176
179 {
180 return asImp().mat_access(i);
181 }
182
184 {
185 return asImp().mat_access(i);
186 }
187
190 {
191 return rows();
192 }
193
194 //===== iterator interface to rows of the matrix
202 typedef typename std::remove_reference<row_reference>::type::Iterator ColIterator;
203
206 {
207 return Iterator(*this,0);
208 }
209
212 {
213 return Iterator(*this,rows());
214 }
215
219 {
220 return Iterator(*this,rows()-1);
221 }
222
226 {
227 return Iterator(*this,-1);
228 }
229
237 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator ConstColIterator;
238
241 {
242 return ConstIterator(*this,0);
243 }
244
247 {
248 return ConstIterator(*this,rows());
249 }
250
254 {
255 return ConstIterator(*this,rows()-1);
256 }
257
261 {
262 return ConstIterator(*this,-1);
263 }
264
265 //===== assignment
266
267 template< class RHS >
268 DenseMatrix &operator= ( const RHS &rhs )
269 {
271 return *this;
272 }
273
274 //===== vector space arithmetic
275
277 template <class Other>
279 {
280 DUNE_ASSERT_BOUNDS(rows() == y.rows());
281 for (size_type i=0; i<rows(); i++)
282 (*this)[i] += y[i];
283 return *this;
284 }
285
287 template <class Other>
289 {
290 DUNE_ASSERT_BOUNDS(rows() == y.rows());
291 for (size_type i=0; i<rows(); i++)
292 (*this)[i] -= y[i];
293 return *this;
294 }
295
298 {
299 for (size_type i=0; i<rows(); i++)
300 (*this)[i] *= k;
301 return *this;
302 }
303
306 {
307 for (size_type i=0; i<rows(); i++)
308 (*this)[i] /= k;
309 return *this;
310 }
311
313 template <class Other>
315 {
316 DUNE_ASSERT_BOUNDS(rows() == y.rows());
317 for( size_type i = 0; i < rows(); ++i )
318 (*this)[ i ].axpy( k, y[ i ] );
319 return *this;
320 }
321
323 template <class Other>
324 bool operator== (const DenseMatrix<Other>& y) const
325 {
326 DUNE_ASSERT_BOUNDS(rows() == y.rows());
327 for (size_type i=0; i<rows(); i++)
328 if ((*this)[i]!=y[i])
329 return false;
330 return true;
331 }
333 template <class Other>
334 bool operator!= (const DenseMatrix<Other>& y) const
335 {
336 return !operator==(y);
337 }
338
339
340 //===== linear maps
341
343 template<class X, class Y>
344 void mv (const X& x, Y& y) const
345 {
346 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
347 DUNE_ASSERT_BOUNDS(x.N() == M());
348 DUNE_ASSERT_BOUNDS(y.N() == N());
349
350 using field_type = typename FieldTraits<Y>::field_type;
351 for (size_type i=0; i<rows(); ++i)
352 {
353 y[i] = field_type(0);
354 for (size_type j=0; j<cols(); j++)
355 y[i] += (*this)[i][j] * x[j];
356 }
357 }
358
360 template< class X, class Y >
361 void mtv ( const X &x, Y &y ) const
362 {
363 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
364 DUNE_ASSERT_BOUNDS(x.N() == N());
365 DUNE_ASSERT_BOUNDS(y.N() == M());
366
367 using field_type = typename FieldTraits<Y>::field_type;
368 for(size_type i = 0; i < cols(); ++i)
369 {
370 y[i] = field_type(0);
371 for(size_type j = 0; j < rows(); ++j)
372 y[i] += (*this)[j][i] * x[j];
373 }
374 }
375
377 template<class X, class Y>
378 void umv (const X& x, Y& y) const
379 {
380 DUNE_ASSERT_BOUNDS(x.N() == M());
381 DUNE_ASSERT_BOUNDS(y.N() == N());
382 for (size_type i=0; i<rows(); i++)
383 for (size_type j=0; j<cols(); j++)
384 y[i] += (*this)[i][j] * x[j];
385 }
386
388 template<class X, class Y>
389 void umtv (const X& x, Y& y) const
390 {
391 DUNE_ASSERT_BOUNDS(x.N() == N());
392 DUNE_ASSERT_BOUNDS(y.N() == M());
393 for (size_type i=0; i<rows(); i++)
394 for (size_type j=0; j<cols(); j++)
395 y[j] += (*this)[i][j]*x[i];
396 }
397
399 template<class X, class Y>
400 void umhv (const X& x, Y& y) const
401 {
402 DUNE_ASSERT_BOUNDS(x.N() == N());
403 DUNE_ASSERT_BOUNDS(y.N() == M());
404 for (size_type i=0; i<rows(); i++)
405 for (size_type j=0; j<cols(); j++)
406 y[j] += conjugateComplex((*this)[i][j])*x[i];
407 }
408
410 template<class X, class Y>
411 void mmv (const X& x, Y& y) const
412 {
413 DUNE_ASSERT_BOUNDS(x.N() == M());
414 DUNE_ASSERT_BOUNDS(y.N() == N());
415 for (size_type i=0; i<rows(); i++)
416 for (size_type j=0; j<cols(); j++)
417 y[i] -= (*this)[i][j] * x[j];
418 }
419
421 template<class X, class Y>
422 void mmtv (const X& x, Y& y) const
423 {
424 DUNE_ASSERT_BOUNDS(x.N() == N());
425 DUNE_ASSERT_BOUNDS(y.N() == M());
426 for (size_type i=0; i<rows(); i++)
427 for (size_type j=0; j<cols(); j++)
428 y[j] -= (*this)[i][j]*x[i];
429 }
430
432 template<class X, class Y>
433 void mmhv (const X& x, Y& y) const
434 {
435 DUNE_ASSERT_BOUNDS(x.N() == N());
436 DUNE_ASSERT_BOUNDS(y.N() == M());
437 for (size_type i=0; i<rows(); i++)
438 for (size_type j=0; j<cols(); j++)
439 y[j] -= conjugateComplex((*this)[i][j])*x[i];
440 }
441
443 template<class X, class Y>
444 void usmv (const typename FieldTraits<Y>::field_type & alpha,
445 const X& x, Y& y) const
446 {
447 DUNE_ASSERT_BOUNDS(x.N() == M());
448 DUNE_ASSERT_BOUNDS(y.N() == N());
449 for (size_type i=0; i<rows(); i++)
450 for (size_type j=0; j<cols(); j++)
451 y[i] += alpha * (*this)[i][j] * x[j];
452 }
453
455 template<class X, class Y>
456 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
457 const X& x, Y& y) const
458 {
459 DUNE_ASSERT_BOUNDS(x.N() == N());
460 DUNE_ASSERT_BOUNDS(y.N() == M());
461 for (size_type i=0; i<rows(); i++)
462 for (size_type j=0; j<cols(); j++)
463 y[j] += alpha*(*this)[i][j]*x[i];
464 }
465
467 template<class X, class Y>
468 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
469 const X& x, Y& y) const
470 {
471 DUNE_ASSERT_BOUNDS(x.N() == N());
472 DUNE_ASSERT_BOUNDS(y.N() == M());
473 for (size_type i=0; i<rows(); i++)
474 for (size_type j=0; j<cols(); j++)
475 y[j] += alpha*conjugateComplex((*this)[i][j])*x[i];
476 }
477
478 //===== norms
479
482 {
483 typename FieldTraits<value_type>::real_type sum=(0.0);
484 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
485 return fvmeta::sqrt(sum);
486 }
487
490 {
491 typename FieldTraits<value_type>::real_type sum=(0.0);
492 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
493 return sum;
494 }
495
497 template <typename vt = value_type,
498 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
500 using real_type = typename FieldTraits<vt>::real_type;
501 using std::max;
502
503 real_type norm = 0;
504 for (auto const &x : *this) {
505 real_type const a = x.one_norm();
506 norm = max(a, norm);
507 }
508 return norm;
509 }
510
512 template <typename vt = value_type,
513 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
515 using real_type = typename FieldTraits<vt>::real_type;
516 using std::max;
517
518 real_type norm = 0;
519 for (auto const &x : *this) {
520 real_type const a = x.one_norm_real();
521 norm = max(a, norm);
522 }
523 return norm;
524 }
525
527 template <typename vt = value_type,
528 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
530 using real_type = typename FieldTraits<vt>::real_type;
531 using std::max;
532
533 real_type norm = 0;
534 real_type isNaN = 1;
535 for (auto const &x : *this) {
536 real_type const a = x.one_norm();
537 norm = max(a, norm);
538 isNaN += a;
539 }
540 isNaN /= isNaN;
541 return norm * isNaN;
542 }
543
545 template <typename vt = value_type,
546 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
548 using real_type = typename FieldTraits<vt>::real_type;
549 using std::max;
550
551 real_type norm = 0;
552 real_type isNaN = 1;
553 for (auto const &x : *this) {
554 real_type const a = x.one_norm_real();
555 norm = max(a, norm);
556 isNaN += a;
557 }
558 isNaN /= isNaN;
559 return norm * isNaN;
560 }
561
562 //===== solve
563
568 template <class V>
569 void solve (V& x, const V& b) const;
570
575 void invert();
576
579
581 template<typename M2>
583 {
584 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
585 DUNE_ASSERT_BOUNDS(M.rows() == rows());
586 MAT C(asImp());
587
588 for (size_type i=0; i<rows(); i++)
589 for (size_type j=0; j<cols(); j++) {
590 (*this)[i][j] = 0;
591 for (size_type k=0; k<rows(); k++)
592 (*this)[i][j] += M[i][k]*C[k][j];
593 }
594
595 return asImp();
596 }
597
599 template<typename M2>
601 {
602 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
603 DUNE_ASSERT_BOUNDS(M.cols() == cols());
604 MAT C(asImp());
605
606 for (size_type i=0; i<rows(); i++)
607 for (size_type j=0; j<cols(); j++) {
608 (*this)[i][j] = 0;
609 for (size_type k=0; k<cols(); k++)
610 (*this)[i][j] += C[i][k]*M[k][j];
611 }
612 return asImp();
613 }
614
615#if 0
617 template<int l>
618 DenseMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
619 {
621
622 for (size_type i=0; i<l; i++) {
623 for (size_type j=0; j<cols(); j++) {
624 C[i][j] = 0;
625 for (size_type k=0; k<rows(); k++)
626 C[i][j] += M[i][k]*(*this)[k][j];
627 }
628 }
629 return C;
630 }
631
633 template<int l>
634 FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
635 {
636 FieldMatrix<K,rows,l> C;
637
638 for (size_type i=0; i<rows(); i++) {
639 for (size_type j=0; j<l; j++) {
640 C[i][j] = 0;
641 for (size_type k=0; k<cols(); k++)
642 C[i][j] += (*this)[i][k]*M[k][j];
643 }
644 }
645 return C;
646 }
647#endif
648
649 //===== sizes
650
652 size_type N () const
653 {
654 return rows();
655 }
656
658 size_type M () const
659 {
660 return cols();
661 }
662
665 {
666 return asImp().mat_rows();
667 }
668
671 {
672 return asImp().mat_cols();
673 }
674
675 //===== query
676
678 bool exists (size_type i, size_type j) const
679 {
682 DUNE_ASSERT_BOUNDS(i >= 0 && i < rows());
683 DUNE_ASSERT_BOUNDS(j >= 0 && j < cols());
684 return true;
685 }
686
687 private:
688
689#ifndef DOXYGEN
690 struct ElimPivot
691 {
692 ElimPivot(std::vector<size_type> & pivot);
693
694 void swap(int i, int j);
695
696 template<typename T>
697 void operator()(const T&, int, int)
698 {}
699
700 std::vector<size_type> & pivot_;
701 };
702
703 template<typename V>
704 struct Elim
705 {
706 Elim(V& rhs);
707
708 void swap(int i, int j);
709
710 void operator()(const typename V::field_type& factor, int k, int i);
711
712 V* rhs_;
713 };
714
715 struct ElimDet
716 {
717 ElimDet(field_type& sign) : sign_(sign)
718 { sign_ = 1; }
719
720 void swap(int, int)
721 { sign_ *= -1; }
722
723 void operator()(const field_type&, int, int)
724 {}
725
726 field_type& sign_;
727 };
728#endif // DOXYGEN
729
730 template<class Func>
731 void luDecomposition(DenseMatrix<MAT>& A, Func func) const;
732 };
733
734#ifndef DOXYGEN
735 template<typename MAT>
736 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
737 : pivot_(pivot)
738 {
739 typedef typename std::vector<size_type>::size_type size_type;
740 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
741 }
742
743 template<typename MAT>
744 void DenseMatrix<MAT>::ElimPivot::swap(int i, int j)
745 {
746 pivot_[i]=j;
747 }
748
749 template<typename MAT>
750 template<typename V>
751 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
752 : rhs_(&rhs)
753 {}
754
755 template<typename MAT>
756 template<typename V>
757 void DenseMatrix<MAT>::Elim<V>::swap(int i, int j)
758 {
759 std::swap((*rhs_)[i], (*rhs_)[j]);
760 }
761
762 template<typename MAT>
763 template<typename V>
764 void DenseMatrix<MAT>::
765 Elim<V>::operator()(const typename V::field_type& factor, int k, int i)
766 {
767 (*rhs_)[k] -= factor*(*rhs_)[i];
768 }
769 template<typename MAT>
770 template<typename Func>
771 inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func) const
772 {
774 real_type;
775 real_type norm = A.infinity_norm_real(); // for relative thresholds
778
779 // LU decomposition of A in A
780 for (size_type i=0; i<rows(); i++) // loop over all rows
781 {
782 typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]);
783
784 // pivoting ?
785 if (pivmax<pivthres)
786 {
787 // compute maximum of column
788 size_type imax=i;
789 typename FieldTraits<value_type>::real_type abs(0.0);
790 for (size_type k=i+1; k<rows(); k++)
791 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
792 {
793 pivmax = abs; imax = k;
794 }
795 // swap rows
796 if (imax!=i) {
797 for (size_type j=0; j<rows(); j++)
798 std::swap(A[i][j],A[imax][j]);
799 func.swap(i, imax); // swap the pivot or rhs
800 }
801 }
802
803 // singular ?
804 if (pivmax<singthres)
805 DUNE_THROW(FMatrixError,"matrix is singular");
806
807 // eliminate
808 for (size_type k=i+1; k<rows(); k++)
809 {
810 field_type factor = A[k][i]/A[i][i];
811 A[k][i] = factor;
812 for (size_type j=i+1; j<rows(); j++)
813 A[k][j] -= factor*A[i][j];
814 func(factor, k, i);
815 }
816 }
817 }
818
819 template<typename MAT>
820 template <class V>
821 inline void DenseMatrix<MAT>::solve(V& x, const V& b) const
822 {
823 // never mind those ifs, because they get optimized away
824 if (rows()!=cols())
825 DUNE_THROW(FMatrixError, "Can't solve for a " << rows() << "x" << cols() << " matrix!");
826
827 if (rows()==1) {
828
829#ifdef DUNE_FMatrix_WITH_CHECKING
830 if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit())
831 DUNE_THROW(FMatrixError,"matrix is singular");
832#endif
833 x[0] = b[0]/(*this)[0][0];
834
835 }
836 else if (rows()==2) {
837
838 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
839#ifdef DUNE_FMatrix_WITH_CHECKING
840 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
841 DUNE_THROW(FMatrixError,"matrix is singular");
842#endif
843 detinv = 1.0/detinv;
844
845 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
846 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
847
848 }
849 else if (rows()==3) {
850
851 field_type d = determinant();
852#ifdef DUNE_FMatrix_WITH_CHECKING
853 if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit())
854 DUNE_THROW(FMatrixError,"matrix is singular");
855#endif
856
857 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
858 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
859 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
860
861 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
862 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
863 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
864
865 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
866 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
867 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
868
869 }
870 else {
871
872 V& rhs = x; // use x to store rhs
873 rhs = b; // copy data
874 Elim<V> elim(rhs);
875 MAT A(asImp());
876
877 luDecomposition(A, elim);
878
879 // backsolve
880 for(int i=rows()-1; i>=0; i--) {
881 for (size_type j=i+1; j<rows(); j++)
882 rhs[i] -= A[i][j]*x[j];
883 x[i] = rhs[i]/A[i][i];
884 }
885 }
886 }
887
888 template<typename MAT>
889 inline void DenseMatrix<MAT>::invert()
890 {
891 // never mind those ifs, because they get optimized away
892 if (rows()!=cols())
893 DUNE_THROW(FMatrixError, "Can't invert a " << rows() << "x" << cols() << " matrix!");
894
895 if (rows()==1) {
896
897#ifdef DUNE_FMatrix_WITH_CHECKING
898 if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit())
899 DUNE_THROW(FMatrixError,"matrix is singular");
900#endif
901 (*this)[0][0] = field_type( 1 ) / (*this)[0][0];
902
903 }
904 else if (rows()==2) {
905
906 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
907#ifdef DUNE_FMatrix_WITH_CHECKING
908 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
909 DUNE_THROW(FMatrixError,"matrix is singular");
910#endif
911 detinv = field_type( 1 ) / detinv;
912
913 field_type temp=(*this)[0][0];
914 (*this)[0][0] = (*this)[1][1]*detinv;
915 (*this)[0][1] = -(*this)[0][1]*detinv;
916 (*this)[1][0] = -(*this)[1][0]*detinv;
917 (*this)[1][1] = temp*detinv;
918
919 }
920 else {
921
922 MAT A(asImp());
923 std::vector<size_type> pivot(rows());
924 luDecomposition(A, ElimPivot(pivot));
925 DenseMatrix<MAT>& L=A;
926 DenseMatrix<MAT>& U=A;
927
928 // initialize inverse
929 *this=field_type();
930
931 for(size_type i=0; i<rows(); ++i)
932 (*this)[i][i]=1;
933
934 // L Y = I; multiple right hand sides
935 for (size_type i=0; i<rows(); i++)
936 for (size_type j=0; j<i; j++)
937 for (size_type k=0; k<rows(); k++)
938 (*this)[i][k] -= L[i][j]*(*this)[j][k];
939
940 // U A^{-1} = Y
941 for (size_type i=rows(); i>0;) {
942 --i;
943 for (size_type k=0; k<rows(); k++) {
944 for (size_type j=i+1; j<rows(); j++)
945 (*this)[i][k] -= U[i][j]*(*this)[j][k];
946 (*this)[i][k] /= U[i][i];
947 }
948 }
949
950 for(size_type i=rows(); i>0; ) {
951 --i;
952 if(i!=pivot[i])
953 for(size_type j=0; j<rows(); ++j)
954 std::swap((*this)[j][pivot[i]], (*this)[j][i]);
955 }
956 }
957 }
958
959 // implementation of the determinant
960 template<typename MAT>
961 inline typename DenseMatrix<MAT>::field_type
963 {
964 // never mind those ifs, because they get optimized away
965 if (rows()!=cols())
966 DUNE_THROW(FMatrixError, "There is no determinant for a " << rows() << "x" << cols() << " matrix!");
967
968 if (rows()==1)
969 return (*this)[0][0];
970
971 if (rows()==2)
972 return (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
973
974 if (rows()==3) {
975 // code generated by maple
976 field_type t4 = (*this)[0][0] * (*this)[1][1];
977 field_type t6 = (*this)[0][0] * (*this)[1][2];
978 field_type t8 = (*this)[0][1] * (*this)[1][0];
979 field_type t10 = (*this)[0][2] * (*this)[1][0];
980 field_type t12 = (*this)[0][1] * (*this)[2][0];
981 field_type t14 = (*this)[0][2] * (*this)[2][0];
982
983 return (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
984 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
985
986 }
987
988 MAT A(asImp());
989 field_type det;
990 try
991 {
992 luDecomposition(A, ElimDet(det));
993 }
994 catch (FMatrixError&)
995 {
996 return 0;
997 }
998 for (size_type i = 0; i < rows(); ++i)
999 det *= A[i][i];
1000 return det;
1001 }
1002
1003#endif // DOXYGEN
1004
1005 namespace DenseMatrixHelp {
1006#if 0
1008 template <typename K>
1009 static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
1010 {
1011 inverse[0][0] = 1.0/matrix[0][0];
1012 return matrix[0][0];
1013 }
1014
1016 template <typename K>
1017 static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
1018 {
1019 return invertMatrix(matrix,inverse);
1020 }
1021
1022
1024 template <typename K>
1025 static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
1026 {
1027 // code generated by maple
1028 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1029 field_type det_1 = 1.0/det;
1030 inverse[0][0] = matrix[1][1] * det_1;
1031 inverse[0][1] = - matrix[0][1] * det_1;
1032 inverse[1][0] = - matrix[1][0] * det_1;
1033 inverse[1][1] = matrix[0][0] * det_1;
1034 return det;
1035 }
1036
1039 template <typename K>
1040 static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
1041 {
1042 // code generated by maple
1043 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1044 field_type det_1 = 1.0/det;
1045 inverse[0][0] = matrix[1][1] * det_1;
1046 inverse[1][0] = - matrix[0][1] * det_1;
1047 inverse[0][1] = - matrix[1][0] * det_1;
1048 inverse[1][1] = matrix[0][0] * det_1;
1049 return det;
1050 }
1051
1053 template <typename K>
1054 static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
1055 {
1056 // code generated by maple
1057 field_type t4 = matrix[0][0] * matrix[1][1];
1058 field_type t6 = matrix[0][0] * matrix[1][2];
1059 field_type t8 = matrix[0][1] * matrix[1][0];
1060 field_type t10 = matrix[0][2] * matrix[1][0];
1061 field_type t12 = matrix[0][1] * matrix[2][0];
1062 field_type t14 = matrix[0][2] * matrix[2][0];
1063
1064 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1065 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1066 field_type t17 = 1.0/det;
1067
1068 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1069 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1070 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1071 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1072 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1073 inverse[1][2] = -(t6-t10) * t17;
1074 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1075 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1076 inverse[2][2] = (t4-t8) * t17;
1077
1078 return det;
1079 }
1080
1082 template <typename K>
1083 static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
1084 {
1085 // code generated by maple
1086 field_type t4 = matrix[0][0] * matrix[1][1];
1087 field_type t6 = matrix[0][0] * matrix[1][2];
1088 field_type t8 = matrix[0][1] * matrix[1][0];
1089 field_type t10 = matrix[0][2] * matrix[1][0];
1090 field_type t12 = matrix[0][1] * matrix[2][0];
1091 field_type t14 = matrix[0][2] * matrix[2][0];
1092
1093 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1094 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1095 field_type t17 = 1.0/det;
1096
1097 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1098 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1099 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1100 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1101 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1102 inverse[2][1] = -(t6-t10) * t17;
1103 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1104 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1105 inverse[2][2] = (t4-t8) * t17;
1106
1107 return det;
1108 }
1109
1111 template< class K, int m, int n, int p >
1112 static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
1113 const FieldMatrix< K, n, p > &B,
1115 {
1116 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
1117
1118 for( size_type i = 0; i < m; ++i )
1119 {
1120 for( size_type j = 0; j < p; ++j )
1121 {
1122 ret[ i ][ j ] = K( 0 );
1123 for( size_type k = 0; k < n; ++k )
1124 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
1125 }
1126 }
1127 }
1128
1130 template <typename K, int rows, int cols>
1131 static inline void multTransposedMatrix(const FieldMatrix<K,rows,cols> &matrix, FieldMatrix<K,cols,cols>& ret)
1132 {
1133 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
1134
1135 for(size_type i=0; i<cols(); i++)
1136 for(size_type j=0; j<cols(); j++)
1137 {
1138 ret[i][j]=0.0;
1139 for(size_type k=0; k<rows(); k++)
1140 ret[i][j]+=matrix[k][i]*matrix[k][j];
1141 }
1142 }
1143#endif
1144
1146 template <typename MAT, typename V1, typename V2>
1147 static inline void multAssign(const DenseMatrix<MAT> &matrix, const DenseVector<V1> & x, DenseVector<V2> & ret)
1148 {
1149 DUNE_ASSERT_BOUNDS(x.size() == matrix.cols());
1150 DUNE_ASSERT_BOUNDS(ret.size() == matrix.rows());
1151 typedef typename DenseMatrix<MAT>::size_type size_type;
1152
1153 for(size_type i=0; i<matrix.rows(); ++i)
1154 {
1155 ret[i] = 0.0;
1156 for(size_type j=0; j<matrix.cols(); ++j)
1157 {
1158 ret[i] += matrix[i][j]*x[j];
1159 }
1160 }
1161 }
1162
1163#if 0
1165 template <typename K, int rows, int cols>
1166 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
1167 {
1168 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
1169
1170 for(size_type i=0; i<cols(); ++i)
1171 {
1172 ret[i] = 0.0;
1173 for(size_type j=0; j<rows(); ++j)
1174 ret[i] += matrix[j][i]*x[j];
1175 }
1176 }
1177
1179 template <typename K, int rows, int cols>
1180 static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
1181 {
1182 FieldVector<K,rows> ret;
1183 multAssign(matrix,x,ret);
1184 return ret;
1185 }
1186
1188 template <typename K, int rows, int cols>
1189 static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x)
1190 {
1191 FieldVector<K,cols> ret;
1192 multAssignTransposed( matrix, x, ret );
1193 return ret;
1194 }
1195#endif
1196
1197 } // end namespace DenseMatrixHelp
1198
1200 template<typename MAT>
1201 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1202 {
1203 for (typename DenseMatrix<MAT>::size_type i=0; i<a.rows(); i++)
1204 s << a[i] << std::endl;
1205 return s;
1206 }
1207
1210} // end namespace Dune
1211
1212#endif
#define DUNE_ASSERT_BOUNDS(cond)
Definition boundschecking.hh:20
Implements a vector constructed from a given type representing a field and a compile-time given size.
Various precision settings for calculations with FieldMatrix and FieldVector.
Some useful basic math stuff.
A few common exception classes.
A free function to provide the demangled class name of a given object or type as a string.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition unused.hh:18
std::ostream & operator<<(std::ostream &s, const std::array< T, N > &e)
Output operator for array.
Definition array.hh:28
#define DUNE_THROW(E, m)
Definition exceptions.hh:216
Dune namespace.
Definition alignment.hh:11
int sign(const T &val)
Return the sign of the value.
Definition math.hh:119
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition math.hh:103
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition densematrix.hh:1147
constexpr auto size(const Dune::FieldVector< T, i > *, const PriorityTag< 5 > &) -> decltype(std::integral_constant< std::size_t, i >())
Definition hybridutilities.hh:22
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition apply.hh:54
A dense n x m matrix.
Definition densematrix.hh:135
DenseMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition densematrix.hh:305
ConstIterator const_iterator
typedef for stl compliant access
Definition densematrix.hh:233
bool operator!=(const DenseMatrix< Other > &y) const
Binary matrix incomparison.
Definition densematrix.hh:334
void mv(const X &x, Y &y) const
y = A x
Definition densematrix.hh:344
Traits::value_type field_type
export the type representing the field
Definition densematrix.hh:152
void solve(V &x, const V &b) const
Solve system A x = b.
size_type cols() const
number of columns
Definition densematrix.hh:670
ConstIterator beforeEnd() const
Definition densematrix.hh:253
bool operator==(const DenseMatrix< Other > &y) const
Binary matrix comparison.
Definition densematrix.hh:324
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition densematrix.hh:422
DenseMatrix & operator=(const RHS &rhs)
Definition densematrix.hh:268
DenseMatrix & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition densematrix.hh:278
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition densematrix.hh:514
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition densematrix.hh:237
ConstIterator beforeBegin() const
Definition densematrix.hh:260
Traits::value_type block_type
export the type representing the components
Definition densematrix.hh:155
void mtv(const X &x, Y &y) const
y = A^T x
Definition densematrix.hh:361
DenseMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition densematrix.hh:297
DenseMatrix & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition densematrix.hh:288
size_type size() const
size method (number of rows)
Definition densematrix.hh:189
size_type M() const
number of columns
Definition densematrix.hh:658
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition densematrix.hh:600
Iterator end()
end iterator
Definition densematrix.hh:211
Iterator beforeBegin()
Definition densematrix.hh:225
Iterator beforeEnd()
Definition densematrix.hh:218
Traits::value_type value_type
export the type representing the field
Definition densematrix.hh:149
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition densematrix.hh:444
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition densematrix.hh:468
size_type rows() const
number of rows
Definition densematrix.hh:664
void mmv(const X &x, Y &y) const
y -= A x
Definition densematrix.hh:411
void invert()
Compute inverse.
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition densematrix.hh:582
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition densematrix.hh:456
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition densematrix.hh:433
Traits::derived_type derived_type
type of derived matrix class
Definition densematrix.hh:146
row_reference operator[](size_type i)
random access
Definition densematrix.hh:178
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition densematrix.hh:678
DenseMatrix & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition densematrix.hh:314
Iterator RowIterator
rename the iterators for easier access
Definition densematrix.hh:200
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition densematrix.hh:481
void umv(const X &x, Y &y) const
y += A x
Definition densematrix.hh:378
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition densematrix.hh:231
size_type N() const
number of rows
Definition densematrix.hh:652
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition densematrix.hh:499
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition densematrix.hh:161
Traits::size_type size_type
The type used for the index access and size operation.
Definition densematrix.hh:158
Traits::const_row_reference const_row_reference
The type used to represent a reference to a constant row (usually const row_type &)
Definition densematrix.hh:167
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition densematrix.hh:489
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition densematrix.hh:202
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition densematrix.hh:164
Iterator iterator
typedef for stl compliant access
Definition densematrix.hh:198
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition densematrix.hh:235
field_type determinant() const
calculates the determinant of this matrix
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition densematrix.hh:196
void umtv(const X &x, Y &y) const
y += A^T x
Definition densematrix.hh:389
ConstIterator begin() const
begin iterator
Definition densematrix.hh:240
@ blocklevel
The number of block levels we contain. This is 1.
Definition densematrix.hh:172
Iterator begin()
begin iterator
Definition densematrix.hh:205
void umhv(const X &x, Y &y) const
y += A^H x
Definition densematrix.hh:400
ConstIterator end() const
end iterator
Definition densematrix.hh:246
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::real_type real_type
Definition densematrix.hh:30
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::field_type field_type
Definition densematrix.hh:29
A dense n x m matrix.
Definition fmatrix.hh:67
Base::size_type size_type
Definition fmatrix.hh:80
vector space out of a tensor product of fields.
Definition fvector.hh:93
constexpr FieldVector()
Constructor making default-initialized vector.
Definition fvector.hh:113
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition densematrix.hh:73
Error thrown if operations of a FieldMatrix fail.
Definition densematrix.hh:121
Interface for a class of dense vectors over a given field.
Definition densevector.hh:235
size_type size() const
size method
Definition densevector.hh:297
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition densevector.hh:678
Generic iterator class for dense vector and matrix implementations.
Definition densevector.hh:128
Default exception class for mathematical errors.
Definition exceptions.hh:239
Definition ftraits.hh:24
T field_type
export the type representing the field
Definition ftraits.hh:26
T real_type
export the type representing the real type of the field
Definition ftraits.hh:28
Definition matvectraits.hh:29
static ctype pivoting_limit()
return threshold to do pivoting
Definition precision.hh:26
static ctype absolute_limit()
return threshold to declare matrix singular
Definition precision.hh:50
static ctype singular_limit()
return threshold to declare matrix singular
Definition precision.hh:38