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