dune-common  2.9.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 // 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 
16 #include <dune/common/classname.hh>
18 #include <dune/common/fvector.hh>
19 #include <dune/common/math.hh>
20 #include <dune/common/precision.hh>
21 #include <dune/common/simd/simd.hh>
24 
25 namespace Dune
26 {
27 
28  template<typename M> class DenseMatrix;
29 
30  template<typename M>
31  struct FieldTraits< DenseMatrix<M> >
32  {
35  };
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 
200  size_type size() const
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
Macro for wrapping boundary checks.
A free function to provide the demangled class name of a given object or type as a string.
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
Some useful basic math stuff.
Various precision settings for calculations with FieldMatrix and FieldVector.
Implements a scalar vector view wrapper around an existing scalar.
Traits for type conversions and type information.
#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
Rebind< bool, V > Mask
Mask type type of some SIMD type.
Definition: simd/interface.hh:289
Mask< V > mask(ADLTag< 0, std::is_same< V, Mask< V > >::value >, const V &v)
implements Simd::mask()
Definition: defaults.hh:153
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
auto max(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:425
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1169
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:829
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:838
static void multAssignTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x, FieldVector< K, cols > &ret)
calculates ret = matrix^T * x
Definition: fmatrix.hh:815
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Definition: fvector.hh:604
A dense n x m matrix.
Definition: densematrix.hh:140
derived_type & operator=(const RHS &rhs)
Definition: densematrix.hh:279
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
ConstIterator beforeEnd() const
Definition: densematrix.hh:264
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:458
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
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:536
constexpr size_type cols() const
number of columns
Definition: densematrix.hh:715
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:645
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
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:329
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition: densematrix.hh:312
Iterator end()
end iterator
Definition: densematrix.hh:222
Iterator beforeBegin()
Definition: densematrix.hh:236
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:561
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
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:528
constexpr size_type rows() const
number of rows
Definition: densematrix.hh:709
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: densematrix.hh:498
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
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition: densematrix.hh:289
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:211
constexpr static int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: densematrix.hh:178
derived_type & axpy(const field_type &a, const DenseMatrix< Other > &x)
vector space axpy operation (*this += a x)
Definition: densematrix.hh:338
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:321
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
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
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
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:627
bool exists([[maybe_unused]] size_type i, [[maybe_unused]] size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:723
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
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:546
ConstIterator end() const
end iterator
Definition: densematrix.hh:257
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::field_type field_type
Definition: densematrix.hh:33
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::real_type real_type
Definition: densematrix.hh:34
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.