From 93dd050187b4add073c4ac0b16dacfeffdfea67e Mon Sep 17 00:00:00 2001
From: Christian Engwer <christi@dune-project.org>
Date: Tue, 2 Aug 2005 12:24:59 +0000
Subject: [PATCH] latest expression template updates
[[Imported from SVN: r2502]]
---
common/exprtmpl.hh | 505 +++++++++++++++++++++++++++++++++++++--
common/fmatrix.hh | 22 ++
common/test/.gitignore | 4 +-
common/test/Makefile.am | 13 +-
common/test/exprtmpl.cc | 512 +++++++++-------------------------------
common/test/timing.cc | 107 +++++++++
istl/bcrsmatrix.hh | 24 ++
istl/bvector.hh | 7 +-
8 files changed, 771 insertions(+), 423 deletions(-)
create mode 100644 common/test/timing.cc
diff --git a/common/exprtmpl.hh b/common/exprtmpl.hh
index 056ced1c8..a4964dc92 100644
--- a/common/exprtmpl.hh
+++ b/common/exprtmpl.hh
@@ -3,6 +3,22 @@
#ifndef DUNE_EXPRTMPL_HH
#define DUNE_EXPRTMPL_HH
+/*
+ This file is part of DUNE, a Distributed and Unified Numerics Environment
+ It is distributed under the terms of the GNU Lesser General Public License version 2.1
+ See COPYING at the top of the source tree for the full licence.
+ */
+
+/*! @file
+
+ @brief This file provides expression templates for the
+ «Dense Matrix and Vector Template Library» and for the
+ «Iterative Solvers Template Library».
+ @verbatim
+ $Id$
+ @endverbatim
+ */
+
#include <iostream>
#include <iomanip>
#include <cstdlib>
@@ -37,6 +53,11 @@ namespace Dune {
template<class V> class FlatIterator;
template<class K, int N> class FieldVector;
+ template<class K, int N, int M> class FieldMatrix;
+#warning this header should not know about BCRSMatrix and BlockVector
+ class ISTLAllocator;
+ template<class B, class A=ISTLAllocator> class BCRSMatrix;
+ template<class B, class A=ISTLAllocator> class BlockVector;
/**
Type Traits for field_type and block_type
@@ -63,13 +84,13 @@ namespace Dune {
{
typedef int type;
};
- // FieldType specialization for FieldMatrix
+ // FieldType specialization for FieldVector
template <class K, int N>
struct FieldType< FieldVector<K,N> >
{
typedef K type;
};
- // BlockType specialization for FieldMatrix
+ // BlockType specialization for FieldVector
template <class K, int N>
struct BlockType< FieldVector<K,N> >
{
@@ -77,15 +98,15 @@ namespace Dune {
};
// FieldType specialization for const T
template<class T>
- struct BlockType<const T>
+ struct FieldType<const T>
{
- typedef const typename BlockType<T>::type type;
+ typedef const typename FieldType<T>::type type;
};
// BlockType specialization for const T
template<class T>
- struct FieldType<const T>
+ struct BlockType<const T>
{
- typedef const typename FieldType<T>::type type;
+ typedef const typename BlockType<T>::type type;
};
namespace ExprTmpl {
@@ -95,6 +116,7 @@ namespace Dune {
template <class Ex> class Expression;
template <class A, class B, template<class> class Op> class ExBinOp;
template <class I> class Vector;
+ template <class I> class Matrix;
/**
Type Trait for nested Expression
@@ -196,8 +218,7 @@ namespace Dune {
};
/**
- 1 Dimensional Vector Base Class
- for Glommable Expr1 Templates
+ Vector Base Class for Expression Templates
*/
template <class I>
class Vector {
@@ -266,25 +287,22 @@ namespace Dune {
--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);
+ 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);
+ 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);
+ 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);
+ 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();
@@ -363,12 +381,463 @@ namespace Dune {
typedef typename BlockType<I>::type type;
};
- // ---
+ /**
+ Type Traits for row_type of Matrix
+ */
+ template<class M>
+ struct RowType
+ {
+ typedef typename M::row_type type;
+ };
+ template<class K, int N, int M>
+ struct RowType< FieldMatrix<K,N,M> >
+ {
+ typedef FieldVector<K,M> type;
+ };
+ template <class I>
+ struct RowType< ExprTmpl::Matrix<I> >
+ {
+ typedef typename RowType<I>::type type;
+ };
+ // RowType specialization for const T
+ template<class T>
+ struct RowType<const T>
+ {
+ typedef const typename RowType<T>::type type;
+ };
+
+ // Matrix-Vector Multiplication
namespace ExprTmpl {
- // ---
+ /**
+ Matrix Base Class for Expression Templates
+ */
+ template <class I>
+ class Matrix {
+ public:
+ explicit Matrix() {}
+ typedef typename RowType<I>::type row_type;
+ typedef typename FieldType<I>::type field_type;
+ //! dimension of the vector space
+ int N() const {
+ return asImp().N();
+ }
+ int M() const {
+ return asImp().M();
+ }
+ row_type & operator[] (int i) {
+ return asImp()[i];
+ }
+ const row_type & operator[] (int i) const {
+ return asImp()[i];
+ }
+ private:
+ I & asImp() { return static_cast<I&>(*this); }
+ const I & asImp() const { return static_cast<const I&>(*this); }
+ };
+
+ // Trait Structs to extract infos needed for Matrix-Vector Multiplication
+ template<class M>
+ struct NestedDepth
+ {
+ enum { value = NestedDepth<typename BlockType<M>::type>::value + 1 };
+ };
- }
+ template<class K, int N, int M>
+ struct NestedDepth< FieldMatrix<K,N,M> >
+ {
+ enum { value = 1 };
+ };
+
+ template<class Me, class M>
+ struct MyDepth
+ {
+ enum { value = MyDepth<Me,typename BlockType<M>::type>::value+1 };
+ };
+
+ template<class Me>
+ struct MyDepth<Me,Me>
+ {
+ enum { value = 0 };
+ };
+
+ template<class B, int i>
+ struct BlockTypeN
+ {
+ typedef typename BlockTypeN<typename BlockType<B>::type, i-1>::type type;
+ };
+
+ template<class B>
+ struct BlockTypeN<B,0>
+ {
+ typedef B type;
+ };
+
+ template<class B>
+ struct BlockTypeN<B,-1>
+ {
+ typedef B type;
+ };
+
+ //! Type Traits for Vector::Iterator vs (const Vector)::ConstIterator
+ template<class T>
+ struct ColIteratorType
+ {
+ typedef typename T::ColIterator type;
+ };
+
+ template<class T>
+ struct ColIteratorType<const T>
+ {
+ typedef typename T::ConstColIterator type;
+ };
+
+ //! Iterator class for flat sequential access to a nested Matrix Row
+ template<class A>
+ class FlatColIterator :
+ public ForwardIteratorFacade<FlatColIterator<A>,
+ typename FieldType<A>::type,
+ typename FieldType<A>::type&,
+ int>
+ {
+ public:
+ typedef typename ColIteratorType<A>::type ColBlockIterator;
+ typedef std::ptrdiff_t DifferenceType;
+ // typedef typename BlockIterator::DifferenceType DifferenceType;
+ typedef typename BlockType<A>::type block_type;
+ typedef typename FieldType<A>::type field_type;
+ typedef FlatColIterator<block_type> SubBlockIterator;
+ FlatColIterator(const ColBlockIterator & i, const int* _M) :
+ M(_M), it(i),
+ bit((*i)[(*M)].begin(), M+1),
+ bend((*i)[(*M)].end(), M+1) {};
+ void increment ()
+ {
+ ++bit;
+ if (bit == bend)
+ {
+ ++it;
+ bit = (*it)[(*M)].begin();
+ bend = (*it)[(*M)].end();
+ }
+ }
+ bool equals (const FlatColIterator & fit) const
+ {
+ return fit.it == it && fit.bit == bit;
+ }
+ const field_type& dereference() const
+ {
+ return *bit;
+ }
+ template<class V>
+ const field_type& vectorentry(const Dune::ExprTmpl::Vector<V> & v) const
+ {
+ return bit.vectorentry(v[it.index()]);
+ }
+ //! return index
+ DifferenceType index () const
+ {
+ return bit.index();
+ }
+ FlatColIterator operator = (const ColBlockIterator & _i)
+ {
+ it = _i;
+ bit = (*it)[(*M)].begin();
+ bend = (*it)[(*M)].end();
+ return *this;
+ }
+ private:
+ const int* M;
+ ColBlockIterator it;
+ SubBlockIterator bit;
+ SubBlockIterator bend;
+ };
+
+ template<class K, int N, int M>
+ class FlatColIterator<FieldMatrix<K,N,M> > :
+ public ForwardIteratorFacade<
+ FlatColIterator< FieldMatrix<K,N,M> >, K, K&, int>
+ {
+ public:
+ typedef
+ typename ColIteratorType< FieldMatrix<K,N,M> >::type ColBlockIterator;
+ typedef std::ptrdiff_t DifferenceType;
+ typedef K field_type;
+ FlatColIterator(const ColBlockIterator & i, const int*) :
+ it(i) {};
+ void increment ()
+ {
+ ++it;
+ }
+ bool equals (const FlatColIterator & fit) const
+ {
+ return fit.it == it;
+ }
+ field_type& dereference() const
+ {
+ return *it;
+ }
+ const field_type& vectorentry(const FieldVector<K,M> & v) const
+ {
+ return v[it.index()];
+ }
+ //! return index
+ DifferenceType index () const
+ {
+ return it.index();
+ }
+ FlatColIterator operator = (const ColBlockIterator & _i)
+ {
+ it = _i;
+ return *this;
+ }
+ private:
+ ColBlockIterator it;
+ };
+
+ template<class K, int N, int M>
+ class FlatColIterator<const FieldMatrix<K,N,M> > :
+ public ForwardIteratorFacade<
+ FlatColIterator< const FieldMatrix<K,N,M> >, const K, const K&, int>
+ {
+ public:
+ typedef
+ typename ColIteratorType< const FieldMatrix<K,N,M> >::type ColBlockIterator;
+ typedef std::ptrdiff_t DifferenceType;
+ typedef const K field_type;
+ FlatColIterator(const ColBlockIterator & i, const int*) :
+ it(i) {};
+ void increment ()
+ {
+ ++it;
+ }
+ bool equals (const FlatColIterator & fit) const
+ {
+ return fit.it == it;
+ }
+ field_type& dereference() const
+ {
+ return *it;
+ }
+ const field_type& vectorentry(const FieldVector<K,M> & v) const
+ {
+ return v[it.index()];
+ }
+ //! return index
+ DifferenceType index () const
+ {
+ return it.index();
+ }
+ FlatColIterator operator = (const ColBlockIterator & _i)
+ {
+ it = _i;
+ return *this;
+ }
+ private:
+ ColBlockIterator it;
+ };
+
+ /**
+ B: BlockMatrix type -> indicated the current level.
+ M: ,,global'' Matrix type
+ V: ,,global'' Vector type
+ */
+ template <class B, class Mat, class Vec>
+ class MatrixMulVector
+ {
+ public:
+ typedef typename
+ BlockTypeN<MatrixMulVector<Mat,Mat,Vec>, MyDepth<B,Mat>::value-1>::type
+ ParentBlockType;
+ typedef
+ MatrixMulVector<typename BlockType<B>::type,Mat,Vec> SubMatrixMulVector;
+ typedef typename
+ Dune::ExprTmpl::BlockExpression< MatrixMulVector<B,Mat,Vec> >::type
+ BlockExpr;
+ typedef
+ typename BlockType<Mat>::type::ColIterator SubColIterator;
+ typedef typename Dune::FieldType<Vec>::type field_type;
+ /* constructor */
+ MatrixMulVector(const Mat & _A, const Vec & _v, int* _M,
+ const ParentBlockType & _parent) :
+ parent(_parent), M(_M), A(_A), v(_v) {};
+ BlockExpr operator[] (int i) const {
+ M[MyDepth<B,Mat>::value] = i;
+ return SubMatrixMulVector(A,v,M,*this);
+ }
+ int N() const { return -1; }; //r.begin()->N(); }
+ const ParentBlockType & parent;
+ private:
+ mutable int* M;
+ const Mat & A;
+ const Vec & v;
+ };
+
+ template <class Mat, class Vec>
+ class MatrixMulVector<Mat,Mat,Vec>
+ {
+ public:
+ typedef
+ MatrixMulVector<typename BlockType<Mat>::type,Mat,Vec> SubMatrixMulVector;
+ typedef
+ typename BlockType<Mat>::type::ColIterator SubColIterator;
+ typedef typename
+ Dune::ExprTmpl::BlockExpression< MatrixMulVector<Mat,Mat,Vec> >::type
+ BlockExpr;
+ typedef typename Dune::FieldType<Vec>::type field_type;
+ /* constructor */
+ MatrixMulVector(const Mat & _A, const Vec & _v, int* _M) :
+ M(_M), A(_A), v(_v) {};
+ BlockExpr operator[] (int i) const {
+ M[0] = i;
+ return SubMatrixMulVector(A,v,M,*this);
+ }
+ int N() const { return -1; }; // { parent.begin().N(); }
+ private:
+ mutable int* M;
+ const Mat & A;
+ const Vec & v;
+ };
+
+ template <class K, int iN, int iM, class Mat, class Vec>
+ class MatrixMulVector< FieldMatrix<K,iN,iM>, Mat, Vec >
+ {
+ public:
+ typedef typename
+ BlockTypeN<MatrixMulVector<Mat,Mat,Vec>,
+ MyDepth<FieldMatrix<K,iN,iM>,Mat>::value-1>::type
+ ParentBlockType;
+ /* constructor */
+ MatrixMulVector(const Mat & _A, const Vec & _v, int* _M,
+ const ParentBlockType & _parent) :
+ parent(_parent), M(_M), A(_A), v(_v ) {};
+ K operator[] (int i) const {
+ K x=0;
+ M[MyDepth<FieldMatrix<K,iN,iM>,Mat>::value] = i;
+
+ FlatColIterator<const Mat> j(A[*M].begin(),M+1);
+ FlatColIterator<const Mat> endj(A[*M].end(),M+1);
+ for (; j!=endj; ++j)
+ {
+ x += (*j) * j.vectorentry(v);
+ }
+ return x;
+ }
+ int N() const { iN; };
+ const ParentBlockType & parent;
+ private:
+ mutable int* M;
+ const Mat & A;
+ const Vec & v;
+ };
+
+ template <class K, int iN, int iM>
+ class MatrixMulVector< FieldMatrix<K,iN,iM>, FieldMatrix<K,iN,iM>,
+ FieldVector<K,iM> >
+ {
+ public:
+ typedef FieldMatrix<K,iN,iM> Mat;
+ typedef FieldVector<K,iM> Vec;
+ MatrixMulVector(const Mat & _A, const Vec & _v) :
+ A(_A), v(_v ){};
+ K operator[] (int i) const {
+ K x=0;
+ typename Mat::ColIterator j = A[i].begin();
+ typename Mat::ColIterator endj = A[i].end();
+ for (; j!=endj; ++j)
+ {
+ x += (*j) * j.vectorentry(v);
+ }
+ return x;
+ }
+ int N() const { iN; };
+ private:
+ const Mat & A;
+ const Vec & v;
+ };
+
+ template <class M, class A, class B>
+ struct BlockExpression< MatrixMulVector< M, A, B > >
+ {
+ typedef Expression< MatrixMulVector<typename BlockType<M>::type,A,B> > type;
+ };
+
+ template <class K, int N, int M, class A, class B>
+ struct BlockExpression< MatrixMulVector< FieldMatrix<K,N,M>, A, B > >
+ {
+ typedef K type;
+ };
+
+ template<class K, int N, int M>
+ ExprTmpl::Expression<
+ MatrixMulVector<FieldMatrix<K,N,M>, FieldMatrix<K,N,M>, FieldVector<K,M> > >
+ operator * ( const FieldMatrix<K,N,M> & A, const FieldVector<K,M> & v )
+ {
+ return
+ ExprTmpl::Expression<
+ MatrixMulVector<FieldMatrix<K,N,M>, FieldMatrix<K,N,M>, FieldVector<K,M> > >
+ (
+ MatrixMulVector<FieldMatrix<K,N,M>, FieldMatrix<K,N,M>, FieldVector<K,M> >
+ (A, v)
+ );
+ }
+
+ // template<class BM, class BV>
+ // ExprTmpl::Expression<
+ // MatrixMulVector<BCRSMatrix<BM>, BCRSMatrix<BM>, BlockVector<BV> > >
+ // operator * ( const BCRSMatrix<BM> & A, const BlockVector<BV> & v )
+ // {
+ // static int indizes[20];
+ // return
+ // Expression<
+ // MatrixMulVector<BCRSMatrix<BM>, BCRSMatrix<BM>, BlockVector<BV> > >
+ // (
+ // MatrixMulVector<BCRSMatrix<BM>, BCRSMatrix<BM>, BlockVector<BV> >(A, v, indizes)
+ // );
+ // }
+
+ template<class M, class V>
+ ExprTmpl::Expression<
+ MatrixMulVector<Matrix<M>, Matrix<M>, Vector<V> > >
+ operator * ( const Matrix<M> & A, const Vector<V> & v )
+ {
+ static int indizes[20];
+ return
+ Expression<
+ MatrixMulVector<Matrix<M>, Matrix<M>, Vector<V> > >
+ (
+ MatrixMulVector<Matrix<M>, Matrix<M>, Vector<V> >(A, v, indizes)
+ );
+ }
+
+ template<class I>
+ struct ColIteratorType< Matrix<I> >
+ {
+ typedef typename I::ColIterator type;
+ };
+ template<class I>
+ struct ColIteratorType< const Matrix<I> >
+ {
+ typedef typename I::ConstColIterator type;
+ };
+
+ } // namespace ExprTmpl
+
+ template <class B, class A, class V>
+ struct FieldType< ExprTmpl::MatrixMulVector<B,A,V> >
+ {
+ typedef typename FieldType<V>::type type;
+ };
+ template <class I>
+ struct BlockType< ExprTmpl::Matrix<I> >
+ {
+ typedef typename BlockType<I>::type type;
+ };
+ template <class I>
+ struct FieldType< ExprTmpl::Matrix<I> >
+ {
+ typedef typename FieldType<I>::type type;
+ };
// OPERATORS
diff --git a/common/fmatrix.hh b/common/fmatrix.hh
index 778303cd5..9f9480741 100644
--- a/common/fmatrix.hh
+++ b/common/fmatrix.hh
@@ -155,6 +155,7 @@ namespace Dune {
fmmeta_umv<I-1>::template umv<Mat,X,Y,c>(A,x,y);
}
};
+
template<>
struct fmmeta_umv<0> {
template<class Mat, class X, class Y, int c>
@@ -502,8 +503,13 @@ namespace Dune {
Implementation of all members uses template meta programs where appropriate
*/
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ template<class K, int n, int m>
+ class FieldMatrix : ExprTmpl::Matrix< FieldMatrix<K,n,m> >
+#else
template<class K, int n, int m>
class FieldMatrix
+#endif
{
public:
// standard constructor and everything is sufficient ...
@@ -684,6 +690,9 @@ namespace Dune {
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
fmmeta_umv<n-1>::template umv<FieldMatrix,X,Y,m-1>(*this,x,y);
+ // for (int i=0; i<n; i++)
+ // for (int j=0; j<m; j++)
+ // y[i] += (*this)[i][j] * x[j];
}
//! y += A^T x
@@ -1417,6 +1426,19 @@ namespace Dune {
} // end namespace FMatrixHelp
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ template <class K, int N, int M>
+ struct BlockType< FieldMatrix<K,N,M> >
+ {
+ typedef K type;
+ };
+
+ template <class K, int N, int M>
+ struct FieldType< FieldMatrix<K,N,M> >
+ {
+ typedef K type;
+ };
+#endif // DUNE_EXPRESSIONTEMPLATES
/** @} end documentation */
diff --git a/common/test/.gitignore b/common/test/.gitignore
index 14efccf9b..cc769b2c2 100644
--- a/common/test/.gitignore
+++ b/common/test/.gitignore
@@ -18,4 +18,6 @@ poolallocatortest
gmon.out
gcdlcdtest
streamtest
-exprtmpl
\ No newline at end of file
+exprtmpl
+timing_xpr
+timing_old
\ No newline at end of file
diff --git a/common/test/Makefile.am b/common/test/Makefile.am
index b749532c5..db449a202 100644
--- a/common/test/Makefile.am
+++ b/common/test/Makefile.am
@@ -11,8 +11,9 @@ TESTS = $(TESTPROGS)
check_PROGRAMS = $(TESTPROGS)
# output code coverage
-AM_CXXFLAGS = -fprofile-arcs -ftest-coverage
-AM_LDFLAGS = $(top_builddir)/lib/libdune.la
+#AM_CXXFLAGS = -fprofile-arcs -ftest-coverage
+#AM_LDFLAGS = $(top_builddir)/lib/libdune.la
+AM_LDFLAGS = $(top_builddir)/common/libcommon.la
# define the programs
@@ -42,9 +43,13 @@ settest_SOURCES=settest.cc
gcdlcdtest_SOURCES = gcdlcdtest.cc
+bin_PROGRAMS = exprtmpl timing_old timing_xpr
exprtmpl_SOURCES = exprtmpl.cc
-exprtmpl_CXXFLAGS = -DDUNE_EXPRESSIONTEMPLATES
+exprtmpl_CXXFLAGS = -DDUNE_EXPRESSIONTEMPLATES -DNOPRINT
# -DDUNE_ISTL_WITH_CHECKING
-exprtmpl_LDFLAGS = $(top_builddir)/common/libcommon.la
+
+timing_old_SOURCES = timing.cc
+timing_xpr_SOURCES = timing.cc
+timing_xpr_CXXFLAGS = -DDUNE_EXPRESSIONTEMPLATES -g
include $(top_srcdir)/am/global-rules
diff --git a/common/test/exprtmpl.cc b/common/test/exprtmpl.cc
index 0952fb79a..50b02c052 100644
--- a/common/test/exprtmpl.cc
+++ b/common/test/exprtmpl.cc
@@ -7,11 +7,14 @@
int *M;
- fix RowBlock::N()
- remove second template parameter of FlatColIterator
+ - vectorentry -> exrpressionentry
+ - FlatColIterator<Matrix> does not work if Matrix is mutable
*/
#include <iostream>
#include <fstream>
#include <dune/common/fvector.hh>
+#include <dune/common/timer.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/io.hh>
#include <dune/common/iteratorfacades.hh>
@@ -58,11 +61,12 @@ void test_blockvector()
void test_blockblockvector()
{
- typedef Dune::FieldVector<double,2> VB;
+ const int bs = 2;
+ const int sz = 3;
+ typedef Dune::FieldVector<double,bs> 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;
@@ -77,6 +81,12 @@ void test_blockblockvector()
bbv[0] = Dune::ExprTmpl::Expression<RBV>(rbv1);
bbv[1].resize(bv2.N());
bbv[1] = Dune::ExprTmpl::Expression<RBV>(rbv2);
+
+ Dune::Timer stopwatch;
+ stopwatch.reset();
+ for (int i=0; i<100; i++) bbv *= 2;
+ std::cout << "Time bbv*2: " << stopwatch.elapsed() << std::endl;
+#ifndef NOPRINT
Dune::FlatIterator<BBV> fit(bbv.begin());
Dune::FlatIterator<BBV> fend(bbv.end());
int index = 0;
@@ -90,380 +100,66 @@ void test_blockblockvector()
printvector (std::cout, bv1, "bv1", "r");
printvector (std::cout, bv2, "bv1", "r");
printvector (std::cout, bbv, "bbv", "r");
-}
-
-namespace Dune {
-
- //! Type Traits for Vector::Iterator vs (const Vector)::ConstIterator
- template<class T>
- struct ColIteratorType
- {
- typedef typename T::ColIterator type;
- };
-
- template<class T>
- struct ColIteratorType<const T>
- {
- typedef typename T::ConstColIterator type;
- };
-
- //! Iterator class for flat sequential access to a nested Matrix Row
- template<class A, class V>
- class FlatColIterator :
- public ForwardIteratorFacade<FlatColIterator<A,V>,
- typename FieldType<A>::type,
- typename FieldType<A>::type&,
- int>
- {
- public:
- typedef typename ColIteratorType<A>::type ColBlockIterator;
- typedef std::ptrdiff_t DifferenceType;
- // typedef typename BlockIterator::DifferenceType DifferenceType;
- typedef typename BlockType<A>::type block_type;
- typedef typename BlockType<V>::type vblock_type;
- typedef typename FieldType<A>::type field_type;
- typedef FlatColIterator<block_type, vblock_type> SubBlockIterator;
- FlatColIterator(const ColBlockIterator & i, const int* _M) :
- M(_M), it(i),
- bit((*i)[(*M)].begin(), M+1),
- bend((*i)[(*M)].end(), M+1) {};
- void increment ()
- {
- ++bit;
- if (bit == bend)
- {
- ++it;
- bit = (*it)[(*M)].begin();
- bend = (*it)[(*M)].end();
- }
- }
- bool equals (const FlatColIterator & fit) const
- {
- return fit.it == it && fit.bit == bit;
- }
- const field_type& dereference() const
- {
- return *bit;
- }
- const field_type& vectorentry(const V & v) const
- {
- return bit.vectorentry(v[it.index()]);
- }
- //! return index
- DifferenceType index () const
- {
- return bit.index();
- }
- FlatColIterator operator = (const ColBlockIterator & _i)
- {
- it = _i;
- bit = (*it)[(*M)].begin();
- bend = (*it)[(*M)].end();
- return *this;
- }
- private:
- const int* M;
- ColBlockIterator it;
- SubBlockIterator bit;
- SubBlockIterator bend;
- };
-
- template<class K, int N, int M>
- class FlatColIterator<FieldMatrix<K,N,M>, FieldVector<K,M> > :
- public ForwardIteratorFacade<
- FlatColIterator< FieldMatrix<K,N,M>, FieldVector<K,M> >, K, K&, int>
- {
- public:
- typedef
- typename ColIteratorType< FieldMatrix<K,N,M> >::type ColBlockIterator;
- typedef std::ptrdiff_t DifferenceType;
- typedef K field_type;
- FlatColIterator(const ColBlockIterator & i, const int*) :
- it(i) {};
- void increment ()
- {
- ++it;
- }
- bool equals (const FlatColIterator & fit) const
- {
- return fit.it == it;
- }
- field_type& dereference() const
- {
- return *it;
- }
- const field_type& vectorentry(const FieldVector<K,M> & v) const
- {
- return v[it.index()];
- }
- //! return index
- DifferenceType index () const
- {
- return it.index();
- }
- FlatColIterator operator = (const ColBlockIterator & _i)
- {
- it = _i;
- return *this;
- }
- private:
- ColBlockIterator it;
- };
-
- template<class K, int N, int M>
- class FlatColIterator<const FieldMatrix<K,N,M>, const FieldVector<K,M> > :
- public ForwardIteratorFacade<
- FlatColIterator< const FieldMatrix<K,N,M>, const FieldVector<K,M> >,
- const K, const K&, int>
- {
- public:
- typedef
- typename ColIteratorType< const FieldMatrix<K,N,M> >::type ColBlockIterator;
- typedef std::ptrdiff_t DifferenceType;
- typedef const K field_type;
- FlatColIterator(const ColBlockIterator & i, const int*) :
- it(i) {};
- void increment ()
- {
- ++it;
- }
- bool equals (const FlatColIterator & fit) const
- {
- return fit.it == it;
- }
- field_type& dereference() const
- {
- return *it;
- }
- const field_type& vectorentry(const FieldVector<K,M> & v) const
- {
- return v[it.index()];
- }
- //! return index
- DifferenceType index () const
- {
- return it.index();
- }
- FlatColIterator operator = (const ColBlockIterator & _i)
- {
- it = _i;
- return *this;
- }
- private:
- ColBlockIterator it;
- };
-
- template<class M>
- struct NestedDepth
- {
- enum { value = NestedDepth<typename BlockType<M>::type>::value + 1 };
- };
-
- template<class K, int N, int M>
- struct NestedDepth< FieldMatrix<K,N,M> >
- {
- enum { value = 1 };
- };
-
- template<class Me, class M>
- struct MyDepth
- {
- enum { value = MyDepth<Me,typename BlockType<M>::type>::value+1 };
- };
-
- template<class Me>
- struct MyDepth<Me,Me>
- {
- enum { value = 0 };
- };
-
- /**
- B: BlockMatrix type -> indicated the current level.
- M: ,,global'' Matrix type
- V: ,,global'' Vector type
- */
- template <class B, class Mat, class Vec>
- class RowBlock
- {
- public:
- typedef RowBlock<typename BlockType<B>::type,Mat,Vec> SubRowBlock;
- typedef typename Dune::FieldType<Vec>::type field_type;
- RowBlock(const Mat & _A, const Vec & _v, int* _M) :
- M(_M), A(_A), v(_v) {};
- SubRowBlock operator[] (int i) const {
- M[MyDepth<B,Mat>::value] = i;
- return SubRowBlock(A,v,M);
- }
- int N() const { return -1; }; //r.begin()->N(); }
- private:
- mutable int* M;
- const Mat & A;
- const Vec & v;
- };
-
- // template <class K, int n, class M, class V>
- // class BlockRowSum< FieldVector<K,n>, M, V >
- template <class K, int iN, int iM, class Mat, class Vec>
- class RowBlock< FieldMatrix<K,iN,iM>, Mat, Vec >
- {
- public:
- RowBlock(const Mat & _A, const Vec & _v, int* _M) :
- M(_M), A(_A), v(_v) {};
- K operator[] (int i) const {
- K x=0;
- M[MyDepth<FieldMatrix<K,iN,iM>,Mat>::value] = i;
-
- FlatColIterator<const Mat, const Vec> j(A[*M].begin(),M+1);
- FlatColIterator<const Mat, const Vec> endj(A[*M].end(),M+1);
- for (; j!=endj; ++j)
- {
- x += (*j) * j.vectorentry(v);
- }
- return x;
- }
- int N() const { return -1; }; //r.begin()->N(); }
- private:
- mutable int* M;
- const Mat & A;
- const Vec & v;
- };
-
- template <class M, class V>
- class MV
- {
- public:
- MV(const M & _A, const V & _v) : A(_A), v(_v) {};
- RowBlock<typename BlockType<M>::type,M,V> operator[] (int i) const {
- indizes[0] = i;
- return RowBlock<typename BlockType<M>::type,M,V>(A, v, indizes);
- }
- int N() const { return A.N(); }
- private:
- mutable int indizes[256];
- 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 {
- std::cout << "ARGH\n";
- 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< RowBlock<B,A,V> >
- {
- typedef typename FieldType<V>::type type;
- };
-
- namespace ExprTmpl {
-
- template <class A, class B>
- struct BlockExpression< MV< BCRSMatrix<A>,BlockVector<B> > >
- {
- typedef ExprTmpl::Expression< RowBlock< A, BCRSMatrix<A>,BlockVector<B> > > type;
- };
-
-#if 0
- // !TODO!
- template <class A, class B>
- struct BlockExpression< RowBlock< A, BCRSMatrix<A>,BlockVector<B> > >
- {
- typedef ExprTmpl::Expression< RowBlock< typename BlockType<A>::type, RowBlock< A, BCRSMatrix<A>,BlockVector<B> >,BlockVector<B> > > type;
- };
#endif
+}
- template <class K, int N, int M, class A, class B>
- struct BlockExpression< RowBlock< FieldMatrix<K,N,M>, A, B > >
- {
- typedef K type;
- };
-
- template <class K, int N, int M>
- struct BlockExpression< MV< FieldMatrix<K,N,M>, FieldVector<K,M> > >
- {
- typedef K type;
- };
-
- } // NS ExpreTmpl
-
- template <class T>
- struct BlockType< BCRSMatrix<T> >
- {
- typedef T type;
- };
-
- template <class K, int N, int M>
- struct BlockType< FieldMatrix<K,N,M> >
- {
- typedef K type;
- };
-
- template <>
- struct BlockType< double >
- {
- typedef double type;
- };
-
- template <class T>
- struct FieldType< BCRSMatrix<T> >
- {
- typedef typename FieldType<T>::type type;
- };
-
- template <class K, int N, int M>
- struct FieldType< FieldMatrix<K,N,M> >
- {
- typedef K type;
- };
-
-} // NS Dune
-
+// namespace Dune {
+
+// namespace ExprTmpl {
+
+// template <class K, int iN, int iM>
+// class MatrixMulVector< FieldMatrix<K,iN,iM>,
+// BCRSMatrix< FieldMatrix<K,iN,iM> >,
+// BlockVector< FieldVector<K,iM> > >
+// {
+// public:
+// typedef BCRSMatrix< FieldMatrix<K,iN,iM> > Mat;
+// typedef BlockVector< FieldVector<K,iM> > Vec;
+// typedef typename
+// BlockTypeN<MatrixMulVector<Mat,Mat,Vec>,
+// MyDepth<FieldMatrix<K,iN,iM>,Mat>::value-1>::type
+// ParentBlockType;
+// /* constructor */
+// MatrixMulVector(const Mat & _A, const Vec & _v, int* _M,
+// const ParentBlockType & _parent) :
+// parent(_parent), M(_M), A(_A), v(_v )
+// {
+// int parent_i = M[0];
+// typename Mat::ConstColIterator it = A[parent_i].begin();
+// typename Mat::ConstColIterator end = A[parent_i].end();
+// tmp = 0;
+// for (; it!=end; ++it)
+// {
+// it->umv(tmp,v[it.index()]);
+// }
+// };
+// K operator[] (int i) const {
+// return tmp[i];
+// }
+// int N() const { iN; };
+// const ParentBlockType & parent;
+// private:
+// FieldVector<K,iN> tmp;
+// mutable int* M;
+// const Mat & A;
+// const Vec & v;
+// };
+
+// } // NS ExpreTmpl
+
+// } // NS Dune
+
+//template<int BlockSize, int N, int M>
+template<int BN, int BM, int N, int M>
void test_matrix()
{
- static const int BlockSize = 2;
+ std::cout << "test_matrix<" << BN << ", " << BM << ", "
+ << N << ", " << M << ">\n";
+
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::FieldVector<matvec_t,BN> LVB;
+ typedef Dune::FieldVector<matvec_t,BM> VB;
+ typedef Dune::FieldMatrix<matvec_t,BN,BM> MB;
typedef Dune::BlockVector<LVB> LeftVector;
typedef Dune::BlockVector<VB> Vector;
typedef Dune::BCRSMatrix<MB> Matrix;
@@ -476,20 +172,24 @@ void test_matrix()
// a += M * b
_M.umv(b,a);
+#ifndef NOPRINT
printmatrix (std::cout, _M, "Matrix", "r");
printvector (std::cout, a, "Vector", "r");
+#endif
// a = M * b
+#if 0
a = _M*b;
+#endif
+#ifndef NOPRINT
printvector (std::cout, a, "Vector", "r");
-
- int N = 4;
- int M = 5;
+#endif
Matrix A(N,M,Matrix::row_wise);
- Matrix::CreateIterator i=A.createbegin();
- Matrix::CreateIterator end=A.createend();
+ typename Matrix::CreateIterator i=A.createbegin();
+ typename Matrix::CreateIterator end=A.createend();
+ std::cout << "Building matrix structure\n";
// build up the matrix structure
int c=0;
for (; i!=end; ++i)
@@ -500,27 +200,31 @@ void test_matrix()
i.insert(M-1);
c++;
}
- A = 0.0;
+ std::cout << "...done\n";
+#ifndef NOPRINT
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 rit=A.begin();
- Matrix::Iterator rend=A.end();
+ std::cout << "Assembling matrix\n";
+ typename Matrix::Iterator rit=A.begin();
+ typename Matrix::Iterator rend=A.end();
for (; rit!=rend; ++rit)
{
- Matrix::ColIterator cit=rit->begin();
- Matrix::ColIterator cend=rit->end();
+ typename Matrix::ColIterator cit=rit->begin();
+ typename Matrix::ColIterator cend=rit->end();
for (; cit!=cend; ++cit)
{
// *rit = rit.index();
*cit = 10*cit.index()+rit.index();
}
}
+ std::cout << "...done\n";
printmatrix (std::cout, A, "Matrix", "r");
+#endif
LeftVector v(N);
LeftVector v2(N);
@@ -533,24 +237,28 @@ void test_matrix()
for (; fit!=fend; ++fit)
*fit=c++;
- std::cout << A.M() << " " << x.N() << " " << v.N() << std::endl;
-
+ Dune::Timer stopwatch;
+ stopwatch.reset();
A.umv(x,v);
+ std::cout << "Time umv: " << stopwatch.elapsed() << std::endl;
using namespace Dune;
+#ifndef NOPRINT
printvector (std::cout, x, "Vector X", "r");
printvector (std::cout, v, "Vector", "r");
- v=1;
- MV<Matrix,Vector> eximp (A,x);
- Dune::ExprTmpl::Expression< MV<Matrix,Vector> > ex ( eximp );
- v = ex;
- printvector (std::cout, v, "Vector", "r");
- v2 = A * x;
+#endif
+ v2 = 0;
+ stopwatch.reset();
+ v2 += A * x;
+ std::cout << "Time v2+=A*x: " << stopwatch.elapsed() << std::endl;
+#ifndef NOPRINT
printvector (std::cout, v2, "Vector2", "r");
-#if 0
+#endif
+
+#ifndef NOPRINT
int rowIndex[]={1};
- FlatColIterator<Matrix,Vector> it(A[2].begin(),rowIndex);
+ FlatColIterator<const Matrix> it(A[2].begin(),rowIndex);
for (int i=0; i<5; i++)
{
std::cout << *it << " ";
@@ -558,6 +266,7 @@ void test_matrix()
}
std::cout << std::endl;
#endif
+ std::cout << std::endl;
}
int main()
@@ -565,16 +274,23 @@ int main()
// Dune::dvverb.attach(std::cout);
try
{
- std::ofstream mylog("/dev/null");
- Dune::dvverb.attach(mylog);
- test_fvector();
- test_blockvector();
- test_blockblockvector();
- test_matrix();
- exit(0);
+ // test_fvector();
+ // test_blockvector();
+ // test_blockblockvector();
+ test_matrix<2,3,3,4>();
+#ifdef NOPRINT
+ test_matrix<3,6,400000,500000>();
+ test_matrix<6,3,400000,500000>();
+ test_matrix<30,60,4000,5000>();
+ test_matrix<150,150,500,4000>();
+ test_matrix<150,150,1000,2000>();
+#endif
+ // test_matrix<150,150,2000,1000>(); // fails in fmeta_something
+ // test_matrix<150,150,4000,500>(); // fails in fmeta_something
}
catch (Dune::Exception & e)
{
std::cout << e << std::endl;
}
+ exit(0);
}
diff --git a/common/test/timing.cc b/common/test/timing.cc
new file mode 100644
index 000000000..f9c07acf7
--- /dev/null
+++ b/common/test/timing.cc
@@ -0,0 +1,107 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <iostream>
+#include <fstream>
+#include <dune/common/fvector.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/io.hh>
+#include <dune/common/iteratorfacades.hh>
+
+#ifdef DUNE_EXPRESSIONTEMPLATES
+Indent INDENT;
+#endif
+
+template<int bs, int sz>
+void timing_vector()
+{
+ std::cout << "timing_vector<" << bs << ", " << sz << ">\n";
+ typedef Dune::FieldVector<double,bs> VB;
+ typedef Dune::BlockVector<VB> BV;
+ typedef Dune::BlockVector<BV> BBV;
+ BV bv1(sz), bv2(sz);
+ bv1 = 1;
+ bv2 = 0;
+ bv2[1][0]=1;
+ bv2[1][1]=2;
+
+ BBV bbv(2);
+ bbv[0].resize(bv1.N());
+ bbv[1].resize(bv2.N());
+
+ Dune::Timer stopwatch;
+ stopwatch.reset();
+ for (int i=0; i<100; i++) bbv *= 2;
+ std::cout << "Time [bbv*=2] " << stopwatch.elapsed() << std::endl;
+}
+
+//template<int BlockSize, int N, int M>
+template<int BN, int BM, int N, int M>
+void timing_matrix()
+{
+ std::cout << "timing_matrix<" << BN << ", " << BM << ", "
+ << N << ", " << M << ">\n";
+
+ typedef double matvec_t;
+ typedef Dune::FieldVector<matvec_t,BN> LVB;
+ typedef Dune::FieldVector<matvec_t,BM> VB;
+ typedef Dune::FieldMatrix<matvec_t,BN,BM> MB;
+ typedef Dune::BlockVector<LVB> LeftVector;
+ typedef Dune::BlockVector<VB> Vector;
+ typedef Dune::BCRSMatrix<MB> Matrix;
+
+ Matrix A(N,M,Matrix::row_wise);
+ typename Matrix::CreateIterator i=A.createbegin();
+ typename Matrix::CreateIterator end=A.createend();
+ std::cout << "Building matrix structure\n";
+ // 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++;
+ }
+ std::cout << "...done\n";
+
+ LeftVector v(N);
+ v = 0;
+ Vector x(M);
+ x = 1;
+
+ Dune::Timer stopwatch;
+ stopwatch.reset();
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ v += A * x;
+#else
+ A.umv(x,v);
+#endif
+ std::cout << "Time [v+=A*x] " << stopwatch.elapsed() << std::endl;
+
+ std::cout << std::endl;
+}
+
+int main ()
+{
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ std::cout << "Expression Templates\n";
+#else
+ std::cout << "Old Version\n";
+#endif
+
+ // timing_vector<1,1000000>();
+ // timing_vector<10,100000>();
+ // timing_vector<100,10000>();
+
+ // timing_matrix<150,150,500,4000>();
+ // timing_matrix<150,150,1000,2000>();
+ timing_matrix<1,18,400000,500000>();
+ // timing_matrix<6,3,400000,500000>();
+ // timing_matrix<3,6,400000,500000>();
+ // timing_matrix<18,1,400000,500000>();
+ // timing_matrix<50,50,9000,10000>();
+
+ std::cout << std::endl;
+}
diff --git a/istl/bcrsmatrix.hh b/istl/bcrsmatrix.hh
index c88efcf54..22e569dc7 100644
--- a/istl/bcrsmatrix.hh
+++ b/istl/bcrsmatrix.hh
@@ -52,8 +52,13 @@ namespace Dune {
Setting the compile time switch DUNE_ISTL_WITH_CHECKING
enables error checking.
*/
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ template<class B, class A>
+ class BCRSMatrix : public ExprTmpl::Matrix< BCRSMatrix<B,A> >
+#else
template<class B, class A=ISTLAllocator>
class BCRSMatrix
+#endif
{
public:
@@ -1282,6 +1287,25 @@ namespace Dune {
};
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ template <class B, class A>
+ struct FieldType< BCRSMatrix<B,A> >
+ {
+ typedef typename FieldType<B>::type type;
+ };
+
+ template <class B, class A>
+ struct BlockType< BCRSMatrix<B,A> >
+ {
+ typedef B type;
+ };
+ template <class B, class A>
+ struct RowType< BCRSMatrix<B,A> >
+ {
+ typedef CompressedBlockVectorWindow<B,A> type;
+ };
+#endif
+
/** @} end documentation */
} // end namespace
diff --git a/istl/bvector.hh b/istl/bvector.hh
index 07da9b6a8..9b6eec2be 100644
--- a/istl/bvector.hh
+++ b/istl/bvector.hh
@@ -222,10 +222,13 @@ namespace Dune {
Setting the compile time switch DUNE_ISTL_WITH_CHECKING
enables error checking.
*/
+#ifdef DUNE_EXPRESSIONTEMPLATES
+ template<class B, class A>
+ class BlockVector : public block_vector_unmanaged<B,A> ,
+ public Dune::ExprTmpl::Vector< Dune::BlockVector<B,A> >
+#else
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:
--
GitLab