* Expression Templates for Vectors (default: disabled)

* Test program
* experiments for Expression Templates for Matrizes

[[Imported from SVN: r2438]]

common/exprtmpl.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_EXPRTMPL_HH
+#define DUNE_EXPRTMPL_HH
+
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+#include <cmath>
+#include "stdstreams.hh"
+
+struct Indent
+{
+  int i;
+  Indent() {
+    i = 0;
+  }
+  void operator ++ ()
+  {
+    i += 3;
+  };
+  void operator -- ()
+  {
+    i -= 3;
+  };
+};
+
+extern Indent INDENT;
+
+std::ostream & operator << (std::ostream & s, const Indent & i)
+{
+  for (int n = 0; n < i.i; n++) s << " ";
+  return s;
+}
+
+namespace Dune {
+
+  template<class V> class FlatIterator;
+  template<class K, int N> class FieldVector;
+
+  /**
+     Type Traits for field_type and block_type
+   */
+  template<class T>
+  struct BlockType
+  {};
+  template<class T>
+  struct FieldType
+  {};
+  // FieldType specializations for basic data types
+  template<>
+  struct FieldType<double>
+  {
+    typedef double type;
+  };
+  template<>
+  struct FieldType<float>
+  {
+    typedef float type;
+  };
+  template<>
+  struct FieldType<int>
+  {
+    typedef int type;
+  };
+  // FieldType specialization for FieldMatrix
+  template <class K, int N>
+  struct FieldType< FieldVector<K,N> >
+  {
+    typedef K type;
+  };
+  // BlockType specialization for FieldMatrix
+  template <class K, int N>
+  struct BlockType< FieldVector<K,N> >
+  {
+    typedef K type;
+  };
+  // FieldType specialization for const T
+  template<class T>
+  struct BlockType<const T>
+  {
+    typedef const typename BlockType<T>::type type;
+  };
+  // BlockType specialization for const T
+  template<class T>
+  struct FieldType<const T>
+  {
+    typedef const typename FieldType<T>::type type;
+  };
+
+  namespace ExprTmpl {
+
+    /* Important Classes */
+    template <class V> class ConstRef;
+    template <class Ex> class Expression;
+    template <class A, class B, template<class> class Op> class ExBinOp;
+    template <class I> class Vector;
+
+    /**
+       Type Trait for nested Expression
+     */
+    template<class T>
+    struct BlockExpression
+    {};
+    /**
+       Type Trait for Implementation of an Expression
+     */
+    template<class T>
+    struct ExpressionImp
+    {};
+
+    template <class V>
+    struct ExpressionImp< ExprTmpl::ConstRef<V> >
+    {
+      typedef ExprTmpl::ConstRef<V> type;
+    };
+
+    // TypeTraits
+    template <>
+    struct ExpressionImp<double>
+    {
+      typedef double type;
+    };
+    template <>
+    struct ExpressionImp<float>
+    {
+      typedef float type;
+    };
+    template <>
+    struct ExpressionImp<int>
+    {
+      typedef int type;
+    };
+    template <class Ex>
+    struct ExpressionImp< Expression<Ex> >
+    {
+      typedef Ex type;
+    };
+
+    /** Identify End of Expression Recursion */
+    template <class V>
+    struct isEndOfExpressionRecusion
+    {
+      enum { value=false };
+    };
+
+    template <>
+    struct isEndOfExpressionRecusion<double>
+    {
+      enum { value=true };
+    };
+
+    template <>
+    struct isEndOfExpressionRecusion<int>
+    {
+      enum { value=true };
+    };
+
+    template <>
+    struct isEndOfExpressionRecusion<float>
+    {
+      enum { value=true };
+    };
+
+    /** Nested Expression in a ConstRef<V> */
+    template <class V>
+    struct BlockExpression< ConstRef<V> >
+    {
+      typedef ExprTmpl::Expression<
+          ExprTmpl::ConstRef<typename BlockType<V>::type> > type;
+    };
+
+    /** No nested Expression in ConstRef<FieldVector> */
+    template<class K, int N>
+    struct BlockExpression< ConstRef< FieldVector<K,N> > >
+    {
+      typedef K type;
+    };
+
+    /**
+       Wrapper for an expression template
+     */
+    template <class Ex>
+    class Expression
+    {
+    public:
+      Expression(const Ex & x) : ex(x) {}
+      typedef typename FieldType<Ex>::type field_type;
+      typedef typename BlockExpression<Ex>::type BlockExpr;
+      BlockExpr operator[] ( int i ) const {
+        return ex[i];
+      }
+      int N() const { return ex.N(); }
+    private:
+      Ex ex;
+    };
+
+    /**
+       1 Dimensional Vector Base Class
+       for Glommable Expr1 Templates
+     */
+    template <class I>
+    class Vector {
+    public:
+      explicit Vector() {}
+      typedef typename BlockType<I>::type block_type;
+      typedef typename FieldType<I>::type field_type;
+      //! dimension of the vector space
+      int N() const {
+        return asImp().N();
+      }
+      block_type & operator[] (int i) {
+        return asImp()[i];
+      }
+      const block_type & operator[] (int i) const {
+        return asImp()[i];
+      }
+      //! begin for FlatIterator
+      FlatIterator<I> fbegin() { return FlatIterator<I>(asImp().begin()); }
+      //! end for FlatIterator
+      FlatIterator<I> fend() { return FlatIterator<I>(asImp().end()); }
+      //! begin for ConstFlatIterator
+      FlatIterator<const I> fbegin() const {
+        return FlatIterator<const I>(asImp().begin());
+      }
+      //! end for ConstFlatIterator
+      FlatIterator<const I> fend() const {
+        return FlatIterator<const I>(asImp().end());
+      }
+      //! assign Vector from Expression
+      template <class E> I& assignFrom(Expression<E>& x) {
+#warning should there be a resize?
+#ifdef DUNE_ISTL_WITH_CHECKING
+        assert(N() == x.N());
+#endif
+        Dune::dvverb << INDENT << "Assign Vector from Expression\n";
+        ++INDENT;
+        for (int i=0; i<N(); ++i) { asImp()[i] = x[i]; }
+        --INDENT;
+        return asImp();
+      }
+      //! assign Vector from Vector
+      template <class V> I& assignFrom(const Vector<V>& v) {
+#ifdef DUNE_ISTL_WITH_CHECKING
+        assert(N() == v.N());
+#endif
+        Dune::dvverb << INDENT << "Assign Vector from Vector\n";
+        ++INDENT;
+        for (int i=0; i<N(); ++i) { asImp()[i] = v[i]; }
+        --INDENT;
+        return asImp();
+      }
+      /*
+         I& assignFrom(const Vector<block_type> & x) {
+         Dune::dvverb << INDENT << "Assign Vector block_type\n";
+         ++INDENT;
+         for (int i=0; i < asImp().N(); i++) asImp()[i] = x;
+         --INDENT;
+         return asImp();
+         }
+       */
+      I& assignFrom(field_type x) {
+        Dune::dvverb << INDENT << "Assign Vector from field_type\n";
+        ++INDENT;
+        for (int i=0; i<N(); ++i) { asImp()[i] = x; }
+        --INDENT;
+        return asImp();
+      }
+#if 0
+#warning rewrite these!
+      template <class E> Vector<I>& operator+=(const Expression<E>& x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] += x(i);
+        return asImp();
+      }
+      template <class V> Vector<I>& operator+=(const Vector<V>& x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] += x(i);
+        return asImp();
+      }
+      template <class E> Vector<I>& operator-=(const Expression<E>& x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] -= x(i);
+        return asImp();
+      }
+      template <class V> Vector<I>& operator-=(const Vector<V>& x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] -= x(i);
+        return asImp();
+      }
+#endif
+      Vector<I>& operator+=(field_type x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] += x;
+        return asImp();
+      }
+      Vector<I>& operator-=(field_type x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] -= x;
+        return asImp();
+      }
+      Vector<I>& operator*=(field_type x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] *= x;
+        return asImp();
+      }
+      Vector<I>& operator/=(field_type x) {
+        for (int i=0; i < asImp().N(); i++) asImp()[i] /= x;
+        return asImp();
+      }
+    private:
+      I & asImp() { return static_cast<I&>(*this); }
+      const I & asImp() const { return static_cast<const I&>(*this); }
+    };
+
+    template <class V>
+    class ConstRef
+    {
+    public:
+      typedef typename FieldType<V>::type field_type;
+      typedef typename BlockType<V>::type block_type;
+      typedef typename BlockExpression< ConstRef<V> >::type BlockExpr;
+      typedef typename ExpressionImp<BlockExpr>::type BlockExprImp;
+      ConstRef (const Vector<V> & _v) : v(_v) {}
+      BlockExpr operator[] (int i) const {
+        Dune::dvverb << INDENT << "ConstRef->dereference " << v[i] << std::endl;
+        return BlockExprImp(v[i]);
+      }
+      int N() const { return v.N(); };
+    private:
+      const Vector<V> & v;
+    };
+
+  } // namespace ExprTmpl
+
+  template <class A>
+  struct FieldType< ExprTmpl::Expression<A> >
+  {
+    typedef typename FieldType<A>::type type;
+  };
+
+  template <class V>
+  struct FieldType< ExprTmpl::ConstRef<V> >
+  {
+    typedef typename FieldType<V>::type type;
+  };
+
+  template <class I>
+  struct FieldType< ExprTmpl::Vector<I> >
+  {
+    typedef typename FieldType<I>::type type;
+  };
+
+  // ----
+  template <class A>
+  struct BlockType< ExprTmpl::Expression<A> >
+  {
+    typedef typename BlockType<A>::type type;
+  };
+
+  template <class V>
+  struct BlockType< ExprTmpl::ConstRef<V> >
+  {
+    typedef typename BlockType<V>::type type;
+  };
+
+  template <class I>
+  struct BlockType< ExprTmpl::Vector<I> >
+  {
+    typedef typename BlockType<I>::type type;
+  };
+
+  // ---
+  namespace ExprTmpl {
+
+    // ---
+
+  }
+
+  // OPERATORS
+
+  /* Scalar-Expression Operator */
+#define OP *
+#define ExpressionOpScalar ExpressionMulScalar
+#define ScalarOpExpression ScalarMulExpression
+#include "exprtmpl/scalar.inc"
+
+#define OP /
+#define ExpressionOpScalar ExpressionDivScalar
+#include "exprtmpl/scalar.inc"
+
+#define OP +
+#define ExpressionOpScalar ExpressionAddScalar
+#define ScalarOpExpression ScalarAdExpression
+#include "exprtmpl/scalar.inc"
+
+#define OP -
+#define ExpressionOpScalar ExpressionMinScalar
+#define ScalarOpExpression ScalarMinExpression
+#include "exprtmpl/scalar.inc"
+
+  /* Expression-Expression Operator */
+#define OP +
+#define ExpressionOpExpression ExpressionAddExpression
+#include "exprtmpl/exprexpr.inc"
+
+#define OP -
+#define ExpressionOpExpression ExpressionMinExpression
+#include "exprtmpl/exprexpr.inc"
+
+} // namespace Dune
+
+#endif // DUNE_EXPRTMPL_HH

