Dune Core Modules (2.4.1)

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