dune-common 2.9.0
Loading...
Searching...
No Matches
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// SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_DENSEMATRIX_HH
6#define DUNE_DENSEMATRIX_HH
7
8#include <cmath>
9#include <cstddef>
10#include <iostream>
11#include <type_traits>
12#include <utility>
13#include <vector>
14
19#include <dune/common/math.hh>
24
25namespace Dune
26{
27
28 template<typename M> class DenseMatrix;
29
30 template<typename M>
36
37 template<class K, int N, int M> class FieldMatrix;
38 template<class K, int N> class FieldVector;
39
58 template< class DenseMatrix, class RHS >
60
61#ifndef DOXYGEN
62 namespace Impl
63 {
64
65 template< class DenseMatrix, class RHS, class = void >
67 {};
68
69 template< class DenseMatrix, class RHS >
70 class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< Dune::IsNumber< RHS >::value > >
71 {
72 public:
73 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
74 {
75 typedef typename DenseMatrix::field_type field_type;
76 std::fill( denseMatrix.begin(), denseMatrix.end(), static_cast< field_type >( rhs ) );
77 }
78 };
79
80 template< class DenseMatrix, class RHS >
81 class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< !std::is_same< typename RHS::const_iterator, void >::value
82 && std::is_convertible< typename RHS::const_iterator::value_type, typename DenseMatrix::iterator::value_type >::value > >
83 {
84 public:
85 static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
86 {
87 DUNE_ASSERT_BOUNDS(rhs.N() == denseMatrix.N());
88 DUNE_ASSERT_BOUNDS(rhs.M() == denseMatrix.M());
89 typename DenseMatrix::iterator tIt = std::begin(denseMatrix);
90 typename RHS::const_iterator sIt = std::begin(rhs);
91 for(; sIt != std::end(rhs); ++tIt, ++sIt)
92 std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt));
93 }
94 };
95
96 } // namespace Impl
97
98
99
100 template< class DenseMatrix, class RHS >
101 struct DenseMatrixAssigner
102 : public Impl::DenseMatrixAssigner< DenseMatrix, RHS >
103 {};
104
105
106 namespace Impl
107 {
108
109 template< class DenseMatrix, class RHS >
110 std::true_type hasDenseMatrixAssigner ( DenseMatrix &, const RHS &, decltype( Dune::DenseMatrixAssigner< DenseMatrix, RHS >::apply( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) ) * = nullptr );
111
112 std::false_type hasDenseMatrixAssigner ( ... );
113
114 } // namespace Impl
115
116 template< class DenseMatrix, class RHS >
117 struct HasDenseMatrixAssigner
118 : public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) )
119 {};
120
121#endif // #ifndef DOXYGEN
122
123
124
126 class FMatrixError : public MathError {};
127
138 template<typename MAT>
140 {
142
143 // Curiously recurring template pattern
144 constexpr MAT & asImp() { return static_cast<MAT&>(*this); }
145 constexpr const MAT & asImp() const { return static_cast<const MAT&>(*this); }
146
147 template <class>
148 friend class DenseMatrix;
149
150 public:
151 //===== type definitions and constants
152
154 typedef typename Traits::derived_type derived_type;
155
157 typedef typename Traits::value_type value_type;
158
160 typedef typename Traits::value_type field_type;
161
163 typedef typename Traits::value_type block_type;
164
166 typedef typename Traits::size_type size_type;
167
169 typedef typename Traits::row_type row_type;
170
172 typedef typename Traits::row_reference row_reference;
173
175 typedef typename Traits::const_row_reference const_row_reference;
176
178 constexpr static int blocklevel = 1;
179
180 private:
183 using simd_index_type = Simd::Rebind<std::size_t, value_type>;
184
185 public:
186 //===== access to components
187
190 {
191 return asImp().mat_access(i);
192 }
193
195 {
196 return asImp().mat_access(i);
197 }
198
201 {
202 return rows();
203 }
204
205 //===== iterator interface to rows of the matrix
213 typedef typename std::remove_reference<row_reference>::type::Iterator ColIterator;
214
217 {
218 return Iterator(*this,0);
219 }
220
223 {
224 return Iterator(*this,rows());
225 }
226
230 {
231 return Iterator(*this,rows()-1);
232 }
233
237 {
238 return Iterator(*this,-1);
239 }
240
248 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator ConstColIterator;
249
252 {
253 return ConstIterator(*this,0);
254 }
255
258 {
259 return ConstIterator(*this,rows());
260 }
261
265 {
266 return ConstIterator(*this,rows()-1);
267 }
268
272 {
273 return ConstIterator(*this,-1);
274 }
275
276 //===== assignment
277
278 template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > >
279 derived_type &operator= ( const RHS &rhs )
280 {
282 return asImp();
283 }
284
285 //===== vector space arithmetic
286
288 template <class Other>
290 {
291 DUNE_ASSERT_BOUNDS(rows() == x.rows());
292 for (size_type i=0; i<rows(); i++)
293 (*this)[i] += x[i];
294 return asImp();
295 }
296
299 {
300 MAT result;
301 using idx_type = typename decltype(result)::size_type;
302
303 for (idx_type i = 0; i < rows(); ++i)
304 for (idx_type j = 0; j < cols(); ++j)
305 result[i][j] = - asImp()[i][j];
306
307 return result;
308 }
309
311 template <class Other>
313 {
314 DUNE_ASSERT_BOUNDS(rows() == x.rows());
315 for (size_type i=0; i<rows(); i++)
316 (*this)[i] -= x[i];
317 return asImp();
318 }
319
322 {
323 for (size_type i=0; i<rows(); i++)
324 (*this)[i] *= k;
325 return asImp();
326 }
327
330 {
331 for (size_type i=0; i<rows(); i++)
332 (*this)[i] /= k;
333 return asImp();
334 }
335
337 template <class Other>
339 {
340 DUNE_ASSERT_BOUNDS(rows() == x.rows());
341 for( size_type i = 0; i < rows(); ++i )
342 (*this)[ i ].axpy( a, x[ i ] );
343 return asImp();
344 }
345
347 template <class Other>
348 bool operator== (const DenseMatrix<Other>& x) const
349 {
350 DUNE_ASSERT_BOUNDS(rows() == x.rows());
351 for (size_type i=0; i<rows(); i++)
352 if ((*this)[i]!=x[i])
353 return false;
354 return true;
355 }
357 template <class Other>
358 bool operator!= (const DenseMatrix<Other>& x) const
359 {
360 return !operator==(x);
361 }
362
363
364 //===== linear maps
365
367 template<class X, class Y>
368 void mv (const X& x, Y& y) const
369 {
370 auto&& xx = Impl::asVector(x);
371 auto&& yy = Impl::asVector(y);
372 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
373 DUNE_ASSERT_BOUNDS(xx.N() == M());
374 DUNE_ASSERT_BOUNDS(yy.N() == N());
375
376 using y_field_type = typename FieldTraits<Y>::field_type;
377 for (size_type i=0; i<rows(); ++i)
378 {
379 yy[i] = y_field_type(0);
380 for (size_type j=0; j<cols(); j++)
381 yy[i] += (*this)[i][j] * xx[j];
382 }
383 }
384
386 template< class X, class Y >
387 void mtv ( const X &x, Y &y ) const
388 {
389 auto&& xx = Impl::asVector(x);
390 auto&& yy = Impl::asVector(y);
391 DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
392 DUNE_ASSERT_BOUNDS(xx.N() == N());
393 DUNE_ASSERT_BOUNDS(yy.N() == M());
394
395 using y_field_type = typename FieldTraits<Y>::field_type;
396 for(size_type i = 0; i < cols(); ++i)
397 {
398 yy[i] = y_field_type(0);
399 for(size_type j = 0; j < rows(); ++j)
400 yy[i] += (*this)[j][i] * xx[j];
401 }
402 }
403
405 template<class X, class Y>
406 void umv (const X& x, Y& y) const
407 {
408 auto&& xx = Impl::asVector(x);
409 auto&& yy = Impl::asVector(y);
410 DUNE_ASSERT_BOUNDS(xx.N() == M());
411 DUNE_ASSERT_BOUNDS(yy.N() == N());
412 for (size_type i=0; i<rows(); ++i)
413 for (size_type j=0; j<cols(); j++)
414 yy[i] += (*this)[i][j] * xx[j];
415 }
416
418 template<class X, class Y>
419 void umtv (const X& x, Y& y) const
420 {
421 auto&& xx = Impl::asVector(x);
422 auto&& yy = Impl::asVector(y);
423 DUNE_ASSERT_BOUNDS(xx.N() == N());
424 DUNE_ASSERT_BOUNDS(yy.N() == M());
425 for(size_type i = 0; i<rows(); ++i)
426 for (size_type j=0; j<cols(); j++)
427 yy[j] += (*this)[i][j]*xx[i];
428 }
429
431 template<class X, class Y>
432 void umhv (const X& x, Y& y) const
433 {
434 auto&& xx = Impl::asVector(x);
435 auto&& yy = Impl::asVector(y);
436 DUNE_ASSERT_BOUNDS(xx.N() == N());
437 DUNE_ASSERT_BOUNDS(yy.N() == M());
438 for (size_type i=0; i<rows(); i++)
439 for (size_type j=0; j<cols(); j++)
440 yy[j] += conjugateComplex((*this)[i][j])*xx[i];
441 }
442
444 template<class X, class Y>
445 void mmv (const X& x, Y& y) const
446 {
447 auto&& xx = Impl::asVector(x);
448 auto&& yy = Impl::asVector(y);
449 DUNE_ASSERT_BOUNDS(xx.N() == M());
450 DUNE_ASSERT_BOUNDS(yy.N() == N());
451 for (size_type i=0; i<rows(); i++)
452 for (size_type j=0; j<cols(); j++)
453 yy[i] -= (*this)[i][j] * xx[j];
454 }
455
457 template<class X, class Y>
458 void mmtv (const X& x, Y& y) const
459 {
460 auto&& xx = Impl::asVector(x);
461 auto&& yy = Impl::asVector(y);
462 DUNE_ASSERT_BOUNDS(xx.N() == N());
463 DUNE_ASSERT_BOUNDS(yy.N() == M());
464 for (size_type i=0; i<rows(); i++)
465 for (size_type j=0; j<cols(); j++)
466 yy[j] -= (*this)[i][j]*xx[i];
467 }
468
470 template<class X, class Y>
471 void mmhv (const X& x, Y& y) const
472 {
473 auto&& xx = Impl::asVector(x);
474 auto&& yy = Impl::asVector(y);
475 DUNE_ASSERT_BOUNDS(xx.N() == N());
476 DUNE_ASSERT_BOUNDS(yy.N() == M());
477 for (size_type i=0; i<rows(); i++)
478 for (size_type j=0; j<cols(); j++)
479 yy[j] -= conjugateComplex((*this)[i][j])*xx[i];
480 }
481
483 template<class X, class Y>
484 void usmv (const typename FieldTraits<Y>::field_type & alpha,
485 const X& x, Y& y) const
486 {
487 auto&& xx = Impl::asVector(x);
488 auto&& yy = Impl::asVector(y);
489 DUNE_ASSERT_BOUNDS(xx.N() == M());
490 DUNE_ASSERT_BOUNDS(yy.N() == N());
491 for (size_type i=0; i<rows(); i++)
492 for (size_type j=0; j<cols(); j++)
493 yy[i] += alpha * (*this)[i][j] * xx[j];
494 }
495
497 template<class X, class Y>
498 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
499 const X& x, Y& y) const
500 {
501 auto&& xx = Impl::asVector(x);
502 auto&& yy = Impl::asVector(y);
503 DUNE_ASSERT_BOUNDS(xx.N() == N());
504 DUNE_ASSERT_BOUNDS(yy.N() == M());
505 for (size_type i=0; i<rows(); i++)
506 for (size_type j=0; j<cols(); j++)
507 yy[j] += alpha*(*this)[i][j]*xx[i];
508 }
509
511 template<class X, class Y>
512 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
513 const X& x, Y& y) const
514 {
515 auto&& xx = Impl::asVector(x);
516 auto&& yy = Impl::asVector(y);
517 DUNE_ASSERT_BOUNDS(xx.N() == N());
518 DUNE_ASSERT_BOUNDS(yy.N() == M());
519 for (size_type i=0; i<rows(); i++)
520 for (size_type j=0; j<cols(); j++)
521 yy[j] +=
522 alpha*conjugateComplex((*this)[i][j])*xx[i];
523 }
524
525 //===== norms
526
529 {
530 typename FieldTraits<value_type>::real_type sum=(0.0);
531 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
532 return fvmeta::sqrt(sum);
533 }
534
537 {
538 typename FieldTraits<value_type>::real_type sum=(0.0);
539 for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
540 return sum;
541 }
542
544 template <typename vt = value_type,
545 typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
547 using real_type = typename FieldTraits<vt>::real_type;
548 using std::max;
549
550 real_type norm = 0;
551 for (auto const &x : *this) {
552 real_type const a = x.one_norm();
553 norm = max(a, norm);
554 }
555 return norm;
556 }
557
559 template <typename vt = value_type,
560 typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
562 using real_type = typename FieldTraits<vt>::real_type;
563 using std::max;
564
565 real_type norm = 0;
566 for (auto const &x : *this) {
567 real_type const a = x.one_norm_real();
568 norm = max(a, norm);
569 }
570 return norm;
571 }
572
574 template <typename vt = value_type,
575 typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
577 using real_type = typename FieldTraits<vt>::real_type;
578 using std::max;
579
580 real_type norm = 0;
581 real_type isNaN = 1;
582 for (auto const &x : *this) {
583 real_type const a = x.one_norm();
584 norm = max(a, norm);
585 isNaN += a;
586 }
587 return norm * (isNaN / isNaN);
588 }
589
591 template <typename vt = value_type,
592 typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
594 using real_type = typename FieldTraits<vt>::real_type;
595 using std::max;
596
597 real_type norm = 0;
598 real_type isNaN = 1;
599 for (auto const &x : *this) {
600 real_type const a = x.one_norm_real();
601 norm = max(a, norm);
602 isNaN += a;
603 }
604 return norm * (isNaN / isNaN);
605 }
606
607 //===== solve
608
613 template <class V1, class V2>
614 void solve (V1& x, const V2& b, bool doPivoting = true) const;
615
620 void invert(bool doPivoting = true);
621
623 field_type determinant (bool doPivoting = true) const;
624
626 template<typename M2>
628 {
629 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
630 DUNE_ASSERT_BOUNDS(M.rows() == rows());
631 AutonomousValue<MAT> C(asImp());
632
633 for (size_type i=0; i<rows(); i++)
634 for (size_type j=0; j<cols(); j++) {
635 (*this)[i][j] = 0;
636 for (size_type k=0; k<rows(); k++)
637 (*this)[i][j] += M[i][k]*C[k][j];
638 }
639
640 return asImp();
641 }
642
644 template<typename M2>
646 {
647 DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
648 DUNE_ASSERT_BOUNDS(M.cols() == cols());
649 AutonomousValue<MAT> C(asImp());
650
651 for (size_type i=0; i<rows(); i++)
652 for (size_type j=0; j<cols(); j++) {
653 (*this)[i][j] = 0;
654 for (size_type k=0; k<cols(); k++)
655 (*this)[i][j] += C[i][k]*M[k][j];
656 }
657 return asImp();
658 }
659
660#if 0
662 template<int l>
663 DenseMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
664 {
666
667 for (size_type i=0; i<l; i++) {
668 for (size_type j=0; j<cols(); j++) {
669 C[i][j] = 0;
670 for (size_type k=0; k<rows(); k++)
671 C[i][j] += M[i][k]*(*this)[k][j];
672 }
673 }
674 return C;
675 }
676
678 template<int l>
679 FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
680 {
681 FieldMatrix<K,rows,l> C;
682
683 for (size_type i=0; i<rows(); i++) {
684 for (size_type j=0; j<l; j++) {
685 C[i][j] = 0;
686 for (size_type k=0; k<cols(); k++)
687 C[i][j] += (*this)[i][k]*M[k][j];
688 }
689 }
690 return C;
691 }
692#endif
693
694 //===== sizes
695
697 constexpr size_type N () const
698 {
699 return rows();
700 }
701
703 constexpr size_type M () const
704 {
705 return cols();
706 }
707
709 constexpr size_type rows() const
710 {
711 return asImp().mat_rows();
712 }
713
715 constexpr size_type cols() const
716 {
717 return asImp().mat_cols();
718 }
719
720 //===== query
721
723 bool exists ([[maybe_unused]] size_type i, [[maybe_unused]] size_type j) const
724 {
725 DUNE_ASSERT_BOUNDS(i >= 0 && i < rows());
726 DUNE_ASSERT_BOUNDS(j >= 0 && j < cols());
727 return true;
728 }
729
730 protected:
731
732#ifndef DOXYGEN
733 struct ElimPivot
734 {
735 ElimPivot(std::vector<simd_index_type> & pivot);
736
737 void swap(std::size_t i, simd_index_type j);
738
739 template<typename T>
740 void operator()(const T&, int, int)
741 {}
742
743 std::vector<simd_index_type> & pivot_;
744 };
745
746 template<typename V>
747 struct Elim
748 {
749 Elim(V& rhs);
750
751 void swap(std::size_t i, simd_index_type j);
752
753 void operator()(const typename V::field_type& factor, int k, int i);
754
755 V* rhs_;
756 };
757
758 struct ElimDet
759 {
760 ElimDet(field_type& sign) : sign_(sign)
761 { sign_ = 1; }
762
763 void swap(std::size_t i, simd_index_type j)
764 {
765 sign_ *=
766 Simd::cond(simd_index_type(i) == j, field_type(1), field_type(-1));
767 }
768
769 void operator()(const field_type&, int, int)
770 {}
771
772 field_type& sign_;
773 };
774#endif // DOXYGEN
775
777
815 template<class Func, class Mask>
816 static void luDecomposition(DenseMatrix<MAT>& A, Func func,
817 Mask &nonsingularLanes, bool throwEarly, bool doPivoting);
818 };
819
820#ifndef DOXYGEN
821 template<typename MAT>
822 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<simd_index_type> & pivot)
823 : pivot_(pivot)
824 {
825 typedef typename std::vector<size_type>::size_type size_type;
826 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
827 }
828
829 template<typename MAT>
830 void DenseMatrix<MAT>::ElimPivot::swap(std::size_t i, simd_index_type j)
831 {
832 pivot_[i] =
833 Simd::cond(Simd::Scalar<simd_index_type>(i) == j, pivot_[i], j);
834 }
835
836 template<typename MAT>
837 template<typename V>
838 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
839 : rhs_(&rhs)
840 {}
841
842 template<typename MAT>
843 template<typename V>
844 void DenseMatrix<MAT>::Elim<V>::swap(std::size_t i, simd_index_type j)
845 {
846 using std::swap;
847
848 // see the comment in luDecomposition()
849 for(std::size_t l = 0; l < Simd::lanes(j); ++l)
850 swap(Simd::lane(l, (*rhs_)[ i ]),
851 Simd::lane(l, (*rhs_)[Simd::lane(l, j)]));
852 }
853
854 template<typename MAT>
855 template<typename V>
856 void DenseMatrix<MAT>::
857 Elim<V>::operator()(const typename V::field_type& factor, int k, int i)
858 {
859 (*rhs_)[k] -= factor*(*rhs_)[i];
860 }
861
862 template<typename MAT>
863 template<typename Func, class Mask>
864 inline void DenseMatrix<MAT>::
865 luDecomposition(DenseMatrix<MAT>& A, Func func, Mask &nonsingularLanes,
866 bool throwEarly, bool doPivoting)
867 {
868 using std::max;
869 using std::swap;
870
871 typedef typename FieldTraits<value_type>::real_type real_type;
872
873 // LU decomposition of A in A
874 for (size_type i=0; i<A.rows(); i++) // loop over all rows
875 {
876 real_type pivmax = fvmeta::absreal(A[i][i]);
877
878 if (doPivoting)
879 {
880 // compute maximum of column
881 simd_index_type imax=i;
882 for (size_type k=i+1; k<A.rows(); k++)
883 {
884 auto abs = fvmeta::absreal(A[k][i]);
885 auto mask = abs > pivmax;
886 pivmax = Simd::cond(mask, abs, pivmax);
887 imax = Simd::cond(mask, simd_index_type(k), imax);
888 }
889 // swap rows
890 for (size_type j=0; j<A.rows(); j++)
891 {
892 // This is a swap operation where the second operand is scattered,
893 // and on top of that is also extracted from deep within a
894 // moderately complicated data structure (a DenseMatrix), where we
895 // can't assume much on the memory layout. On intel processors,
896 // the only instruction that might help us here is vgather, but it
897 // is unclear whether that is even faster than a software
898 // implementation, and we would also need vscatter which does not
899 // exist. So break vectorization here and do it manually.
900 for(std::size_t l = 0; l < Simd::lanes(A[i][j]); ++l)
901 swap(Simd::lane(l, A[i][j]),
902 Simd::lane(l, A[Simd::lane(l, imax)][j]));
903 }
904 func.swap(i, imax); // swap the pivot or rhs
905 }
906
907 // singular ?
908 nonsingularLanes = nonsingularLanes && (pivmax != real_type(0));
909 if (throwEarly) {
910 if(!Simd::allTrue(nonsingularLanes))
911 DUNE_THROW(FMatrixError, "matrix is singular");
912 }
913 else { // !throwEarly
914 if(!Simd::anyTrue(nonsingularLanes))
915 return;
916 }
917
918 // eliminate
919 for (size_type k=i+1; k<A.rows(); k++)
920 {
921 // in the simd case, A[i][i] may be close to zero in some lanes. Pray
922 // that the result is no worse than a quiet NaN.
923 field_type factor = A[k][i]/A[i][i];
924 A[k][i] = factor;
925 for (size_type j=i+1; j<A.rows(); j++)
926 A[k][j] -= factor*A[i][j];
927 func(factor, k, i);
928 }
929 }
930 }
931
932 template<typename MAT>
933 template <class V1, class V2>
934 inline void DenseMatrix<MAT>::solve(V1& x, const V2& b, bool doPivoting) const
935 {
936 using real_type = typename FieldTraits<value_type>::real_type;
937 // never mind those ifs, because they get optimized away
938 if (rows()!=cols())
939 DUNE_THROW(FMatrixError, "Can't solve for a " << rows() << "x" << cols() << " matrix!");
940
941 if (rows()==1) {
942
943#ifdef DUNE_FMatrix_WITH_CHECKING
944 if (Simd::anyTrue(fvmeta::absreal((*this)[0][0])
946 DUNE_THROW(FMatrixError,"matrix is singular");
947#endif
948 x[0] = b[0]/(*this)[0][0];
949
950 }
951 else if (rows()==2) {
952
953 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
954#ifdef DUNE_FMatrix_WITH_CHECKING
955 if (Simd::anyTrue(fvmeta::absreal(detinv)
957 DUNE_THROW(FMatrixError,"matrix is singular");
958#endif
959 detinv = real_type(1.0)/detinv;
960
961 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
962 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
963
964 }
965 else if (rows()==3) {
966
967 field_type d = determinant(doPivoting);
968#ifdef DUNE_FMatrix_WITH_CHECKING
969 if (Simd::anyTrue(fvmeta::absreal(d)
971 DUNE_THROW(FMatrixError,"matrix is singular");
972#endif
973
974 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
975 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
976 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
977
978 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
979 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
980 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
981
982 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
983 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
984 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
985
986 }
987 else {
988
989 V1& rhs = x; // use x to store rhs
990 rhs = b; // copy data
991 Elim<V1> elim(rhs);
992 AutonomousValue<MAT> A(asImp());
993 Simd::Mask<typename FieldTraits<value_type>::real_type>
994 nonsingularLanes(true);
995
996 AutonomousValue<MAT>::luDecomposition(A, elim, nonsingularLanes, true, doPivoting);
997
998 // backsolve
999 for(int i=rows()-1; i>=0; i--) {
1000 for (size_type j=i+1; j<rows(); j++)
1001 rhs[i] -= A[i][j]*x[j];
1002 x[i] = rhs[i]/A[i][i];
1003 }
1004 }
1005 }
1006
1007 template<typename MAT>
1008 inline void DenseMatrix<MAT>::invert(bool doPivoting)
1009 {
1010 using real_type = typename FieldTraits<MAT>::real_type;
1011 using std::swap;
1012
1013 // never mind those ifs, because they get optimized away
1014 if (rows()!=cols())
1015 DUNE_THROW(FMatrixError, "Can't invert a " << rows() << "x" << cols() << " matrix!");
1016
1017 if (rows()==1) {
1018
1019#ifdef DUNE_FMatrix_WITH_CHECKING
1020 if (Simd::anyTrue(fvmeta::absreal((*this)[0][0])
1022 DUNE_THROW(FMatrixError,"matrix is singular");
1023#endif
1024 (*this)[0][0] = real_type( 1 ) / (*this)[0][0];
1025
1026 }
1027 else if (rows()==2) {
1028
1029 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
1030#ifdef DUNE_FMatrix_WITH_CHECKING
1031 if (Simd::anyTrue(fvmeta::absreal(detinv)
1033 DUNE_THROW(FMatrixError,"matrix is singular");
1034#endif
1035 detinv = real_type( 1 ) / detinv;
1036
1037 field_type temp=(*this)[0][0];
1038 (*this)[0][0] = (*this)[1][1]*detinv;
1039 (*this)[0][1] = -(*this)[0][1]*detinv;
1040 (*this)[1][0] = -(*this)[1][0]*detinv;
1041 (*this)[1][1] = temp*detinv;
1042
1043 }
1044 else if (rows()==3)
1045 {
1046 using K = field_type;
1047 // code generated by maple
1048 K t4 = (*this)[0][0] * (*this)[1][1];
1049 K t6 = (*this)[0][0] * (*this)[1][2];
1050 K t8 = (*this)[0][1] * (*this)[1][0];
1051 K t10 = (*this)[0][2] * (*this)[1][0];
1052 K t12 = (*this)[0][1] * (*this)[2][0];
1053 K t14 = (*this)[0][2] * (*this)[2][0];
1054
1055 K det = (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
1056 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
1057 K t17 = K(1.0)/det;
1058
1059 K matrix01 = (*this)[0][1];
1060 K matrix00 = (*this)[0][0];
1061 K matrix10 = (*this)[1][0];
1062 K matrix11 = (*this)[1][1];
1063
1064 (*this)[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1])*t17;
1065 (*this)[0][1] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1])*t17;
1066 (*this)[0][2] = (matrix01 * (*this)[1][2] - (*this)[0][2] * (*this)[1][1])*t17;
1067 (*this)[1][0] = -((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0])*t17;
1068 (*this)[1][1] = (matrix00 * (*this)[2][2] - t14) * t17;
1069 (*this)[1][2] = -(t6-t10) * t17;
1070 (*this)[2][0] = (matrix10 * (*this)[2][1] - matrix11 * (*this)[2][0]) * t17;
1071 (*this)[2][1] = -(matrix00 * (*this)[2][1] - t12) * t17;
1072 (*this)[2][2] = (t4-t8) * t17;
1073 }
1074 else {
1075 using std::swap;
1076
1077 AutonomousValue<MAT> A(asImp());
1078 std::vector<simd_index_type> pivot(rows());
1079 Simd::Mask<typename FieldTraits<value_type>::real_type>
1080 nonsingularLanes(true);
1081 AutonomousValue<MAT>::luDecomposition(A, ElimPivot(pivot), nonsingularLanes, true, doPivoting);
1082 auto& L=A;
1083 auto& U=A;
1084
1085 // initialize inverse
1086 *this=field_type();
1087
1088 for(size_type i=0; i<rows(); ++i)
1089 (*this)[i][i]=1;
1090
1091 // L Y = I; multiple right hand sides
1092 for (size_type i=0; i<rows(); i++)
1093 for (size_type j=0; j<i; j++)
1094 for (size_type k=0; k<rows(); k++)
1095 (*this)[i][k] -= L[i][j]*(*this)[j][k];
1096
1097 // U A^{-1} = Y
1098 for (size_type i=rows(); i>0;) {
1099 --i;
1100 for (size_type k=0; k<rows(); k++) {
1101 for (size_type j=i+1; j<rows(); j++)
1102 (*this)[i][k] -= U[i][j]*(*this)[j][k];
1103 (*this)[i][k] /= U[i][i];
1104 }
1105 }
1106
1107 for(size_type i=rows(); i>0; ) {
1108 --i;
1109 for(std::size_t l = 0; l < Simd::lanes((*this)[0][0]); ++l)
1110 {
1111 std::size_t pi = Simd::lane(l, pivot[i]);
1112 if(i!=pi)
1113 for(size_type j=0; j<rows(); ++j)
1114 swap(Simd::lane(l, (*this)[j][pi]),
1115 Simd::lane(l, (*this)[j][ i]));
1116 }
1117 }
1118 }
1119 }
1120
1121 // implementation of the determinant
1122 template<typename MAT>
1123 inline typename DenseMatrix<MAT>::field_type
1124 DenseMatrix<MAT>::determinant(bool doPivoting) const
1125 {
1126 // never mind those ifs, because they get optimized away
1127 if (rows()!=cols())
1128 DUNE_THROW(FMatrixError, "There is no determinant for a " << rows() << "x" << cols() << " matrix!");
1129
1130 if (rows()==1)
1131 return (*this)[0][0];
1132
1133 if (rows()==2)
1134 return (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
1135
1136 if (rows()==3) {
1137 // code generated by maple
1138 field_type t4 = (*this)[0][0] * (*this)[1][1];
1139 field_type t6 = (*this)[0][0] * (*this)[1][2];
1140 field_type t8 = (*this)[0][1] * (*this)[1][0];
1141 field_type t10 = (*this)[0][2] * (*this)[1][0];
1142 field_type t12 = (*this)[0][1] * (*this)[2][0];
1143 field_type t14 = (*this)[0][2] * (*this)[2][0];
1144
1145 return (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
1146 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
1147
1148 }
1149
1150 AutonomousValue<MAT> A(asImp());
1151 field_type det;
1152 Simd::Mask<typename FieldTraits<value_type>::real_type>
1153 nonsingularLanes(true);
1154
1155 AutonomousValue<MAT>::luDecomposition(A, ElimDet(det), nonsingularLanes, false, doPivoting);
1156 det = Simd::cond(nonsingularLanes, det, field_type(0));
1157
1158 for (size_type i = 0; i < rows(); ++i)
1159 det *= A[i][i];
1160 return det;
1161 }
1162
1163#endif // DOXYGEN
1164
1165 namespace DenseMatrixHelp {
1166
1168 template <typename MAT, typename V1, typename V2>
1169 static inline void multAssign(const DenseMatrix<MAT> &matrix, const DenseVector<V1> & x, DenseVector<V2> & ret)
1170 {
1171 DUNE_ASSERT_BOUNDS(x.size() == matrix.cols());
1172 DUNE_ASSERT_BOUNDS(ret.size() == matrix.rows());
1173 typedef typename DenseMatrix<MAT>::size_type size_type;
1174
1175 for(size_type i=0; i<matrix.rows(); ++i)
1176 {
1177 ret[i] = 0.0;
1178 for(size_type j=0; j<matrix.cols(); ++j)
1179 {
1180 ret[i] += matrix[i][j]*x[j];
1181 }
1182 }
1183 }
1184
1185#if 0
1187 template <typename K, int rows, int cols>
1188 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
1189 {
1190 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
1191
1192 for(size_type i=0; i<cols(); ++i)
1193 {
1194 ret[i] = 0.0;
1195 for(size_type j=0; j<rows(); ++j)
1196 ret[i] += matrix[j][i]*x[j];
1197 }
1198 }
1199
1201 template <typename K, int rows, int cols>
1202 static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
1203 {
1204 FieldVector<K,rows> ret;
1205 multAssign(matrix,x,ret);
1206 return ret;
1207 }
1208
1210 template <typename K, int rows, int cols>
1211 static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x)
1212 {
1213 FieldVector<K,cols> ret;
1214 multAssignTransposed( matrix, x, ret );
1215 return ret;
1216 }
1217#endif
1218
1219 } // end namespace DenseMatrixHelp
1220
1222 template<typename MAT>
1223 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1224 {
1225 for (typename DenseMatrix<MAT>::size_type i=0; i<a.rows(); i++)
1226 s << a[i] << std::endl;
1227 return s;
1228 }
1229
1232} // end namespace Dune
1233
1234#endif
Various precision settings for calculations with FieldMatrix and FieldVector.
Implements a vector constructed from a given type representing a field and a compile-time given size.
Traits for type conversions and type information.
Implements a scalar vector view wrapper around an existing scalar.
A free function to provide the demangled class name of a given object or type as a string.
A few common exception classes.
Some useful basic math stuff.
Macro for wrapping boundary checks.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition boundschecking.hh:30
typename AutonomousValueType< T >::type AutonomousValue
Type free of internal references that T can be converted to.
Definition typetraits.hh:558
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition bigunsignedint.hh:278
#define DUNE_THROW(E, m)
Definition exceptions.hh:218
bool anyTrue(const Mask &mask)
Whether any entry is true
Definition simd/interface.hh:429
V cond(M &&mask, const V &ifTrue, const V &ifFalse)
Like the ?: operator.
Definition simd/interface.hh:386
bool allTrue(const Mask &mask)
Whether all entries are true
Definition simd/interface.hh:439
typename Overloads::RebindType< std::decay_t< S >, std::decay_t< V > >::type Rebind
Construct SIMD type with different scalar type.
Definition simd/interface.hh:253
constexpr std::size_t lanes()
Number of lanes in a SIMD type.
Definition simd/interface.hh:305
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition simd/interface.hh:324
Mask< V > mask(ADLTag< 0, std::is_same< V, Mask< V > >::value >, const V &v)
implements Simd::mask()
Definition defaults.hh:153
STL namespace.
Dune namespace.
Definition alignedallocator.hh:13
void swap(T &v1, T &v2, bool mask)
Definition simd.hh:472
int sign(const T &val)
Return the sign of the value.
Definition math.hh:180
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition math.hh:164
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition densematrix.hh:1169
A dense n x m matrix.
Definition densematrix.hh:140
ConstIterator const_iterator
typedef for stl compliant access
Definition densematrix.hh:244
derived_type operator-() const
Matrix negation.
Definition densematrix.hh:298
void solve(V1 &x, const V2 &b, bool doPivoting=true) const
Solve system A x = b.
void mv(const X &x, Y &y) const
y = A x
Definition densematrix.hh:368
Traits::value_type field_type
export the type representing the field
Definition densematrix.hh:160
derived_type & axpy(const field_type &a, const DenseMatrix< Other > &x)
vector space axpy operation (*this += a x)
Definition densematrix.hh:338
ConstIterator beforeEnd() const
Definition densematrix.hh:264
derived_type & operator=(const RHS &rhs)
Definition densematrix.hh:279
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition densematrix.hh:458
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition densematrix.hh:561
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition densematrix.hh:248
ConstIterator beforeBegin() const
Definition densematrix.hh:271
void invert(bool doPivoting=true)
Compute inverse.
static void luDecomposition(DenseMatrix< MAT > &A, Func func, Mask &nonsingularLanes, bool throwEarly, bool doPivoting)
do an LU-Decomposition on matrix A
Traits::value_type block_type
export the type representing the components
Definition densematrix.hh:163
void mtv(const X &x, Y &y) const
y = A^T x
Definition densematrix.hh:387
constexpr size_type cols() const
number of columns
Definition densematrix.hh:715
size_type size() const
size method (number of rows)
Definition densematrix.hh:200
constexpr size_type M() const
number of columns
Definition densematrix.hh:703
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition densematrix.hh:645
Iterator end()
end iterator
Definition densematrix.hh:222
Iterator beforeBegin()
Definition densematrix.hh:236
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition densematrix.hh:329
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition densematrix.hh:321
Iterator beforeEnd()
Definition densematrix.hh:229
Traits::value_type value_type
export the type representing the field
Definition densematrix.hh:157
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition densematrix.hh:484
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition densematrix.hh:512
void mmv(const X &x, Y &y) const
y -= A x
Definition densematrix.hh:445
constexpr size_type rows() const
number of rows
Definition densematrix.hh:709
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition densematrix.hh:627
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition densematrix.hh:498
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition densematrix.hh:312
bool operator!=(const DenseMatrix< Other > &x) const
Binary matrix incomparison.
Definition densematrix.hh:358
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition densematrix.hh:471
Traits::derived_type derived_type
type of derived matrix class
Definition densematrix.hh:154
row_reference operator[](size_type i)
random access
Definition densematrix.hh:189
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition densematrix.hh:723
Iterator RowIterator
rename the iterators for easier access
Definition densematrix.hh:211
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition densematrix.hh:178
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition densematrix.hh:528
void umv(const X &x, Y &y) const
y += A x
Definition densematrix.hh:406
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition densematrix.hh:242
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition densematrix.hh:546
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition densematrix.hh:169
constexpr size_type N() const
number of rows
Definition densematrix.hh:697
Traits::size_type size_type
The type used for the index access and size operation.
Definition densematrix.hh:166
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:175
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition densematrix.hh:536
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition densematrix.hh:213
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition densematrix.hh:172
bool operator==(const DenseMatrix< Other > &x) const
Binary matrix comparison.
Definition densematrix.hh:348
Iterator iterator
typedef for stl compliant access
Definition densematrix.hh:209
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition densematrix.hh:246
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition densematrix.hh:207
void umtv(const X &x, Y &y) const
y += A^T x
Definition densematrix.hh:419
ConstIterator begin() const
begin iterator
Definition densematrix.hh:251
field_type determinant(bool doPivoting=true) const
calculates the determinant of this matrix
Iterator begin()
begin iterator
Definition densematrix.hh:216
void umhv(const X &x, Y &y) const
y += A^H x
Definition densematrix.hh:432
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition densematrix.hh:289
ConstIterator end() const
end iterator
Definition densematrix.hh:257
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::real_type real_type
Definition densematrix.hh:34
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::field_type field_type
Definition densematrix.hh:33
A dense n x m matrix.
Definition fmatrix.hh:117
Base::size_type size_type
Definition fmatrix.hh:127
vector space out of a tensor product of fields.
Definition fvector.hh:95
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition densematrix.hh:59
Error thrown if operations of a FieldMatrix fail.
Definition densematrix.hh:126
Interface for a class of dense vectors over a given field.
Definition densevector.hh:229
size_type size() const
size method
Definition densevector.hh:336
Generic iterator class for dense vector and matrix implementations.
Definition densevector.hh:131
Default exception class for mathematical errors.
Definition exceptions.hh:241
Definition ftraits.hh:26
T field_type
export the type representing the field
Definition ftraits.hh:28
T real_type
export the type representing the real type of the field
Definition ftraits.hh:30
Definition matvectraits.hh:31
static ctype absolute_limit()
return threshold to declare matrix singular
Definition precision.hh:28
Include file for users of the SIMD abstraction layer.