common/exprtmpl/exprexpr.inc

+// -*- C++ -*-
+#ifndef OP
+#error OP undefined
+#endif
+
+namespace ExprTmpl {
+
+template <class A, class B>
+class ExpressionOpExpression
+{
+public:
+//#warning compare field_type of class A and B and choose infomation richest one!
+  typedef typename FieldType<A>::type field_type;
+  typedef typename BlockType<A>::type block_type;
+  typedef typename BlockExpression<ExpressionOpExpression>::type BlockExpr;
+  ExpressionOpExpression (const Expression<A> & _a, const Expression<B> & _b)
+    : a(_a), b(_b) {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      assert(a.N() == b.N());
+#endif
+  }
+  BlockExpr operator[] (int i) const {
+    return a[i] OP b[i];
+  }
+  int N() const { return a.N(); };
+private:
+  Expression<A> a;
+  Expression<B> b;
+};
+
+// Expression op Expression
+template <class A, class B>
+Expression< ExpressionOpExpression<A,B> >
+operator OP (const Expression<A> & a, const Expression<B> & b)
+{
+  ExpressionOpExpression<A,B> ex(a, b);
+  return Expression< ExpressionOpExpression<A,B> >(ex);
+}
+
+// Expression op Vector
+template <class A, class B>
+Expression< ExpressionOpExpression<A, ConstRef<B> > >
+operator OP (const Expression<A> & a, const Vector<B> & b)
+{
+  ConstRef<B> rb(b);
+  ExpressionOpExpression<A, ConstRef<B> > ex(a, rb);
+  return Expression< ExpressionOpExpression<A, ConstRef<B> > >(ex);
+}
+
+// Vector op Expression
+template <class A, class B>
+Expression< ExpressionOpExpression< ConstRef<A>, B> >
+operator OP (const Vector<A> & a, const Expression<B> & b)
+{
+  ConstRef<A> ra(a);
+  ExpressionOpExpression< ConstRef<A>, B> ex(ra, b);
+  return Expression< ExpressionOpExpression<ConstRef<A>, B> >(ex);
+}
+
+// Vector op Vector
+template <class V>
+Expression< ExpressionOpExpression< ConstRef<V>, ConstRef<V> > >
+operator OP (const Vector<V> & a, const Vector<V> & b)
+{
+  ConstRef<V> ra(a);
+  ConstRef<V> rb(b);
+  ExpressionOpExpression< ConstRef<V>, ConstRef<V> > ex(ra, rb);
+  return Expression< ExpressionOpExpression< ConstRef<V>, ConstRef<V> > >(ex);
+}
+
+// TypeTraits
+template<class A, class B>
+struct BlockExpression< ExpressionOpExpression<A,B> >
+{
+  typedef
+  typename SelectType<
+    // decide whether we are at the and of the recusrion
+    isEndOfExpressionRecusion<typename BlockType<A>::type>::value,
+    // select FieldType if where are at the end of the recursion
+    typename FieldType<A>::type,
+    // construct nested expression otherwise
+    ExprTmpl::Expression<
+    ExprTmpl::ExpressionOpExpression<
+      typename ExpressionImp< typename BlockExpression<
+        typename ExpressionImp<A>::type>::type>::type,
+      typename ExpressionImp< typename BlockExpression<
+        typename ExpressionImp<B>::type>::type>::type > >
+  >::Type type;
+};
+
+template <class A, class B>
+struct ExpressionImp< ExpressionOpExpression<A,B> >
+{
+  typedef ExpressionOpExpression<A,B> type;
+};
+
+}
+
+template <class A, class B>
+struct FieldType< ExprTmpl::ExpressionOpExpression<A,B> >
+{
+  typedef typename FieldType<A>::type type;
+};
+  
+template <class A, class B>
+struct BlockType< ExprTmpl::ExpressionOpExpression<A,B> >
+{
+  typedef typename BlockType<A>::type type;
+};
+  
+#undef ExpressionOpExpression
+#undef OP

