diff --git a/common/fmatrix.hh b/common/fmatrix.hh
index 4c33a0edd189f7e83291426208d34213a145af09..51284af6af4181883c3bf5ff853c12c78a92e735 100644
--- a/common/fmatrix.hh
+++ b/common/fmatrix.hh
@@ -566,212 +566,72 @@ namespace Dune {
//===== iterator interface to rows of the matrix
-
- // forward declaration
- class ConstIterator;
-
- //! Iterator access to rows
- class Iterator
- {
- public:
- //! constructor
- Iterator (row_type* _p, int _i)
- {
- p = _p;
- i = _i;
- }
-
- //! empty constructor, use with care!
- Iterator ()
- { }
-
- //! prefix increment
- Iterator& operator++()
- {
- ++i;
- return *this;
- }
-
- //! prefix decrement
- Iterator& operator--()
- {
- --i;
- return *this;
- }
-
- //! equality
- bool operator== (const Iterator& it) const
- {
- return (p+i)==(it.p+it.i);
- }
-
- //! inequality
- bool operator!= (const Iterator& it) const
- {
- return (p+i)!=(it.p+it.i);
- }
-
- //! dereferencing
- row_type& operator* ()
- {
- return p[i];
- }
-
- //! arrow
- row_type* operator-> ()
- {
- return p+i;
- }
-
- //! return index
- int index ()
- {
- return i;
- }
-
- friend class ConstIterator;
-
- private:
- row_type* p;
- int i;
- };
+ //! Iterator class for sequential access
+ typedef Dune::GenericIterator<FieldMatrix<K,n,m>,row_type> Iterator;
+ //! typedef for stl compliant access
+ typedef Iterator iterator;
+ //! rename the iterators for easier access
+ typedef Iterator RowIterator;
+ //! rename the iterators for easier access
+ typedef typename row_type::Iterator ColIterator;
//! begin iterator
Iterator begin ()
{
- return Iterator(p,0);
+ return Iterator(*this,0);
}
//! end iterator
Iterator end ()
{
- return Iterator(p,n);
+ return Iterator(*this,n);
}
//! begin iterator
Iterator rbegin ()
{
- return Iterator(p,n-1);
+ return Iterator(*this,n-1);
}
//! end iterator
Iterator rend ()
{
- return Iterator(p,-1);
+ return Iterator(*this,-1);
}
+ //! Iterator class for sequential access
+ typedef Dune::GenericIterator<const FieldMatrix<K,n,m>,const row_type> ConstIterator;
+ //! typedef for stl compliant access
+ typedef ConstIterator const_iterator;
//! rename the iterators for easier access
- typedef Iterator RowIterator;
- typedef typename row_type::Iterator ColIterator;
-
-
- //! Iterator access to rows
- class ConstIterator
- {
- public:
- //! constructor
- ConstIterator (const row_type* _p, int _i) : p(_p), i(_i)
- { }
-
- //! empty constructor, use with care!
- ConstIterator ()
- {
- p = 0;
- i = 0;
- }
-
- //! prefix increment
- ConstIterator& operator++()
- {
- ++i;
- return *this;
- }
-
- //! prefix decrement
- ConstIterator& operator--()
- {
- --i;
- return *this;
- }
-
- //! equality
- bool operator== (const ConstIterator& it) const
- {
- return (p+i)==(it.p+it.i);
- }
-
- //! inequality
- bool operator!= (const ConstIterator& it) const
- {
- return (p+i)!=(it.p+it.i);
- }
-
- //! equality
- bool operator== (const Iterator& it) const
- {
- return (p+i)==(it.p+it.i);
- }
-
- //! inequality
- bool operator!= (const Iterator& it) const
- {
- return (p+i)!=(it.p+it.i);
- }
-
- //! dereferencing
- const row_type& operator* () const
- {
- return p[i];
- }
-
- //! arrow
- const row_type* operator-> () const
- {
- return p+i;
- }
-
- //! return index
- int index () const
- {
- return i;
- }
-
- friend class Iterator;
-
- private:
- const row_type* p;
- int i;
- };
+ typedef ConstIterator ConstRowIterator;
+ //! rename the iterators for easier access
+ typedef typename row_type::ConstIterator ConstColIterator;
//! begin iterator
ConstIterator begin () const
{
- return ConstIterator(p,0);
+ return ConstIterator(*this,0);
}
//! end iterator
ConstIterator end () const
{
- return ConstIterator(p,n);
+ return ConstIterator(*this,n);
}
//! begin iterator
ConstIterator rbegin () const
{
- return ConstIterator(p,n-1);
+ return ConstIterator(*this,n-1);
}
//! end iterator
ConstIterator rend () const
{
- return ConstIterator(p,-1);
+ return ConstIterator(*this,-1);
}
- //! rename the iterators for easier access
- typedef ConstIterator ConstRowIterator;
- typedef typename row_type::ConstIterator ConstColIterator;
-
-
//===== assignment from scalar
FieldMatrix& operator= (const K& k)
{
@@ -1126,7 +986,7 @@ namespace Dune {
return HelpMat::determinantMatrix(*this);
}
-
+#ifdef USE_DEPRECATED_K1
/** \brief Special type for 1x1 matrices
*/
template<class K>
@@ -1363,7 +1223,336 @@ namespace Dune {
// the data, just a single scalar
K a;
};
+#endif // USE_DEPRECATED_K1
+
+ /** \brief Special type for 1x1 matrices
+ */
+ template<class K>
+ class FieldMatrix<K,1,1>
+ {
+ public:
+ // standard constructor and everything is sufficient ...
+
+ //===== type definitions and constants
+
+ //! export the type representing the field
+ typedef K field_type;
+
+ //! export the type representing the components
+ typedef K block_type;
+
+ //! We are at the leaf of the block recursion
+ enum {
+ //! The number of block levels we contain.
+ //! This is always one for this type.
+ blocklevel = 1
+ };
+
+ //! Each row is implemented by a field vector
+ typedef FieldVector<K,1> row_type;
+
+ //! export size
+ enum {
+ //! \brief The number of rows.
+ //! This is always one for this type.
+ rows = 1,
+ n = 1,
+ //! \brief The number of columns.
+ //! This is always one for this type.
+ cols = 1,
+ m = 1
+ };
+
+ //===== random access interface to rows of the matrix
+
+ //! random access to the rows
+ row_type& operator[] (int i)
+ {
+#ifdef DUNE_FMatrix_WITH_CHECKING
+ if (i<0 || i>=n) DUNE_THROW(FMatrixError,"index out of range");
+#endif
+ return a;
+ }
+ //! same for read only access
+ const row_type& operator[] (int i) const
+ {
+#ifdef DUNE_FMatrix_WITH_CHECKING
+ if (i<0 || i>=n) DUNE_THROW(FMatrixError,"index out of range");
+#endif
+ return a;
+ }
+
+ //===== iterator interface to rows of the matrix
+ //! Iterator class for sequential access
+ typedef Dune::GenericIterator<FieldMatrix<K,n,m>,row_type> Iterator;
+ //! typedef for stl compliant access
+ typedef Iterator iterator;
+ //! rename the iterators for easier access
+ typedef Iterator RowIterator;
+ //! rename the iterators for easier access
+ typedef typename row_type::Iterator ColIterator;
+
+ //! begin iterator
+ Iterator begin ()
+ {
+ return Iterator(*this,0);
+ }
+
+ //! end iterator
+ Iterator end ()
+ {
+ return Iterator(*this,n);
+ }
+
+ //! begin iterator
+ Iterator rbegin ()
+ {
+ return Iterator(*this,n-1);
+ }
+
+ //! end iterator
+ Iterator rend ()
+ {
+ return Iterator(*this,-1);
+ }
+
+ //! Iterator class for sequential access
+ typedef Dune::GenericIterator<const FieldMatrix<K,n,m>,const row_type> ConstIterator;
+ //! typedef for stl compliant access
+ typedef ConstIterator const_iterator;
+ //! rename the iterators for easier access
+ typedef ConstIterator ConstRowIterator;
+ //! rename the iterators for easier access
+ typedef typename row_type::ConstIterator ConstColIterator;
+
+ //! begin iterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(*this,0);
+ }
+
+ //! end iterator
+ ConstIterator end () const
+ {
+ return ConstIterator(*this,n);
+ }
+
+ //! begin iterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(*this,n-1);
+ }
+
+ //! end iterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(*this,-1);
+ }
+
+ //===== assignment from scalar
+
+ FieldMatrix& operator= (const K& k)
+ {
+ a[0] = k;
+ return *this;
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ FieldMatrix& operator+= (const K& y)
+ {
+ a[0] += y;
+ return *this;
+ }
+
+ //! vector space subtraction
+ FieldMatrix& operator-= (const K& y)
+ {
+ a[0] -= y;
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ FieldMatrix& operator*= (const K& k)
+ {
+ a[0] *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ FieldMatrix& operator/= (const K& k)
+ {
+ a[0] /= k;
+ return *this;
+ }
+
+ //===== linear maps
+
+ //! y += A x
+ void umv (const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p += a[0] * x.p;
+ }
+
+ //! y += A^T x
+ void umtv (const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p += a[0] * x.p;
+ }
+
+ //! y += A^H x
+ void umhv (const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p += fm_ck(a[0]) * x.p;
+ }
+
+ //! y -= A x
+ void mmv (const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p -= a[0] * x.p;
+ }
+
+ //! y -= A^T x
+ void mmtv (const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p -= a[0] * x.p;
+ }
+
+ //! y -= A^H x
+ void mmhv (const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p -= fm_ck(a[0]) * x.p;
+ }
+
+ //! y += alpha A x
+ void usmv (const K& alpha, const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p += alpha * a[0] * x.p;
+ }
+
+ //! y += alpha A^T x
+ void usmtv (const K& alpha, const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p += alpha * a[0] * x.p;
+ }
+
+ //! y += alpha A^H x
+ void usmhv (const K& alpha, const FieldVector<K,1>& x, FieldVector<K,1>& y) const
+ {
+ y.p += alpha * fm_ck(a[0]) * x.p;
+ }
+
+ //===== norms
+
+ //! frobenius norm: sqrt(sum over squared values of entries)
+ double frobenius_norm () const
+ {
+ return sqrt(fvmeta_abs2(a[0]));
+ }
+
+ //! square of frobenius norm, need for block recursion
+ double frobenius_norm2 () const
+ {
+ return fvmeta_abs2(a[0]);
+ }
+
+ //! infinity norm (row sum norm, how to generalize for blocks?)
+ double infinity_norm () const
+ {
+ return fvmeta_abs(a[0]);
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ return fvmeta_abs_real(a[0]);
+ }
+
+ //===== solve
+
+ //! Solve system A x = b
+ void solve (FieldVector<K,1>& x, const FieldVector<K,1>& b) const
+ {
+ x.p = b.p/a[0];
+ }
+
+ //! compute inverse
+ void invert ()
+ {
+ a[0] = 1/a[0];
+ }
+
+ //! left multiplication
+ FieldMatrix& leftmultiply (const FieldMatrix& M)
+ {
+ a[0] *= M.a[0];
+ return *this;
+ }
+
+ //! left multiplication
+ FieldMatrix& rightmultiply (const FieldMatrix& M)
+ {
+ a[0] *= M.a[0];
+ return *this;
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in row direction
+ int N () const
+ {
+ return 1;
+ }
+
+ //! number of blocks in column direction
+ int M () const
+ {
+ return 1;
+ }
+
+ //! row dimension of block r
+ int rowdim (int r) const
+ {
+ return 1;
+ }
+
+ //! col dimension of block c
+ int coldim (int c) const
+ {
+ return 1;
+ }
+
+ //! dimension of the destination vector space
+ int rowdim () const
+ {
+ return 1;
+ }
+
+ //! dimension of the source vector space
+ int coldim () const
+ {
+ return 1;
+ }
+
+ //===== query
+
+ //! return true when (i,j) is in pattern
+ bool exists (int i, int j) const
+ {
+ return i==0 && j==0;
+ }
+
+ //===== conversion operator
+
+ operator K () const {return a[0];}
+
+ private:
+ // the data, just a single row with a single scalar
+ row_type a;
+ };
namespace FMatrixHelp {
diff --git a/common/fvector.hh b/common/fvector.hh
index a0310fa79d8774b0bae964ced12ea232c2a7bea4..3151cb47f2d1e9c585cbd03b249114eab37f8d7d 100644
--- a/common/fvector.hh
+++ b/common/fvector.hh
@@ -4,10 +4,14 @@
#ifndef DUNE_FVECTOR_HH
#define DUNE_FVECTOR_HH
+// this for backwards compatiblity
+#define USE_DEPRECATED_K1
+
#include <math.h>
#include <complex>
#include "exceptions.hh"
+#include "genericiterator.hh"
namespace Dune {
@@ -405,221 +409,80 @@ namespace Dune {
return p[i];
}
- // forward declaration
- class ConstIterator;
-
//! Iterator class for sequential access
- class Iterator
- {
- public:
- //! constructor
- Iterator (K* _p, int _i)
- {
- p = _p;
- i = _i;
- }
-
- //! empty constructor, use with care
- Iterator ()
- { }
-
- //! prefix increment
- Iterator& operator++()
- {
- ++i;
- return *this;
- }
-
- //! prefix decrement
- Iterator& operator--()
- {
- --i;
- return *this;
- }
-
- //! equality
- bool operator== (const Iterator& it) const
- {
- return (p+i)==(it.p+it.i);
- }
-
- //! inequality
- bool operator!= (const Iterator& it) const
- {
- return (p+i)!=(it.p+it.i);
- }
-
- //! dereferencing
- K& operator* ()
- {
- return p[i];
- }
-
- //! arrow
- K* operator-> ()
- {
- return p+i;
- }
-
- //! return index
- int index ()
- {
- return i;
- }
-
- friend class ConstIterator;
-
- private:
- K* p;
- int i;
- };
+ typedef Dune::GenericIterator<FieldVector<K,n>,K> Iterator;
+ //! typedef for stl compliant access
+ typedef Iterator iterator;
//! begin iterator
Iterator begin ()
{
- return Iterator(p,0);
+ return Iterator(*this,0);
}
//! end iterator
Iterator end ()
{
- return Iterator(p,n);
+ return Iterator(*this,n);
}
//! begin iterator
Iterator rbegin ()
{
- return Iterator(p,n-1);
+ return Iterator(*this,n-1);
}
//! end iterator
Iterator rend ()
{
- return Iterator(p,-1);
+ return Iterator(*this,-1);
}
//! return iterator to given element or end()
Iterator find (int i)
{
if (i>=0 && i<n)
- return Iterator(p,i);
+ return Iterator(*this,i);
else
- return Iterator(p,n);
+ return Iterator(*this,n);
}
//! ConstIterator class for sequential access
- class ConstIterator
- {
- public:
- //! constructor
- ConstIterator ()
- {
- p = 0;
- i = 0;
- }
-
- //! constructor from pointer
- ConstIterator (const K* _p, int _i) : p(_p), i(_i)
- { }
-
- //! constructor from non_const iterator
- ConstIterator (const Iterator& it) : p(it.p), i(it.i)
- { }
-
- //! prefix increment
- ConstIterator& operator++()
- {
- ++i;
- return *this;
- }
-
- //! prefix decrement
- ConstIterator& operator--()
- {
- --i;
- return *this;
- }
-
- //! equality
- bool operator== (const ConstIterator& it) const
- {
- return (p+i)==(it.p+it.i);
- }
-
- //! inequality
- bool operator!= (const ConstIterator& it) const
- {
- return (p+i)!=(it.p+it.i);
- }
-
- //! equality
- bool operator== (const Iterator& it) const
- {
- return (p+i)==(it.p+it.i);
- }
-
- //! inequality
- bool operator!= (const Iterator& it) const
- {
- return (p+i)!=(it.p+it.i);
- }
-
- //! dereferencing
- const K& operator* () const
- {
- return p[i];
- }
-
- //! arrow
- const K* operator-> () const
- {
- return p+i;
- }
-
- //! return index corresponding to pointer
- int index () const
- {
- return i;
- }
-
- friend class Iterator;
-
- private:
- const K* p;
- int i;
- };
+ typedef Dune::GenericIterator<const FieldVector<K,n>,const K> ConstIterator;
+ //! typedef for stl compliant access
+ typedef ConstIterator const_iterator;
//! begin ConstIterator
ConstIterator begin () const
{
- return ConstIterator(p,0);
+ return ConstIterator(*this,0);
}
//! end ConstIterator
ConstIterator end () const
{
- return ConstIterator(p,n);
+ return ConstIterator(*this,n);
}
//! begin ConstIterator
ConstIterator rbegin () const
{
- return ConstIterator(p,n-1);
+ return ConstIterator(*this,n-1);
}
//! end ConstIterator
ConstIterator rend () const
{
- return ConstIterator(p,-1);
+ return ConstIterator(*this,-1);
}
//! return iterator to given element or end()
ConstIterator find (int i) const
{
if (i>=0 && i<n)
- return ConstIterator(p,i);
+ return ConstIterator(*this,i);
else
- return ConstIterator(p,n);
+ return ConstIterator(*this,n);
}
//===== vector space arithmetic
@@ -779,9 +642,9 @@ namespace Dune {
return s;
}
+#ifdef USE_DEPRECATED_K1
// forward declarations
template<class K> class K1Vector;
- template<class K> class K11Matrix;
/**! \brief Vectors containing only one component
*/
@@ -958,7 +821,279 @@ namespace Dune {
// the data
K p;
};
+#endif // USE_DEPRECATED_K1
+
+ // forward declarations
+ template<class K, int n, int m> class FieldMatrix;
+
+ /**! Vectors containing only one component
+ */
+ template<class K>
+ class FieldVector<K,1>
+ {
+ enum { n = 1 };
+ public:
+ friend class FieldMatrix<K,1,1>;
+
+ //===== type definitions and constants
+
+ //! export the type representing the field
+ typedef K field_type;
+
+ //! export the type representing the components
+ typedef K block_type;
+
+ //! We are at the leaf of the block recursion
+ enum {blocklevel = 1};
+
+ //! export size
+ enum {size = 1};
+
+ //===== construction
+
+ /** \brief Default constructor */
+ FieldVector ()
+ { }
+
+ /** \brief Constructor with a given scalar */
+ FieldVector (const K& k)
+ {
+ p = k;
+ }
+
+ /** \brief Assignment from the base type */
+ FieldVector& operator= (const K& k)
+ {
+ p = k;
+ return *this;
+ }
+
+ //===== access to components
+
+ //! random access
+ K& operator[] (int i)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i != 0) DUNE_THROW(MathError,"index out of range");
+#endif
+ return p;
+ }
+
+ //! same for read only access
+ const K& operator[] (int i) const
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i != 0) DUNE_THROW(MathError,"index out of range");
+#endif
+ return p;
+ }
+
+ //! Iterator class for sequential access
+ typedef Dune::GenericIterator<FieldVector<K,n>,K> Iterator;
+ //! typedef for stl compliant access
+ typedef Iterator iterator;
+
+ //! begin iterator
+ Iterator begin ()
+ {
+ return Iterator(*this,0);
+ }
+
+ //! end iterator
+ Iterator end ()
+ {
+ return Iterator(*this,n);
+ }
+
+ //! begin iterator
+ Iterator rbegin ()
+ {
+ return Iterator(*this,n-1);
+ }
+
+ //! end iterator
+ Iterator rend ()
+ {
+ return Iterator(*this,-1);
+ }
+
+ //! return iterator to given element or end()
+ Iterator find (int i)
+ {
+ if (i>=0 && i<n)
+ return Iterator(*this,i);
+ else
+ return Iterator(*this,n);
+ }
+
+ //! ConstIterator class for sequential access
+ typedef Dune::GenericIterator<const FieldVector<K,n>,const K> ConstIterator;
+ //! typedef for stl compliant access
+ typedef ConstIterator const_iterator;
+
+ //! begin ConstIterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(*this,0);
+ }
+
+ //! end ConstIterator
+ ConstIterator end () const
+ {
+ return ConstIterator(*this,n);
+ }
+
+ //! begin ConstIterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(*this,n-1);
+ }
+
+ //! end ConstIterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(*this,-1);
+ }
+
+ //! return iterator to given element or end()
+ ConstIterator find (int i) const
+ {
+ if (i>=0 && i<n)
+ return ConstIterator(*this,i);
+ else
+ return ConstIterator(*this,n);
+ }
+ //===== vector space arithmetic
+
+ //! vector space addition
+ FieldVector& operator+= (const FieldVector& y)
+ {
+ p += y.p;
+ return *this;
+ }
+
+ //! vector space subtraction
+ FieldVector& operator-= (const FieldVector& y)
+ {
+ p -= y.p;
+ return *this;
+ }
+ //! vector space add scalar to each comp
+ FieldVector& operator+= (const K& k)
+ {
+ p += k;
+ return *this;
+ }
+
+ //! vector space subtract scalar from each comp
+ FieldVector& operator-= (const K& k)
+ {
+ p -= k;
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ FieldVector& operator*= (const K& k)
+ {
+ p *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ FieldVector& operator/= (const K& k)
+ {
+ p /= k;
+ return *this;
+ }
+
+ //! vector space axpy operation
+ FieldVector& axpy (const K& a, const FieldVector& y)
+ {
+ p += a*y.p;
+ return *this;
+ }
+
+
+ //===== Euclidean scalar product
+
+ //! scalar product
+ const K operator* (const FieldVector& y) const
+ {
+ return p*y.p;
+ }
+
+
+ //===== norms
+
+ //! one norm (sum over absolute values of entries)
+ double one_norm () const
+ {
+ return fvmeta_abs(p);
+ }
+
+ //! simplified one norm (uses Manhattan norm for complex values)
+ double one_norm_real () const
+ {
+ return fvmeta_abs_real(p);
+ }
+
+ //! two norm sqrt(sum over squared values of entries)
+ double two_norm () const
+ {
+ return sqrt(fvmeta_abs2(p));
+ }
+
+ //! square of two norm (sum over squared values of entries), need for block recursion
+ double two_norm2 () const
+ {
+ return fvmeta_abs2(p);
+ }
+
+ //! infinity norm (maximum of absolute values of entries)
+ double infinity_norm () const
+ {
+ return fvmeta_abs(p);
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ return fvmeta_abs_real(p);
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of size 1 here)
+ int N () const
+ {
+ return 1;
+ }
+
+ //! dimension of the vector space (==1)
+ int dim () const
+ {
+ return 1;
+ }
+
+ //! Send vector to output stream
+ void print (std::ostream& s) const
+ {
+ s << p;
+ }
+ //===== conversion operator
+
+ /** \brief Conversion operator */
+ operator K () {return p;}
+
+ /** \brief Const conversion operator */
+ operator K () const {return p;}
+
+ private:
+ // the data
+ K p;
+ };
/** @} end documentation */
diff --git a/common/test/.gitignore b/common/test/.gitignore
index e3a967229c7c7d063670033e1196b491e5fc523a..bccf541a3937a3e9eb433b5540e1095d6a3eb717 100644
--- a/common/test/.gitignore
+++ b/common/test/.gitignore
@@ -11,6 +11,7 @@ iteratorfacadetest
sllisttest
tuplestest
settest
+fmatrixtest
*.gcda
*.gcno
gmon.out
\ No newline at end of file
diff --git a/common/test/Makefile.am b/common/test/Makefile.am
index 507d4daf1aa357174540c917e5c74abdcbd3d3fc..7cfdb95320fcf750b4198b782a354c95ee0c4cc6 100644
--- a/common/test/Makefile.am
+++ b/common/test/Makefile.am
@@ -1,7 +1,7 @@
# $Id$
TESTPROGS = parsetest test-stack arraylisttest smartpointertest \
- sllisttest iteratorfacadetest tuplestest
+ sllisttest iteratorfacadetest tuplestest fmatrixtest
# which tests to run
TESTS = $(TESTPROGS)
@@ -29,3 +29,5 @@ tuplestest_SOURCES = tuplestest.cc
# mention headers so that they are distributed too
iteratorfacadetest_SOURCES = iteratorfacadetest.cc iteratorfacadetest.hh \
iteratortest.hh
+
+fmatrixtest_SOURCES = fmatrixtest.cc
diff --git a/common/test/fmatrixtest.cc b/common/test/fmatrixtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c27779d3e666acb3f789740d64d3dee481aa767b
--- /dev/null
+++ b/common/test/fmatrixtest.cc
@@ -0,0 +1,76 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#define DUNE_ISTL_WITH_CHECKING
+#include "../fmatrix.hh"
+#include <iostream>
+
+using namespace Dune;
+
+template<class K, int n, int m>
+void test_matrix()
+{
+ FieldMatrix<K,n,m> A;
+ FieldVector<K,n> f;
+ FieldVector<K,m> v;
+
+ // assign matrix
+ A=0;
+ // random access matrix
+ for (int i=0; i<A.rowdim(); i++)
+ for (int j=0; j<A.coldim(); j++)
+ A[i][j] = i*j;
+ // iterator matrix
+ typename FieldMatrix<K,n,m>::RowIterator rit = A.begin();
+ for (; rit!=A.end(); ++rit)
+ {
+ typename FieldMatrix<K,n,m>::ColIterator cit = rit->begin();
+ for (; cit!=rit->end(); ++cit)
+ (*cit) *= 2;
+ }
+
+ // assign vector
+ f = 1;
+ // random access vector
+ for (int i=0; i<v.dim(); i++)
+ v[i] = i;
+ // iterator vector
+ typename FieldVector<K,m>::iterator it = v.begin();
+ typename FieldVector<K,m>::ConstIterator end = v.end();
+ for (; it!=end; ++it)
+ (*it) *= 2;
+ // reverse iterator vector
+ it = v.rbegin();
+ end = v.rend();
+ for (; it!=end; --it)
+ (*it) /= 2;
+ // find vector
+ for (int i=0; i<v.dim(); i++)
+ {
+ it = v.find(i);
+ (*it) += 1;
+ }
+
+ // matrix vector product
+ A.umv(v,f);
+
+ A.infinity_norm();
+
+ // print matrix
+ std::cout << A << std::endl;
+ // print vector
+ std::cout << f << std::endl;
+}
+
+int main()
+{
+ try {
+ test_matrix<float, 1, 1>();
+ test_matrix<double, 1, 1>();
+ test_matrix<int, 10, 5>();
+ test_matrix<double, 5, 10>();
+ }
+ catch (Dune::Exception & e)
+ {
+ std::cerr << "Exception: " << e << std::endl;
+ }
+}