# Dune Core Modules (2.6.0)

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 <utility>
11#include <vector>
12
17#include <dune/common/math.hh>
19#include <dune/common/simd.hh>
21#include <dune/common/unused.hh>
22
23namespace Dune
24{
25
26 template<typename M> class DenseMatrix;
27
28 template<typename M>
29 struct FieldTraits< DenseMatrix<M> >
30 {
31 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::field_type field_type;
32 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::real_type real_type;
33 };
34
35 /*
36 work around a problem of FieldMatrix/FieldVector,
37 there is no unique way to obtain the size of a class
38 */
39 template<class K, int N, int M> class FieldMatrix;
40 template<class K, int N> class FieldVector;
41 namespace {
42 template<class V>
43 struct VectorSize
44 {
45 static typename V::size_type size(const V & v) { return v.size(); }
46 };
47
48 template<class K, int N>
49 struct VectorSize< const FieldVector<K,N> >
50 {
51 typedef FieldVector<K,N> V;
52 static typename V::size_type size(const V & v)
53 {
55 return N;
56 }
57 };
58 }
59
78 template< class DenseMatrix, class RHS >
80
81#ifndef DOXYGEN
82 namespace Impl
83 {
84
85 template< class DenseMatrix, class RHS, class = void >
87 {};
88
89 template< class DenseMatrix, class RHS >
90 class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< Dune::IsNumber< RHS >::value > >
91 {
92 public:
93 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
94 {
95 typedef typename DenseMatrix::field_type field_type;
96 std::fill( denseMatrix.begin(), denseMatrix.end(), static_cast< field_type >( rhs ) );
97 }
98 };
99
100 template< class DenseMatrix, class RHS >
101 class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< !std::is_same< typename RHS::const_iterator, void >::value
102 && std::is_convertible< typename RHS::const_iterator::value_type, typename DenseMatrix::iterator::value_type >::value > >
103 {
104 public:
105 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
106 {
107 DUNE_ASSERT_BOUNDS(rhs.N() == denseMatrix.N());
108 DUNE_ASSERT_BOUNDS(rhs.M() == denseMatrix.M());
109 typename DenseMatrix::iterator tIt = std::begin(denseMatrix);
110 typename RHS::const_iterator sIt = std::begin(rhs);
111 for(; sIt != std::end(rhs); ++tIt, ++sIt)
112 std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt));
113 }
114 };
115
116 } // namespace Impl
117
118
119
120 template< class DenseMatrix, class RHS >
121 struct DenseMatrixAssigner
122 : public Impl::DenseMatrixAssigner< DenseMatrix, RHS >
123 {};
124
125
126 namespace Impl
127 {
128
129 template< class DenseMatrix, class RHS >
130 std::true_type hasDenseMatrixAssigner ( DenseMatrix &, const RHS &, decltype( Dune::DenseMatrixAssigner< DenseMatrix, RHS >::apply( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) ) * = nullptr );
131
132 std::false_type hasDenseMatrixAssigner ( ... );
133
134 } // namespace Impl
135
136 template< class DenseMatrix, class RHS >
137 struct HasDenseMatrixAssigner
138 : public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) )
139 {};
140
141#endif // #ifndef DOXYGEN
142
143
144
146 class FMatrixError : public MathError {};
147
158 template<typename MAT>
160 {
161 typedef DenseMatVecTraits<MAT> Traits;
162
163 // Curiously recurring template pattern
164 MAT & asImp() { return static_cast<MAT&>(*this); }
165 const MAT & asImp() const { return static_cast<const MAT&>(*this); }
166
167 public:
168 //===== type definitions and constants
169
171 typedef typename Traits::derived_type derived_type;
172
174 typedef typename Traits::value_type value_type;
175
177 typedef typename Traits::value_type field_type;
178
180 typedef typename Traits::value_type block_type;
181
183 typedef typename Traits::size_type size_type;
184
186 typedef typename Traits::row_type row_type;
187
189 typedef typename Traits::row_reference row_reference;
190
192 typedef typename Traits::const_row_reference const_row_reference;
193
195 enum {
197 blocklevel = 1
198 };
199
200 private:
203
206 using simd_index_type = SimdIndex<value_type>;
207
208 public:
210
213 {
214 return asImp().mat_access(i);
215 }
216
218 {
219 return asImp().mat_access(i);
220 }
221
224 {
225 return rows();
226 }
227
228 //===== iterator interface to rows of the matrix
236 typedef typename std::remove_reference<row_reference>::type::Iterator ColIterator;
237
240 {
241 return Iterator(*this,0);
242 }
243
246 {
247 return Iterator(*this,rows());
248 }
249
253 {
254 return Iterator(*this,rows()-1);
255 }
256
260 {
261 return Iterator(*this,-1);
262 }
263
271 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator ConstColIterator;
272
275 {
276 return ConstIterator(*this,0);
277 }
278
281 {
282 return ConstIterator(*this,rows());
283 }
284
288 {
289 return ConstIterator(*this,rows()-1);
290 }
291
295 {
296 return ConstIterator(*this,-1);
297 }
298
299 //===== assignment
300
301 template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > >
302 derived_type &operator= ( const RHS &rhs )
303 {
305 return asImp();
306 }
307
308 //===== vector space arithmetic
309
311 template <class Other>
313 {
314 DUNE_ASSERT_BOUNDS(rows() == y.rows());
315 for (size_type i=0; i<rows(); i++)
316 (*this)[i] += y[i];
317 return asImp();
318 }
319
321 template <class Other>
323 {
324 DUNE_ASSERT_BOUNDS(rows() == y.rows());
325 for (size_type i=0; i<rows(); i++)
326 (*this)[i] -= y[i];
327 return asImp();
328 }
329
332 {
333 for (size_type i=0; i<rows(); i++)
334 (*this)[i] *= k;
335 return asImp();
336 }
337
340 {
341 for (size_type i=0; i<rows(); i++)
342 (*this)[i] /= k;
343 return asImp();
344 }
345
347 template <class Other>
349 {
350 DUNE_ASSERT_BOUNDS(rows() == y.rows());
351 for( size_type i = 0; i < rows(); ++i )
352 (*this)[ i ].axpy( k, y[ i ] );
353 return asImp();
354 }
355
357 template <class Other>
358 bool operator== (const DenseMatrix<Other>& y) const
359 {
360 DUNE_ASSERT_BOUNDS(rows() == y.rows());
361 for (size_type i=0; i<rows(); i++)
362 if ((*this)[i]!=y[i])
363 return false;
364 return true;
365 }
367 template <class Other>
368 bool operator!= (const DenseMatrix<Other>& y) const
369 {
370 return !operator==(y);
371 }
372
373
374 //===== linear maps
375
377 template<class X, class Y>
378 void mv (const X& x, Y& y) const
379 {
380 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
381 DUNE_ASSERT_BOUNDS(x.N() == M());
382 DUNE_ASSERT_BOUNDS(y.N() == N());
383
384 using field_type = typename FieldTraits<Y>::field_type;
385 for (size_type i=0; i<rows(); ++i)
386 {
387 y[i] = field_type(0);
388 for (size_type j=0; j<cols(); j++)
389 y[i] += (*this)[i][j] * x[j];
390 }
391 }
392
394 template< class X, class Y >
395 void mtv ( const X &x, Y &y ) const
396 {
397 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
398 DUNE_ASSERT_BOUNDS(x.N() == N());
399 DUNE_ASSERT_BOUNDS(y.N() == M());
400
401 using field_type = typename FieldTraits<Y>::field_type;
402 for(size_type i = 0; i < cols(); ++i)
403 {
404 y[i] = field_type(0);
405 for(size_type j = 0; j < rows(); ++j)
406 y[i] += (*this)[j][i] * x[j];
407 }
408 }
409
411 template<class X, class Y>
412 void umv (const X& x, Y& y) const
413 {
414 DUNE_ASSERT_BOUNDS(x.N() == M());
415 DUNE_ASSERT_BOUNDS(y.N() == N());
416 for (size_type i=0; i<rows(); i++)
417 for (size_type j=0; j<cols(); j++)
418 y[i] += (*this)[i][j] * x[j];
419 }
420
422 template<class X, class Y>
423 void umtv (const X& x, Y& y) const
424 {
425 DUNE_ASSERT_BOUNDS(x.N() == N());
426 DUNE_ASSERT_BOUNDS(y.N() == M());
427 for (size_type i=0; i<rows(); i++)
428 for (size_type j=0; j<cols(); j++)
429 y[j] += (*this)[i][j]*x[i];
430 }
431
433 template<class X, class Y>
434 void umhv (const X& x, Y& y) const
435 {
436 DUNE_ASSERT_BOUNDS(x.N() == N());
437 DUNE_ASSERT_BOUNDS(y.N() == M());
438 for (size_type i=0; i<rows(); i++)
439 for (size_type j=0; j<cols(); j++)
440 y[j] += conjugateComplex((*this)[i][j])*x[i];
441 }
442
444 template<class X, class Y>
445 void mmv (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] -= (*this)[i][j] * x[j];
452 }
453
455 template<class X, class Y>
456 void mmtv (const X& x, Y& y) const
457 {
458 DUNE_ASSERT_BOUNDS(x.N() == N());
459 DUNE_ASSERT_BOUNDS(y.N() == M());
460 for (size_type i=0; i<rows(); i++)
461 for (size_type j=0; j<cols(); j++)
462 y[j] -= (*this)[i][j]*x[i];
463 }
464
466 template<class X, class Y>
467 void mmhv (const X& x, Y& y) const
468 {
469 DUNE_ASSERT_BOUNDS(x.N() == N());
470 DUNE_ASSERT_BOUNDS(y.N() == M());
471 for (size_type i=0; i<rows(); i++)
472 for (size_type j=0; j<cols(); j++)
473 y[j] -= conjugateComplex((*this)[i][j])*x[i];
474 }
475
477 template<class X, class Y>
478 void usmv (const typename FieldTraits<Y>::field_type & alpha,
479 const X& x, Y& y) const
480 {
481 DUNE_ASSERT_BOUNDS(x.N() == M());
482 DUNE_ASSERT_BOUNDS(y.N() == N());
483 for (size_type i=0; i<rows(); i++)
484 for (size_type j=0; j<cols(); j++)
485 y[i] += alpha * (*this)[i][j] * x[j];
486 }
487
489 template<class X, class Y>
490 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
491 const X& x, Y& y) const
492 {
493 DUNE_ASSERT_BOUNDS(x.N() == N());
494 DUNE_ASSERT_BOUNDS(y.N() == M());
495 for (size_type i=0; i<rows(); i++)
496 for (size_type j=0; j<cols(); j++)
497 y[j] += alpha*(*this)[i][j]*x[i];
498 }
499
501 template<class X, class Y>
502 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
503 const X& x, Y& y) const
504 {
505 DUNE_ASSERT_BOUNDS(x.N() == N());
506 DUNE_ASSERT_BOUNDS(y.N() == M());
507 for (size_type i=0; i<rows(); i++)
508 for (size_type j=0; j<cols(); j++)
509 y[j] += alpha*conjugateComplex((*this)[i][j])*x[i];
510 }
511
512 //===== norms
513
515 typename FieldTraits<value_type>::real_type frobenius_norm () const
516 {
517 typename FieldTraits<value_type>::real_type sum=(0.0);
518 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
519 return fvmeta::sqrt(sum);
520 }
521
523 typename FieldTraits<value_type>::real_type frobenius_norm2 () const
524 {
525 typename FieldTraits<value_type>::real_type sum=(0.0);
526 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
527 return sum;
528 }
529
531 template <typename vt = value_type,
532 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
533 typename FieldTraits<vt>::real_type infinity_norm() const {
534 using real_type = typename FieldTraits<vt>::real_type;
535 using std::max;
536
537 real_type norm = 0;
538 for (auto const &x : *this) {
539 real_type const a = x.one_norm();
540 norm = max(a, norm);
541 }
542 return norm;
543 }
544
546 template <typename vt = value_type,
547 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
548 typename FieldTraits<vt>::real_type infinity_norm_real() const {
549 using real_type = typename FieldTraits<vt>::real_type;
550 using std::max;
551
552 real_type norm = 0;
553 for (auto const &x : *this) {
554 real_type const a = x.one_norm_real();
555 norm = max(a, norm);
556 }
557 return norm;
558 }
559
561 template <typename vt = value_type,
562 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
563 typename FieldTraits<vt>::real_type infinity_norm() const {
564 using real_type = typename FieldTraits<vt>::real_type;
565 using std::max;
566
567 real_type norm = 0;
568 real_type isNaN = 1;
569 for (auto const &x : *this) {
570 real_type const a = x.one_norm();
571 norm = max(a, norm);
572 isNaN += a;
573 }
574 isNaN /= isNaN;
575 return norm * isNaN;
576 }
577
579 template <typename vt = value_type,
580 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
581 typename FieldTraits<vt>::real_type infinity_norm_real() const {
582 using real_type = typename FieldTraits<vt>::real_type;
583 using std::max;
584
585 real_type norm = 0;
586 real_type isNaN = 1;
587 for (auto const &x : *this) {
588 real_type const a = x.one_norm_real();
589 norm = max(a, norm);
590 isNaN += a;
591 }
592 isNaN /= isNaN;
593 return norm * isNaN;
594 }
595
596 //===== solve
597
602 template <class V>
603 void solve (V& x, const V& b) const;
604
609 void invert();
610
613
615 template<typename M2>
617 {
618 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
619 DUNE_ASSERT_BOUNDS(M.rows() == rows());
620 MAT C(asImp());
621
622 for (size_type i=0; i<rows(); i++)
623 for (size_type j=0; j<cols(); j++) {
624 (*this)[i][j] = 0;
625 for (size_type k=0; k<rows(); k++)
626 (*this)[i][j] += M[i][k]*C[k][j];
627 }
628
629 return asImp();
630 }
631
633 template<typename M2>
635 {
636 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
637 DUNE_ASSERT_BOUNDS(M.cols() == cols());
638 MAT C(asImp());
639
640 for (size_type i=0; i<rows(); i++)
641 for (size_type j=0; j<cols(); j++) {
642 (*this)[i][j] = 0;
643 for (size_type k=0; k<cols(); k++)
644 (*this)[i][j] += C[i][k]*M[k][j];
645 }
646 return asImp();
647 }
648
649#if 0
651 template<int l>
652 DenseMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
653 {
655
656 for (size_type i=0; i<l; i++) {
657 for (size_type j=0; j<cols(); j++) {
658 C[i][j] = 0;
659 for (size_type k=0; k<rows(); k++)
660 C[i][j] += M[i][k]*(*this)[k][j];
661 }
662 }
663 return C;
664 }
665
667 template<int l>
668 FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
669 {
670 FieldMatrix<K,rows,l> C;
671
672 for (size_type i=0; i<rows(); i++) {
673 for (size_type j=0; j<l; j++) {
674 C[i][j] = 0;
675 for (size_type k=0; k<cols(); k++)
676 C[i][j] += (*this)[i][k]*M[k][j];
677 }
678 }
679 return C;
680 }
681#endif
682
683 //===== sizes
684
686 size_type N () const
687 {
688 return rows();
689 }
690
692 size_type M () const
693 {
694 return cols();
695 }
696
699 {
700 return asImp().mat_rows();
701 }
702
705 {
706 return asImp().mat_cols();
707 }
708
709 //===== query
710
712 bool exists (size_type i, size_type j) const
713 {
716 DUNE_ASSERT_BOUNDS(i >= 0 && i < rows());
717 DUNE_ASSERT_BOUNDS(j >= 0 && j < cols());
718 return true;
719 }
720
721 private:
722
723#ifndef DOXYGEN
724 struct ElimPivot
725 {
726 ElimPivot(std::vector<simd_index_type> & pivot);
727
728 void swap(std::size_t i, simd_index_type j);
729
730 template<typename T>
731 void operator()(const T&, int, int)
732 {}
733
734 std::vector<simd_index_type> & pivot_;
735 };
736
737 template<typename V>
738 struct Elim
739 {
740 Elim(V& rhs);
741
742 void swap(std::size_t i, simd_index_type 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(std::size_t i, simd_index_type j)
755 {
756 sign_ *= cond(simd_index_type(i) == j, field_type(1), field_type(-1));
757 }
758
759 void operator()(const field_type&, int, int)
760 {}
761
762 field_type& sign_;
763 };
764#endif // DOXYGEN
765
767
805 void luDecomposition(DenseMatrix<MAT>& A, Func func,
806 Mask &nonsingularLanes, bool throwEarly) const;
807 };
808
809#ifndef DOXYGEN
810 template<typename MAT>
811 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<simd_index_type> & pivot)
812 : pivot_(pivot)
813 {
814 typedef typename std::vector<size_type>::size_type size_type;
815 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
816 }
817
818 template<typename MAT>
819 void DenseMatrix<MAT>::ElimPivot::swap(std::size_t i, simd_index_type j)
820 {
821 assign(pivot_[i], j, SimdScalar<simd_index_type>(i) != j);
822 }
823
824 template<typename MAT>
825 template<typename V>
826 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
827 : rhs_(&rhs)
828 {}
829
830 template<typename MAT>
831 template<typename V>
832 void DenseMatrix<MAT>::Elim<V>::swap(std::size_t i, simd_index_type j)
833 {
834 using std::swap;
835
836 // see the comment in luDecomposition()
837 for(std::size_t l = 0; l < lanes(j); ++l)
838 swap(lane(l, (*rhs_)[ i ]),
839 lane(l, (*rhs_)[lane(l, j)]));
840 }
841
842 template<typename MAT>
843 template<typename V>
844 void DenseMatrix<MAT>::
845 Elim<V>::operator()(const typename V::field_type& factor, int k, int i)
846 {
847 (*rhs_)[k] -= factor*(*rhs_)[i];
848 }
849
850 template<typename MAT>
852 inline void DenseMatrix<MAT>::
853 luDecomposition(DenseMatrix<MAT>& A, Func func, Mask &nonsingularLanes,
854 bool throwEarly) const
855 {
856 using std::swap;
857
858 typedef typename FieldTraits<value_type>::real_type real_type;
859
860 real_type norm = A.infinity_norm_real(); // for relative thresholds
861 real_type pivthres =
864 real_type singthres =
867
868 // LU decomposition of A in A
869 for (size_type i=0; i<rows(); i++) // loop over all rows
870 {
871 real_type pivmax = fvmeta::absreal(A[i][i]);
872 auto do_pivot = pivmax<pivthres;
873
874 // pivoting ?
875 if (any_true(do_pivot && nonsingularLanes))
876 {
877 // compute maximum of column
878 simd_index_type imax=i;
879 for (size_type k=i+1; k<rows(); k++)
880 {
881 auto abs = fvmeta::absreal(A[k][i]);
882 auto mask = abs > pivmax && do_pivot;
883 pivmax = cond(mask, abs, pivmax);
884 imax = cond(mask, simd_index_type(k), imax);
885 }
886 // swap rows
887 if (any_true(imax != i && nonsingularLanes)) {
888 for (size_type j=0; j<rows(); j++)
889 {
890 // This is a swap operation where the second operand is scattered,
891 // and on top of that is also extracted from deep within a
892 // moderately complicated data structure (a DenseMatrix), where we
893 // can't assume much on the memory layout. On intel processors,
894 // the only instruction that might help us here is vgather, but it
895 // is unclear whether that is even faster than a software
896 // implementation, and we would also need vscatter which does not
897 // exist. So break vectorization here and do it manually.
898 for(std::size_t l = 0; l < lanes(A[i][j]); ++l)
899 swap(lane(l, A[i][j]), lane(l, A[lane(l, imax)][j]));
900 }
901 func.swap(i, imax); // swap the pivot or rhs
902 }
903 }
904
905 // singular ?
906 nonsingularLanes = nonsingularLanes && !(pivmax<singthres);
907 if (throwEarly) {
908 if(!all_true(nonsingularLanes))
909 DUNE_THROW(FMatrixError, "matrix is singular");
910 }
911 else { // !throwEarly
912 if(!any_true(nonsingularLanes))
913 return;
914 }
915
916 // eliminate
917 for (size_type k=i+1; k<rows(); k++)
918 {
919 // in the simd case, A[i][i] may be close to zero in some lanes. Pray
920 // that the result is no worse than a quiet NaN.
921 field_type factor = A[k][i]/A[i][i];
922 A[k][i] = factor;
923 for (size_type j=i+1; j<rows(); j++)
924 A[k][j] -= factor*A[i][j];
925 func(factor, k, i);
926 }
927 }
928 }
929
930 template<typename MAT>
931 template <class V>
932 inline void DenseMatrix<MAT>::solve(V& x, const V& b) const
933 {
934 // never mind those ifs, because they get optimized away
935 if (rows()!=cols())
936 DUNE_THROW(FMatrixError, "Can't solve for a " << rows() << "x" << cols() << " matrix!");
937
938 if (rows()==1) {
939
940#ifdef DUNE_FMatrix_WITH_CHECKING
941 if (any_true(fvmeta::absreal((*this)[0][0])
943 DUNE_THROW(FMatrixError,"matrix is singular");
944#endif
945 x[0] = b[0]/(*this)[0][0];
946
947 }
948 else if (rows()==2) {
949
950 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
951#ifdef DUNE_FMatrix_WITH_CHECKING
952 if (any_true(fvmeta::absreal(detinv)
954 DUNE_THROW(FMatrixError,"matrix is singular");
955#endif
956 detinv = 1.0/detinv;
957
958 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
959 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
960
961 }
962 else if (rows()==3) {
963
964 field_type d = determinant();
965#ifdef DUNE_FMatrix_WITH_CHECKING
966 if (any_true(fvmeta::absreal(d) < FMatrixPrecision<>::absolute_limit()))
967 DUNE_THROW(FMatrixError,"matrix is singular");
968#endif
969
970 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
971 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
972 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
973
974 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
975 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
976 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
977
978 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
979 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
980 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
981
982 }
983 else {
984
985 V& rhs = x; // use x to store rhs
986 rhs = b; // copy data
987 Elim<V> elim(rhs);
988 MAT A(asImp());
990 nonsingularLanes(true);
991
992 luDecomposition(A, elim, nonsingularLanes, true);
993
994 // backsolve
995 for(int i=rows()-1; i>=0; i--) {
996 for (size_type j=i+1; j<rows(); j++)
997 rhs[i] -= A[i][j]*x[j];
998 x[i] = rhs[i]/A[i][i];
999 }
1000 }
1001 }
1002
1003 template<typename MAT>
1004 inline void DenseMatrix<MAT>::invert()
1005 {
1006 // never mind those ifs, because they get optimized away
1007 if (rows()!=cols())
1008 DUNE_THROW(FMatrixError, "Can't invert a " << rows() << "x" << cols() << " matrix!");
1009
1010 if (rows()==1) {
1011
1012#ifdef DUNE_FMatrix_WITH_CHECKING
1013 if (any_true(fvmeta::absreal((*this)[0][0])
1015 DUNE_THROW(FMatrixError,"matrix is singular");
1016#endif
1017 (*this)[0][0] = field_type( 1 ) / (*this)[0][0];
1018
1019 }
1020 else if (rows()==2) {
1021
1022 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
1023#ifdef DUNE_FMatrix_WITH_CHECKING
1024 if (any_true(fvmeta::absreal(detinv)
1026 DUNE_THROW(FMatrixError,"matrix is singular");
1027#endif
1028 detinv = field_type( 1 ) / detinv;
1029
1030 field_type temp=(*this)[0][0];
1031 (*this)[0][0] = (*this)[1][1]*detinv;
1032 (*this)[0][1] = -(*this)[0][1]*detinv;
1033 (*this)[1][0] = -(*this)[1][0]*detinv;
1034 (*this)[1][1] = temp*detinv;
1035
1036 }
1037 else if (rows()==3)
1038 {
1039 using K = field_type;
1040 // code generated by maple
1041 K t4 = (*this)[0][0] * (*this)[1][1];
1042 K t6 = (*this)[0][0] * (*this)[1][2];
1043 K t8 = (*this)[0][1] * (*this)[1][0];
1044 K t10 = (*this)[0][2] * (*this)[1][0];
1045 K t12 = (*this)[0][1] * (*this)[2][0];
1046 K t14 = (*this)[0][2] * (*this)[2][0];
1047
1048 K det = (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
1049 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
1050 K t17 = K(1.0)/det;
1051
1052 K matrix01 = (*this)[0][1];
1053 K matrix00 = (*this)[0][0];
1054 K matrix10 = (*this)[1][0];
1055 K matrix11 = (*this)[1][1];
1056
1057 (*this)[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1])*t17;
1058 (*this)[0][1] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1])*t17;
1059 (*this)[0][2] = (matrix01 * (*this)[1][2] - (*this)[0][2] * (*this)[1][1])*t17;
1060 (*this)[1][0] = -((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0])*t17;
1061 (*this)[1][1] = (matrix00 * (*this)[2][2] - t14) * t17;
1062 (*this)[1][2] = -(t6-t10) * t17;
1063 (*this)[2][0] = (matrix10 * (*this)[2][1] - matrix11 * (*this)[2][0]) * t17;
1064 (*this)[2][1] = -(matrix00 * (*this)[2][1] - t12) * t17;
1065 (*this)[2][2] = (t4-t8) * t17;
1066 }
1067 else {
1068
1069 MAT A(asImp());
1070 std::vector<simd_index_type> pivot(rows());
1072 nonsingularLanes(true);
1073 luDecomposition(A, ElimPivot(pivot), nonsingularLanes, true);
1074 DenseMatrix<MAT>& L=A;
1075 DenseMatrix<MAT>& U=A;
1076
1077 // initialize inverse
1078 *this=field_type();
1079
1080 for(size_type i=0; i<rows(); ++i)
1081 (*this)[i][i]=1;
1082
1083 // L Y = I; multiple right hand sides
1084 for (size_type i=0; i<rows(); i++)
1085 for (size_type j=0; j<i; j++)
1086 for (size_type k=0; k<rows(); k++)
1087 (*this)[i][k] -= L[i][j]*(*this)[j][k];
1088
1089 // U A^{-1} = Y
1090 for (size_type i=rows(); i>0;) {
1091 --i;
1092 for (size_type k=0; k<rows(); k++) {
1093 for (size_type j=i+1; j<rows(); j++)
1094 (*this)[i][k] -= U[i][j]*(*this)[j][k];
1095 (*this)[i][k] /= U[i][i];
1096 }
1097 }
1098
1099 for(size_type i=rows(); i>0; ) {
1100 --i;
1101 for(std::size_t l = 0; l < lanes((*this)[0][0]); ++l)
1102 {
1103 std::size_t pi = lane(l, pivot[i]);
1104 if(i!=pi)
1105 for(size_type j=0; j<rows(); ++j)
1106 std::swap(lane(l, (*this)[j][pi]),
1107 lane(l, (*this)[j][ i]));
1108 }
1109 }
1110 }
1111 }
1112
1113 // implementation of the determinant
1114 template<typename MAT>
1115 inline typename DenseMatrix<MAT>::field_type
1117 {
1118 // never mind those ifs, because they get optimized away
1119 if (rows()!=cols())
1120 DUNE_THROW(FMatrixError, "There is no determinant for a " << rows() << "x" << cols() << " matrix!");
1121
1122 if (rows()==1)
1123 return (*this)[0][0];
1124
1125 if (rows()==2)
1126 return (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
1127
1128 if (rows()==3) {
1129 // code generated by maple
1130 field_type t4 = (*this)[0][0] * (*this)[1][1];
1131 field_type t6 = (*this)[0][0] * (*this)[1][2];
1132 field_type t8 = (*this)[0][1] * (*this)[1][0];
1133 field_type t10 = (*this)[0][2] * (*this)[1][0];
1134 field_type t12 = (*this)[0][1] * (*this)[2][0];
1135 field_type t14 = (*this)[0][2] * (*this)[2][0];
1136
1137 return (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
1138 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
1139
1140 }
1141
1142 MAT A(asImp());
1143 field_type det;
1145 nonsingularLanes(true);
1146
1147 luDecomposition(A, ElimDet(det), nonsingularLanes, false);
1148 assign(det, field_type(0), !nonsingularLanes);
1149
1150 for (size_type i = 0; i < rows(); ++i)
1151 det *= A[i][i];
1152 return det;
1153 }
1154
1155#endif // DOXYGEN
1156
1157 namespace DenseMatrixHelp {
1158
1160 template <typename MAT, typename V1, typename V2>
1161 static inline void multAssign(const DenseMatrix<MAT> &matrix, const DenseVector<V1> & x, DenseVector<V2> & ret)
1162 {
1163 DUNE_ASSERT_BOUNDS(x.size() == matrix.cols());
1164 DUNE_ASSERT_BOUNDS(ret.size() == matrix.rows());
1165 typedef typename DenseMatrix<MAT>::size_type size_type;
1166
1167 for(size_type i=0; i<matrix.rows(); ++i)
1168 {
1169 ret[i] = 0.0;
1170 for(size_type j=0; j<matrix.cols(); ++j)
1171 {
1172 ret[i] += matrix[i][j]*x[j];
1173 }
1174 }
1175 }
1176
1177#if 0
1179 template <typename K, int rows, int cols>
1180 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
1181 {
1182 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
1183
1184 for(size_type i=0; i<cols(); ++i)
1185 {
1186 ret[i] = 0.0;
1187 for(size_type j=0; j<rows(); ++j)
1188 ret[i] += matrix[j][i]*x[j];
1189 }
1190 }
1191
1193 template <typename K, int rows, int cols>
1194 static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
1195 {
1196 FieldVector<K,rows> ret;
1197 multAssign(matrix,x,ret);
1198 return ret;
1199 }
1200
1202 template <typename K, int rows, int cols>
1203 static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x)
1204 {
1205 FieldVector<K,cols> ret;
1206 multAssignTransposed( matrix, x, ret );
1207 return ret;
1208 }
1209#endif
1210
1211 } // end namespace DenseMatrixHelp
1212
1214 template<typename MAT>
1215 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1216 {
1217 for (typename DenseMatrix<MAT>::size_type i=0; i<a.rows(); i++)
1218 s << a[i] << std::endl;
1219 return s;
1220 }
1221
1224} // end namespace Dune
1225
1226#endif
Macro for wrapping boundary checks.
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:128
A dense n x m matrix.
Definition: densematrix.hh:160
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:267
bool operator!=(const DenseMatrix< Other > &y) const
Binary matrix incomparison.
Definition: densematrix.hh:368
void mv(const X &x, Y &y) const
y = A x
Definition: densematrix.hh:378
Traits::value_type field_type
export the type representing the field
Definition: densematrix.hh:177
void solve(V &x, const V &b) const
Solve system A x = b.
size_type cols() const
number of columns
Definition: densematrix.hh:704
ConstIterator beforeEnd() const
Definition: densematrix.hh:287
bool operator==(const DenseMatrix< Other > &y) const
Binary matrix comparison.
Definition: densematrix.hh:358
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:456
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:548
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:271
ConstIterator beforeBegin() const
Definition: densematrix.hh:294
Traits::value_type block_type
export the type representing the components
Definition: densematrix.hh:180
derived_type & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition: densematrix.hh:348
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:395
size_type size() const
size method (number of rows)
Definition: densematrix.hh:223
size_type M() const
number of columns
Definition: densematrix.hh:692
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:634
Iterator end()
end iterator
Definition: densematrix.hh:245
Iterator beforeBegin()
Definition: densematrix.hh:259
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:339
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:331
Iterator beforeEnd()
Definition: densematrix.hh:252
Traits::value_type value_type
export the type representing the field
Definition: densematrix.hh:174
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: densematrix.hh:478
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: densematrix.hh:502
size_type rows() const
number of rows
Definition: densematrix.hh:698
void mmv(const X &x, Y &y) const
y -= A x
Definition: densematrix.hh:445
void invert()
Compute inverse.
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:616
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: densematrix.hh:490
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: densematrix.hh:467
Traits::derived_type derived_type
type of derived matrix class
Definition: densematrix.hh:171
row_reference operator[](size_type i)
random access
Definition: densematrix.hh:212
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:712
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:234
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:515
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:412
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition: densematrix.hh:265
derived_type & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition: densematrix.hh:322
size_type N() const
number of rows
Definition: densematrix.hh:686
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:533
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:186
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:183
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:192
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:523
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition: densematrix.hh:236
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:189
Iterator iterator
typedef for stl compliant access
Definition: densematrix.hh:232
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: densematrix.hh:269
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:230
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:423
ConstIterator begin() const
begin iterator
Definition: densematrix.hh:274
@ blocklevel
The number of block levels we contain. This is 1.
Definition: densematrix.hh:197
Iterator begin()
begin iterator
Definition: densematrix.hh:239
void umhv(const X &x, Y &y) const
y += A^H x
Definition: densematrix.hh:434
derived_type & operator+=(const DenseMatrix< Other > &y)
Definition: densematrix.hh:312
ConstIterator end() const
end iterator
Definition: densematrix.hh:280
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
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:146
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:68
vector space out of a tensor product of fields.
Definition: fvector.hh:93
constexpr FieldVector()
Constructor making default-initialized vector.
Definition: fvector.hh:113
Default exception class for mathematical errors.
Definition: exceptions.hh:239
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:1006
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:1161
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:28
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:58
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:10
typename SimdIndexTypeTraits< V >::type SimdIndex
An simd vector of indices corresponding to a simd vector V.
Definition: simd.hh:220
T lane(std::size_t l, const T &v)
access a lane of a simd vector (scalar version)
Definition: simd.hh:340
void assign(T &dst, const T &src, bool mask)
Definition: simd.hh:421
const T1 cond(bool b, const T1 &v1, const T2 &v2)
conditional evaluate
Definition: conditional.hh:26
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:121
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:105
std::size_t lanes(const T &)
get the number of lanes of a simd vector (scalar version)
Definition: simd.hh:336
STL namespace.
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:79
Traits for type conversions and type information.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
|  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Aug 11, 22:30, 2024)