common/exprtmpl/scalar.inc

+// -*- C++ -*-
+#ifndef OP
+#error OP undefined
+#endif
+
+#ifdef ExpressionOpScalar
+
+namespace ExprTmpl {
+
+template <class A>
+class ExpressionOpScalar
+{
+public:
+  typedef typename FieldType<A>::type field_type;
+  typedef typename BlockType<A>::type block_type;
+  typedef typename BlockExpression<ExpressionOpScalar>::type BlockExpr;
+  ExpressionOpScalar (const Expression<A> & _a,
+            const typename FieldType<A>::type & _lambda)
+    : a(_a), lambda(_lambda) {}
+  BlockExpr operator[] (int i) const {
+    return a[i] OP lambda;
+  }
+  int N() const { return a.N(); };
+private:
+  Expression<A> a;
+  const field_type lambda;
+};
+
+// Vector op Scalar
+template <class A>
+Expression< ExpressionOpScalar<ConstRef<A> > >
+operator OP (const Vector<A> & a, const typename FieldType<A>::type & lambda)
+{
+  ConstRef<A> ra(a);
+  ExpressionOpScalar< ConstRef<A> > ex(ra, lambda);
+  return Expression< ExpressionOpScalar<ConstRef<A> > >(ex);
+}
+
+// Expression op Scalar
+template <class A>
+Expression< ExpressionOpScalar<A> >
+operator OP (const Expression<A> & a, const typename FieldType<A>::type & lambda)
+{
+  ExpressionOpScalar<A> ex(a, lambda);
+  return Expression< ExpressionOpScalar<A> >(ex);
+}
+
+// TypeTraits
+//#warning Kann hier nicht das innere ExpressionImp weg?
+template <class Ex>
+struct BlockExpression< ExprTmpl::ExpressionOpScalar<Ex> >
+{
+  typedef
+  typename SelectType<
+    isEndOfExpressionRecusion< typename BlockType<Ex>::type >::value,
+    typename FieldType<Ex>::type,
+    ExprTmpl::Expression<
+    ExprTmpl::ExpressionOpScalar<typename ExpressionImp<
+      typename BlockExpression<
+        typename ExpressionImp<Ex>::type>::type>::type > >
+  >::Type type;
+};
+
+template <class Ex>
+struct ExpressionImp< ExprTmpl::ExpressionOpScalar<Ex> >
+{
+  typedef ExprTmpl::ExpressionOpScalar< Ex > type;
+};
+
+} // namespace ExprTmpl
+
+template <class A>
+struct FieldType< ExprTmpl::ExpressionOpScalar<A> >
+{
+  typedef typename FieldType<A>::type type;
+};
+
+template <class A>
+struct BlockType< ExprTmpl::ExpressionOpScalar<A> >
+{
+  typedef typename BlockType<A>::type type;
+};
+  
+#undef ExpressionOpScalar
+
+#endif
+
+#ifdef ScalarOpExpression
+
+namespace ExprTmpl {
+
+template <class A>
+class ScalarOpExpression
+{
+public:
+  typedef typename FieldType<A>::type field_type;
+  typedef typename BlockType<A>::type block_type;
+  typedef typename BlockExpression<ScalarOpExpression>::type BlockExpr;
+  ScalarOpExpression (const Expression<A> & _a,
+            const typename FieldType<A>::type & _lambda)
+    : a(_a), lambda(_lambda) {}
+  BlockExpr operator[] (int i) const {
+    return lambda OP a[i];
+  }
+  int N() const { return a.N(); };
+private:
+  Expression<A> a;
+  const field_type lambda;
+};
+
+// Scalar op Vector
+template <class A>
+Expression< ScalarOpExpression<ConstRef<A> > >
+operator OP (const typename FieldType<A>::type & lambda, const Vector<A> & a)
+{
+  ConstRef<A> ra(a);
+  ScalarOpExpression< ConstRef<A> > ex(ra, lambda);
+  return Expression< ScalarOpExpression<ConstRef<A> > >(ex);
+}
+
+// Scalar op Expression
+template <class A>
+Expression< ScalarOpExpression<A> >
+operator OP (const typename FieldType<A>::type & lambda, const Expression<A> & a)
+{
+  ScalarOpExpression<A> ex(a, lambda);
+  return Expression< ScalarOpExpression<A> >(ex);
+}
+
+// TypeTraits
+/*
+template<class K, int N>
+struct BlockExpression< ExprTmpl::ScalarOpExpression< ExprTmpl::ConstRef< FieldVector<K,N> > > >
+{
+  typedef K type;
+};
+*/
+template <class Ex>
+struct BlockExpression< ExprTmpl::ScalarOpExpression<Ex> >
+{
+  typedef
+  typename SelectType<
+    isEndOfExpressionRecusion< typename BlockType<Ex>::type >::value,
+    typename FieldType<Ex>::type,
+    ExprTmpl::Expression<
+    ExprTmpl::ScalarOpExpression<typename ExpressionImp<
+      typename BlockExpression<
+        typename ExpressionImp<Ex>::type>::type>::type > >
+  >::Type type;
+};
+
+template <class Ex>
+struct ExpressionImp< ExprTmpl::ScalarOpExpression<Ex> >
+{
+  typedef ExprTmpl::ScalarOpExpression< Ex > type;
+};
+
+} // namespace ExprTmpl
+
+template <class A>
+struct FieldType< ExprTmpl::ScalarOpExpression<A> >
+{
+  typedef typename FieldType<A>::type type;
+};
+
+template <class A>
+struct BlockType< ExprTmpl::ScalarOpExpression<A> >
+{
+  typedef typename BlockType<A>::type type;
+};
+  
+#undef ScalarOpExpression
+
+#endif
+
+#undef OP

