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 
14 #include <dune/common/classname.hh>
16 #include <dune/common/fvector.hh>
17 #include <dune/common/math.hh>
18 #include <dune/common/precision.hh>
19 #include <dune/common/simd.hh>
21 #include <dune/common/unused.hh>
22 
23 namespace 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:
209  //===== access to components
210 
213  {
214  return asImp().mat_access(i);
215  }
216 
218  {
219  return asImp().mat_access(i);
220  }
221 
223  size_type size() const
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 
698  size_type rows() const
699  {
700  return asImp().mat_rows();
701  }
702 
704  size_type cols() const
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 
804  template<class Func, class Mask>
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>
851  template<typename Func, class Mask>
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());
989  SimdMask<typename FieldTraits<value_type>::real_type>
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());
1071  SimdMask<typename FieldTraits<value_type>::real_type>
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;
1144  SimdMask<typename FieldTraits<value_type>::real_type>
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
derived_type & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition: densematrix.hh:322
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
derived_type & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition: densematrix.hh:312
derived_type & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition: densematrix.hh:348
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
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:395
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:523
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:634
size_type size() const
size method (number of rows)
Definition: densematrix.hh:223
size_type M() const
number of columns
Definition: densematrix.hh:692
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:339
Iterator end()
end iterator
Definition: densematrix.hh:245
Iterator beforeBegin()
Definition: densematrix.hh:259
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:548
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.
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:515
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
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:331
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
size_type N() const
number of rows
Definition: densematrix.hh:686
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
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
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:616
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
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:533
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)
masked Simd assignment (scalar version)
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
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (Mar 28, 23:30, 2024)