common/fvector.hh

@@ -11,6 +11,10 @@
 #include "exceptions.hh"
 #include "genericiterator.hh"
 
+#ifdef DUNE_EXPRESSIONTEMPLATES
+#include "exprtmpl.hh"
+#endif
+
 namespace Dune {
 
   /** @defgroup DenseMatVec Dense Matrix and Vector Template Library
@@ -353,11 +357,16 @@ namespace Dune {
       : container_(other.container_), position_(other.position_)
     {}
 
-
+#if 0
     FieldIterator(const FieldIterator<const typename Dune::RemoveConst<C>::Type, const typename Dune::RemoveConst<T>::Type >& other)
       : container_(other.container_), position_(other.position_)
     {}
-
+#endif
+#if 0
+    FieldIterator(const FieldIterator<C,T>& other)
+      : container_(other.container_), position_(other.position_)
+    {}
+#endif
     // Methods needed by the forward iterator
     bool equals(const FieldIterator<typename Dune::RemoveConst<C>::Type,typename Dune::RemoveConst<T>::Type>& other) const
     {
@@ -415,6 +424,135 @@ namespace Dune {
     DifferenceType position_;
   };
 
+  //! Type Traits for Vector::Iterator vs (const Vector)::ConstIterator
+  template<class T>
+  struct IteratorType
+  {
+    typedef typename T::Iterator type;
+  };
+
+  template<class T>
+  struct IteratorType<const T>
+  {
+    typedef typename T::ConstIterator type;
+  };
+
+#ifdef DUNE_EXPRESSIONTEMPLATES
+  //! Iterator class for flat sequential access to a nested Vector
+  template<class V>
+  class FlatIterator :
+    public ForwardIteratorFacade<FlatIterator<V>,
+        typename Dune::FieldType<V>::type,
+        typename Dune::FieldType<V>::type&,
+        int>
+  {
+  public:
+    typedef typename IteratorType<V>::type BlockIterator;
+    typedef std::ptrdiff_t DifferenceType;
+    //    typedef typename BlockIterator::DifferenceType DifferenceType;
+    typedef typename BlockType<V>::type block_type;
+    typedef typename FieldType<V>::type field_type;
+    typedef FlatIterator<block_type> SubBlockIterator;
+    FlatIterator(const BlockIterator & i) :
+      it(i), bit(i->begin()), bend(i->end()) {};
+    void increment ()
+    {
+      ++bit;
+      if (bit == bend)
+      {
+        ++it;
+        bit = it->begin();
+        bend = it->end();
+      }
+    }
+    bool equals (const FlatIterator & fit) const
+    {
+      return fit.it == it;
+    }
+    field_type& dereference() const
+    {
+      return *bit;
+    }
+    //! return index
+    DifferenceType index ()
+    {
+      return bit.index();
+    }
+  private:
+    BlockIterator it;
+    SubBlockIterator bit;
+    SubBlockIterator bend;
+  };
+
+  //! Specialization for FieldVector
+  //! acts as the end of the recursive template
+  template<class K, int N>
+  class FlatIterator< FieldVector<K,N> > :
+    public ForwardIteratorFacade<FlatIterator< FieldVector<K,N> >,
+        K, K&, int>
+  {
+  public:
+    typedef typename FieldVector<K,N>::Iterator BlockIterator;
+    typedef std::ptrdiff_t DifferenceType;
+    //    typedef typename BlockIterator::DifferenceType DifferenceType;
+    typedef typename FieldVector<K,N>::field_type field_type;
+    FlatIterator(const BlockIterator & i) : it(i) {};
+    void increment ()
+    {
+      ++it;
+    }
+    bool equals (const FlatIterator & fit) const
+    {
+      return fit.it == it;
+    }
+    field_type& dereference() const
+    {
+      return *it;
+    }
+    //! return index
+    DifferenceType index ()
+    {
+      return it.index();
+    }
+  private:
+    BlockIterator it;
+  };
+
+  //! Specialization for const FieldVector
+  //! acts as the end of the recursive template
+  template<class K, int N>
+  class FlatIterator< const FieldVector<K,N> > :
+    public ForwardIteratorFacade<FlatIterator< const FieldVector<K,N> >,
+        const K, const K&, int>
+  {
+  public:
+    typedef typename FieldVector<K,N>::ConstIterator BlockIterator;
+    typedef std::ptrdiff_t DifferenceType;
+    //    typedef typename BlockIterator::DifferenceType DifferenceType;
+    typedef typename FieldVector<K,N>::field_type field_type;
+    FlatIterator(const BlockIterator & i) : it(i) {};
+    void increment ()
+    {
+      ++it;
+    }
+    bool equals (const FlatIterator & fit) const
+    {
+      return fit.it == it;
+    }
+    const field_type& dereference() const
+    {
+      return *it;
+    }
+    //! return index
+    DifferenceType index ()
+    {
+      return it.index();
+    }
+  private:
+    BlockIterator it;
+  };
+#endif
+
   /** \brief Construct a vector space out of a tensor product of fields.
 
          K is the field type (use float, double, complex, etc) and n
@@ -427,6 +565,9 @@ namespace Dune {
    */
   template<class K, int SIZE>
   class FieldVector
+#ifdef DUNE_EXPRESSIONTEMPLATES
+    : public Dune::ExprTmpl::Vector< FieldVector<K,SIZE> >
+#endif
   {
     //! The actual number of elements that gets allocated.
     //! It's always at least 1.
@@ -479,6 +620,28 @@ namespace Dune {
       return *this;
     }
 
+#ifdef DUNE_EXPRESSIONTEMPLATES
+    template <class E>
+    FieldVector (Dune::ExprTmpl::Expression<E> op) {
+      Dune::dvverb << INDENT << "FieldVector Copy Constructor Expression\n";
+      assignFrom(op);
+    }
+    template <class V>
+    FieldVector (const Dune::ExprTmpl::Vector<V> & op) {
+      Dune::dvverb << INDENT << "FieldVector Copy Operator Vector\n";
+      assignFrom(op);
+    }
+    template <class E>
+    FieldVector&  operator = (Dune::ExprTmpl::Expression<E> op) {
+      Dune::dvverb << INDENT << "FieldVector Assignment Operator Expression\n";
+      return assignFrom(op);
+    }
+    template <class V>
+    FieldVector& operator = (const Dune::ExprTmpl::Vector<V> & op) {
+      Dune::dvverb << INDENT << "FieldVector Assignment Operator Vector\n";
+      return assignFrom(op);
+    }
+#endif
 
     //===== access to components
 
@@ -592,6 +755,7 @@ namespace Dune {
       return *this;
     }
 
+#ifndef DUNE_EXPRESSIONTEMPLATES
     //! Binary vector addition
     FieldVector<K, size> operator+ (const FieldVector<K, size>& b) const
     {
@@ -605,6 +769,7 @@ namespace Dune {
       FieldVector<K, size> z = *this;
       return (z-=b);
     }
+#endif
 
     //! vector space add scalar to all comps
     FieldVector& operator+= (const K& k)
@@ -741,6 +906,9 @@ namespace Dune {
    */
   template<class K>
   class FieldVector<K,1>
+#ifdef DUNE_EXPRESSIONTEMPLATES
+    : public Dune::ExprTmpl::Vector< FieldVector<K,1> >
+#endif
   {
     enum { n = 1 };
   public:
@@ -785,6 +953,29 @@ namespace Dune {
       return *this;
     }
 
+#ifdef DUNE_EXPRESSIONTEMPLATES
+    template <class E>
+    FieldVector (Dune::ExprTmpl::Expression<E> op) {
+      Dune::dvverb << INDENT << "FieldVector<1> Copy Constructor Expression\n";
+      assignFrom(op);
+    }
+    template <class V>
+    FieldVector (const Dune::ExprTmpl::Vector<V> & op) {
+      Dune::dvverb << INDENT << "FieldVector<1> Copy Operator Vector\n";
+      assignFrom(op);
+    }
+    template <class E>
+    FieldVector&  operator = (Dune::ExprTmpl::Expression<E> op) {
+      Dune::dvverb << INDENT << "FieldVector<1> Assignment Operator Expression\n";
+      return assignFrom(op);
+    }
+    template <class V>
+    FieldVector& operator = (const Dune::ExprTmpl::Vector<V> & op) {
+      Dune::dvverb << INDENT << "FieldVector<1> Assignment Operator Vector\n";
+      return assignFrom(op);
+    }
+#endif
+
     //===== access to components
 
     //! random access
@@ -896,6 +1087,7 @@ namespace Dune {
       return *this;
     }
 
+#ifndef DUNE_EXPRESSIONTEMPLATES
     //! Binary vector addition
     FieldVector operator+ (const FieldVector& b) const
     {
@@ -909,6 +1101,8 @@ namespace Dune {
       FieldVector z = *this;
       return (z-=b);
     }
+#endif
+
     //! vector space add scalar to each comp
     FieldVector& operator+= (const K& k)
     {

common/test/.gitignore

@@ -17,4 +17,5 @@ poolallocatortest
 *.gcno
 gmon.out
 gcdlcdtest
-streamtest
+streamtest
\ No newline at end of file
+exprtmpl
\ No newline at end of file

common/test/Makefile.am

@@ -2,7 +2,7 @@
 
 TESTPROGS = parsetest test-stack arraylisttest smartpointertest \
 	sllisttest iteratorfacadetest tuplestest fmatrixtest \
-	poolallocatortest settest gcdlcdtest streamtest
+	poolallocatortest settest gcdlcdtest streamtest exprtmpl
 
 # which tests to run
 TESTS = $(TESTPROGS)
@@ -42,4 +42,9 @@ settest_SOURCES=settest.cc
 
 gcdlcdtest_SOURCES = gcdlcdtest.cc
 
+#bin_PROGRAMS = exprtmpl
+exprtmpl_SOURCES = exprtmpl.cc
+exprtmpl_CXXFLAGS = -DDUNE_ISTL_WITH_CHECKING -DDUNE_EXPRESSIONTEMPLATES
+exprtmpl_LDFLAGS = $(top_builddir)/common/libcommon.la
+
 include $(top_srcdir)/am/global-rules

common/test/exprtmpl.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <iostream>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/io.hh>
+#include <dune/common/iteratorfacades.hh>
+
+Indent INDENT;
+
+void test_fvector()
+{
+  typedef Dune::FieldVector<double,2> VB;
+  VB v1(1);
+  VB v2(2);
+
+  typedef Dune::ExprTmpl::ConstRef<VB> RVB;
+  VB v = 0.5 * (v1 + v2 * 2) + 3 * v1 - v2;
+  std::cout << " 0.5 * ( " << v1 << " + " << v2 << " * 2) + 3 * " << v1 << " - " << v2 << std::endl;
+}
+
+void test_blockvector()
+{
+  Dune::FieldVector<double,2> v(10);
+  typedef Dune::FieldVector<double,2> VB;
+  typedef Dune::BlockVector<VB> BV;
+
+  const int sz = 3;
+  BV bv1(sz), bv2(sz);
+  bv1 = 1;
+  bv2 = 0;
+  bv2[1][0]=1;
+  bv2[1][1]=2;
+
+  BV bv(sz);
+  bv = -17;
+  printvector (std::cout, bv, "bv", "r");
+  //  bv.emptyClone(bv1);
+  std::cout << "Assingn from ConstRef\n";
+  bv = 2 * (bv1 + bv2);
+  bv -= 1;
+
+  printvector (std::cout, bv1, "bv1", "r");
+  printvector (std::cout, bv2, "bv2", "r");
+  printvector (std::cout, bv, "bv", "r");
+
+}
+
+void test_blockblockvector()
+{
+  typedef Dune::FieldVector<double,2> VB;
+  typedef Dune::BlockVector<VB> BV;
+  typedef Dune::BlockVector<BV> BBV;
+  typedef Dune::ExprTmpl::ConstRef<BV> RBV;
+  const int sz = 3;
+  BV bv1(sz), bv2(sz);
+  bv1 = 1;
+  bv2 = 0;
+  bv2[1][0]=1;
+  bv2[1][1]=2;
+
+  Dune::ExprTmpl::ConstRef<BV> rbv1(bv1);
+  Dune::ExprTmpl::ConstRef<BV> rbv2(bv2);
+
+  BBV bbv(2);
+  bbv[0].resize(bv1.N());
+  bbv[0] = Dune::ExprTmpl::Expression<RBV>(rbv1);
+  bbv[1].resize(bv2.N());
+  bbv[1] = Dune::ExprTmpl::Expression<RBV>(rbv2);
+  Dune::FlatIterator<BBV> fit(bbv.begin());
+  Dune::FlatIterator<BBV> fend(bbv.end());
+  int index = 0;
+  for(; fit!=fend; ++fit)
+  {
+    BBV::field_type x;
+    x = *fit;
+    std::cout << index << "\t" << x << std::endl;
+    index++;
+  }
+  printvector (std::cout, bv1, "bv1", "r");
+  printvector (std::cout, bv2, "bv1", "r");
+  printvector (std::cout, bbv, "bbv", "r");
+}
+
+namespace Dune {
+
+  /**
+     B: rhs Vector type (usually a BlockType)
+     M: lhs Matrix type
+     V: lhs Vector type
+   */
+  // !TODO!
+  template <class B, class M, class V>
+  class BlockRowSum
+  {
+  public:
+    typedef typename M::ConstColIterator ConstColIterator;
+    typedef typename M::row_type row_type;
+    typedef BlockRowSum<typename BlockType<B>::type,M,V> SubBlockRowSum;
+    typedef typename Dune::FieldType<V>::type field_type;
+    //  typedef typename
+    BlockRowSum(const row_type & _r, const V & _v) : r(_r), v(_v) {};
+    SubBlockRowSum operator[] (int i) const {
+      return SubBlockRowSum(r[i],v);
+    }
+    int N() const { return r.begin()->N(); }
+  private:
+    const row_type & r;
+    const V & v;
+  };
+
+  // template <class K, int n, class M, class V>
+  // class BlockRowSum< FieldVector<K,n>, M, V >
+  template <class M, class V>
+  class BlockRowSum< void, M, V >
+  {
+  public:
+    typedef typename M::ConstColIterator ConstColIterator;
+    typedef typename M::row_type row_type;
+    typedef typename FieldType<V>::type K;
+    BlockRowSum(const row_type & _r, const V & _v) : r(_r), v(_v) {};
+    K operator[] (int i) const {
+      K x=0;
+      ConstColIterator j=r.begin();
+      ConstColIterator endj = r.end();
+      for (; j!=endj; ++j)
+      {
+        /*
+           x += (MatrixBlock * VectorBlock  )[i]
+         */
+        //      std::cout << "BRS x=" << x << std::endl;
+        x   += (     (*j)   * v[j.index()] )[i];
+      }
+      //    std::cout << "BRS x=" << x << std::endl;
+      return x;
+    }
+    int N() const { return r.begin()->N(); }
+  private:
+    const row_type & r;
+    const V & v;
+  };
+
+  template <class M, class V>
+  class MV
+  {
+  public:
+    MV(const M & _A, const V & _v) : A(_A), v(_v) {};
+    BlockRowSum<void,M,V> operator[] (int i) const {
+      // !TODO!
+      return BlockRowSum<void,M,V>(A[i],v);
+    }
+    int N() const { return A.N(); }
+  private:
+    const M & A;
+    const V & v;
+  };
+
+  template <class K, int iM, int iN>
+  class MV< Dune::FieldMatrix<K,iN,iM>, Dune::FieldVector<K,iM> >
+  {
+  public:
+    typedef Dune::FieldMatrix<K,iN,iM> M;
+    typedef Dune::FieldVector<K,iM> V;
+    typedef typename Dune::FieldType<V>::type field_type;
+    MV(const M & _A, const V & _v) : A(_A), v(_v) {};
+    field_type operator[] (int i) const {
+      field_type x=0;
+      for (int j=0; j<iM; ++j) {
+        x += A[i][j] * v[j];
+      }
+      return x;
+    }
+    int N() const { return A.N(); }
+  private:
+    const M & A;
+    const V & v;
+  };
+
+  template<class K, int N, int M>
+  ExprTmpl::Expression< MV< FieldMatrix<K,N,M>, FieldVector<K,M> > >
+  operator * ( const FieldMatrix<K,N,M> & A, const FieldVector<K,M> & v )
+  {
+    return ExprTmpl::Expression< MV< FieldMatrix<K,N,M>, FieldVector<K,M> > >
+             ( MV< FieldMatrix<K,N,M>, FieldVector<K,M> > (A,v) );
+  }
+
+  template<class BM, class BV>
+  ExprTmpl::Expression< MV< BCRSMatrix<BM>, BlockVector<BV> > >
+  operator * ( const BCRSMatrix<BM> & A, const BlockVector<BV> & v )
+  {
+    return ExprTmpl::Expression< MV< BCRSMatrix<BM>, BlockVector<BV> > >
+             ( MV< BCRSMatrix<BM>, BlockVector<BV> > (A,v) );
+  }
+
+  template <class A, class B>
+  struct FieldType< MV<A,B> >
+  {
+    typedef typename FieldType<B>::type type;
+  };
+
+  template <class B, class A, class V>
+  struct FieldType< BlockRowSum<B,A,V> >
+  {
+    typedef typename FieldType<V>::type type;
+  };
+
+  namespace ExprTmpl {
+
+    // !TODO!
+    template <class A, class B>
+    struct BlockExpression< MV< BCRSMatrix<A>,BlockVector<B> > >
+    {
+      typedef ExprTmpl::Expression< BlockRowSum< void, BCRSMatrix<A>,BlockVector<B> > > type;
+    };
+
+    // !TODO!
+    template <class A, class B>
+    struct BlockExpression< BlockRowSum< void, BCRSMatrix<A>,BlockVector<B> > >
+    {
+      typedef typename B::field_type type;
+    };
+
+    template <class K, int N, int M>
+    struct BlockExpression< MV< FieldMatrix<K,N,M>, FieldVector<K,M> > >
+    {
+      typedef K type;
+    };
+
+  }
+}
+
+void test_matrix()
+{
+  static const int BlockSize = 2;
+  typedef double matvec_t;
+  typedef Dune::FieldVector<matvec_t,BlockSize+1> VB;
+  typedef Dune::FieldVector<matvec_t,BlockSize> LVB;
+  typedef Dune::FieldMatrix<matvec_t,BlockSize,BlockSize+1> MB;
+  typedef Dune::BlockVector<LVB> LeftVector;
+  typedef Dune::BlockVector<VB> Vector;
+  typedef Dune::BCRSMatrix<MB> Matrix;
+
+  LVB a(0);
+  VB b(2);
+  MB _M(1);
+  _M[1][1] = 3;
+
+  // a += M * b
+  _M.umv(b,a);
+
+  printmatrix (std::cout, _M, "Matrix", "r");
+  printvector (std::cout, a, "Vector", "r");
+
+  // a = M * b
+  a = _M*b;
+
+  printvector (std::cout, a, "Vector", "r");
+
+  int N = 4;
+  int M = 5;
+
+  Matrix A(N,M,Matrix::row_wise);
+  Matrix::CreateIterator i=A.createbegin();
+  Matrix::CreateIterator end=A.createend();
+  // build up the matrix structure
+  int c=0;
+  for (; i!=end; ++i)
+  {
+    // insert a non zero entry for myself
+    i.insert(c);
+    // insert index M-1
+    i.insert(M-1);
+    c++;
+  }
+  A = 0.0;
+
+  std::cout << "Matrix coldim=" << A.coldim() << std::endl;
+  std::cout << "Matrix rowdim=" << A.rowdim() << std::endl;
+  std::cout << "Matrix N=" << A.M() << std::endl;
+  std::cout << "Matrix M=" << A.N() << std::endl;
+
+  Matrix::Iterator cit=A.begin();
+  Matrix::Iterator cend=A.end();
+  for (; cit!=cend; ++cit)
+  {
+    *cit = cit.index();
+  }
+
+  printmatrix (std::cout, A, "Matrix", "r");
+
+  LeftVector v(N);
+  LeftVector v2(N);
+  v = 0;
+  Vector x(M);
+  x = 1;
+  x[M-1] = 2;
+
+  std::cout << A.M() << " " << x.N() << " " << v.N() << std::endl;
+
+  A.umv(x,v);
+
+  printvector (std::cout, v, "Vector", "r");
+  /*
+     v=1;
+     MV<Matrix,Vector> eximp (A,x);
+     Dune::ExprTmpl::Expression< MV<Matrix,Vector> > ex ( eximp );
+     v =
+     printvector (std::cout, v, "Vector", "r");
+   */
+  v2 = A * x;
+
+  printvector (std::cout, v2, "Vector", "r");
+
+}
+
+int main()
+{
+  //  Dune::dvverb.attach(std::cout);
+  try
+  {
+    test_fvector();
+    test_blockvector();
+    test_blockblockvector();
+    test_matrix();
+    exit(0);
+  }
+  catch (Dune::Exception & e)
+  {
+    std::cout << e << std::endl;
+  }
+}

istl/bvector.hh

@@ -10,6 +10,10 @@
 #include "allocator.hh"
 #include "basearray.hh"
 
+#ifdef DUNE_EXPRESSIONTEMPLATES
+#include <dune/common/exprtmpl.hh>
+#endif
+
 /*! \file
 
    \brief  This file implements a vector space as a tensor product of
@@ -63,8 +67,8 @@ namespace Dune {
 
 
     //===== assignment from scalar
-
     //! Assignment from a scalar
+
     block_vector_unmanaged& operator= (const field_type& k)
     {
       for (size_type i=0; i<this->n; i++)
@@ -72,9 +76,8 @@ namespace Dune {
       return *this;
     }
 
-
     //===== vector space arithmetic
-
+#ifndef DUNE_EXPRESSIONTEMPLATES
     //! vector space addition
     block_vector_unmanaged& operator+= (const block_vector_unmanaged& y)
     {
@@ -108,7 +111,7 @@ namespace Dune {
       for (int i=0; i<this->n; ++i) (*this)[i] /= k;
       return *this;
     }
-
+#endif
     //! vector space axpy operation
     block_vector_unmanaged& axpy (const field_type& a, const block_vector_unmanaged& y)
     {
@@ -221,6 +224,9 @@ namespace Dune {
    */
   template<class B, class A=ISTLAllocator>
   class BlockVector : public block_vector_unmanaged<B,A>
+#ifdef DUNE_EXPRESSIONTEMPLATES
+                      , public Dune::ExprTmpl::Vector< Dune::BlockVector<B,A> >
+#endif
   {
   public:
 
@@ -269,6 +275,56 @@ namespace Dune {
       }
     }
 
+#ifdef DUNE_EXPRESSIONTEMPLATES
+    //! random access to blocks
+    B& operator[] (size_type i)
+    {
+      return block_vector_unmanaged<B,A>::operator[](i);
+    }
+
+    //! same for read only access
+    const B& operator[] (size_type i) const
+    {
+      return block_vector_unmanaged<B,A>::operator[](i);
+    }
+
+    //! dimension of the vector space
+    int N () const
+    {
+      return block_vector_unmanaged<B,A>::N();
+    }
+
+    BlockVector (const BlockVector& a) {
+      Dune::dvverb << INDENT << "BlockVector Copy Constructor BlockVector\n";
+      assignFrom(a);
+    }
+    template <class V>
+    BlockVector (Dune::ExprTmpl::Expression<V> op) {
+      Dune::dvverb << INDENT << "BlockVector Copy Constructor Expression\n";
+      assignFrom(op);
+    }
+    template <class V>
+    BlockVector (const Dune::ExprTmpl::Vector<V> & op) {
+      Dune::dvverb << INDENT << "BlockVector Copy Constructor Vector\n";
+      assignFrom(op);
+    }
+    BlockVector& operator = (const BlockVector& a) {
+      Dune::dvverb << INDENT << "BlockVector Assignment Operator BlockVector\n";
+      return assignFrom(a);
+    }
+    template <class E>
+    BlockVector&  operator = (Dune::ExprTmpl::Expression<E> op) {
+      Dune::dvverb << INDENT << "BlockVector Assignment Operator Expression\n";
+      return assignFrom(op);
+    }
+    template <class V>
+    BlockVector& operator = (const Dune::ExprTmpl::Vector<V> & op) {
+      Dune::dvverb << INDENT << "BlockVector Assignment Operator Vector\n";
+      return assignFrom(op);
+    }
+#endif
+
+#ifndef DUNE_EXPRESSIONTEMPLATES
     //! copy constructor
     BlockVector (const BlockVector& a) :
       block_vector_unmanaged<B,A>(a)
@@ -306,7 +362,7 @@ namespace Dune {
       // and copy elements
       for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
     }
-
+#endif
 
     //! free dynamic memory
     ~BlockVector ()
@@ -331,6 +387,7 @@ namespace Dune {
     }
 
     //! assignment
+#ifndef DUNE_EXPRESSIONTEMPLATES
     BlockVector& operator= (const BlockVector& a)
     {
       if (&a!=this)     // check if this and a are different objects
@@ -360,6 +417,7 @@ namespace Dune {
       // forward to regular assignement operator
       return this->operator=(static_cast<const BlockVector&>(a));
     }
+#endif
 
     //! assign from scalar
     BlockVector& operator= (const field_type& k)
@@ -883,7 +941,18 @@ namespace Dune {
     }
   };
 
-
+#ifdef DUNE_EXPRESSIONTEMPLATES
+  template <class B, class A>
+  struct FieldType< BlockVector<B,A> >
+  {
+    typedef typename FieldType<B>::type type;
+  };
+  template <class B, class A>
+  struct BlockType< BlockVector<B,A> >
+  {
+    typedef B type;
+  };
+#endif
 
   /** @} end documentation */