Skip to content
Snippets Groups Projects
Commit 0378a976 authored by Christian Engwer's avatar Christian Engwer
Browse files

remove expression template support, until it is properly implemented 

[[Imported from SVN: r6596]]
parent b17f00a0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1 addition and 1534 deletions
configure.ac
+ 0
− 1
View file @ 0378a976
......@@ -25,7 +25,6 @@ AC_CONFIG_FILES([Makefile
dune/Makefile
dune/common/Makefile
dune/common/test/Makefile
dune/common/exprtmpl/Makefile
dune/common/parallel/Makefile
dune/common/parallel/test/Makefile
doc/Makefile
......
dune/common/Makefile.am
+ 1
− 8
View file @ 0378a976
# $Id$
SUBDIRS = . test exprtmpl parallel
SUBDIRS = . test parallel
# the standard debug streams are put into the libdune
noinst_LTLIBRARIES = libcommon.la
......@@ -83,11 +83,4 @@ commoninclude_HEADERS = \
unused.hh \
version.hh
if EXPRESSIONTEMPLATES
commoninclude_HEADERS += exprtmpl.hh exprtmpl/scalar.inc exprtmpl/exprexpr.inc
libcommon_la_SOURCES += exprtmpl.cc
else
sourcescheck_NOSOURCES = exprtmpl.hh exprtmpl.cc
endif
include $(top_srcdir)/am/global-rules
dune/common/exprtmpl.cc deleted 100644 → 0
+ 0
− 6
View file @ b17f00a0
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "exprtmpl.hh"
/** This is only needed for debugging and will be removed in the future */
Indent INDENT;
dune/common/exprtmpl.hh deleted 100644 → 0
+ 0
− 1105
View file @ b17f00a0
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#if ! defined DUNE_EXPRTMPL_HH && DUNE_EXPRESSIONTEMPLATES
#define DUNE_EXPRTMPL_HH
/*! @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>
#include <cmath>
#include <complex>
#include "iteratorfacades.hh"
#include "static_assert.hh"
#include "stdstreams.hh"
struct Indent
{
int i;
Indent() {
i = 0;
}
void operator ++ ()
{
#ifdef DUNE_VVERBOSE
i += 3;
#endif
};
void operator -- ()
{
#ifdef DUNE_VVERBOSE
i -= 3;
#endif
};
};
extern Indent INDENT;
inline std::ostream & operator << (std::ostream & s, const Indent & i)
{
#ifdef DUNE_VVERBOSE
for (int n = 0; n < i.i; n++) s << " ";
#endif
return s;
}
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
*/
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 FieldVector
template <class K, int N>
struct FieldType< FieldVector<K,N> >
{
typedef K type;
};
// BlockType specialization for FieldVector
template <class K, int N>
struct BlockType< FieldVector<K,N> >
{
typedef K type;
};
// FieldType specialization for const T
template<class T>
struct FieldType<const T>
{
typedef const typename FieldType<T>::type type;
};
// BlockType specialization for const T
template<class T>
struct BlockType<const T>
{
typedef const typename BlockType<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;
template <class I> class Matrix;
/**
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 <class K>
struct isEndOfExpressionRecusion< FieldVector<K,1> >
{
enum { value=true };
};
template <class K>
struct isEndOfExpressionRecusion< std::complex<K> >
{
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];
}
size_t N() const { return ex.N(); }
double one_norm() const { return eval_one_norm(*this); }
double one_norm_real() const { return eval_one_norm_real(*this); }
double two_norm() const { return sqrt(eval_two_norm2(*this)); }
double two_norm2() const { return eval_two_norm2(*this); }
double infinity_norm() const { return eval_infinity_norm(*this); }
double infinity_norm_real() const { return eval_infinity_norm_real(*this); }
private:
Ex ex;
};
/**
Vector Base Class for Expression 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
size_t N() const {
return asImp().N();
}
double one_norm() const { return eval_one_norm(*this); }
double one_norm_real() const { return eval_one_norm_real(*this); }
double two_norm() const { return sqrt(eval_two_norm2(*this)); }
double two_norm2() const { return eval_two_norm2(*this); }
double infinity_norm() const { return eval_infinity_norm(*this); }
double infinity_norm_real() const { return eval_infinity_norm_real(*this); }
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
#ifdef DUNE_VVERBOSE
Dune::dvverb << INDENT << "Assign Vector from Expression\n";
#endif
++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
#ifdef DUNE_VVERBOSE
Dune::dvverb << INDENT << "Assign Vector from Vector\n";
#endif
++INDENT;
for (int i=0; i<N(); ++i) { asImp()[i] = v[i]; }
--INDENT;
return asImp();
}
/*
I& assignFrom(const Vector<block_type> & x) {
#ifdef DUNE_VVERBOSE
Dune::dvverb << INDENT << "Assign Vector block_type\n";
#endif
++INDENT;
for (int i=0; i < asImp().N(); i++) asImp()[i] = x;
--INDENT;
return asImp();
}
*/
I& assignFrom(field_type x) {
#ifdef DUNE_VVERBOSE
Dune::dvverb << INDENT << "Assign Vector from field_type\n";
#endif
++INDENT;
for (size_t i=0; i<N(); ++i) { asImp()[i] = x; }
--INDENT;
return asImp();
}
template <class E> Vector<I>& operator+=(const Expression<E>& x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] += x[i];
return asImp();
}
template <class V> Vector<I>& operator+=(const Vector<V>& x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] += x[i];
return asImp();
}
template <class E> Vector<I>& operator-=(const Expression<E>& x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] -= x[i];
return asImp();
}
template <class V> Vector<I>& operator-=(const Vector<V>& x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] -= x[i];
return asImp();
}
Vector<I>& operator+=(field_type x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] += x;
return asImp();
}
Vector<I>& operator-=(field_type x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] -= x;
return asImp();
}
Vector<I>& operator*=(field_type x) {
for (size_t i=0; i < asImp().N(); i++) asImp()[i] *= x;
return asImp();
}
Vector<I>& operator/=(field_type x) {
for (size_t 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 {
#ifdef DUNE_VVERBOSE
Dune::dvverb << INDENT << "ConstRef->dereference " << v[i] << std::endl;
#endif
return BlockExprImp(v[i]);
}
size_t N() const { return v.N(); };
double one_norm() const { return eval_one_norm(*this); }
double one_norm_real() const { return eval_one_norm_real(*this); }
double two_norm() const { return sqrt(eval_two_norm2(*this)); }
double two_norm2() const { return eval_two_norm2(*this); }
double infinity_norm() const { return eval_infinity_norm(*this); }
double infinity_norm_real() const { return eval_infinity_norm_real(*this); }
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;
};
/**
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
size_t 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);
}
size_t 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);
}
size_t 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;
}
size_t N() const { return 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;
}
size_t N() const { return 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
/* 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"
/* one norm (sum over absolute values of entries) */
#define NORM eval_one_norm
#define NORM_CODE \
{ \
typename FieldType<A>::type val=0; \
Dune::dvverb << INDENT << "Infinity Norm of Expression\n"; \
++INDENT; \
for (size_t i=0; i<a.N(); ++i) { val += eval_one_norm(a[i]); } \
--INDENT; \
return val; \
}
#define VAL_CODE { return std::abs(a); }
#include "exprtmpl/norm.inc"
template<class K>
inline K eval_one_norm (const std::complex<K>& c)
{
sqrt(c.real()*c.real() + c.imag()*c.imag());
}
template <class A>
typename FieldType<A>::type
one_norm (const A & a)
{
return eval_one_norm(a);
}
/* simplified one norm (uses Manhattan norm for complex values) */
#define NORM eval_one_norm_real
#define NORM_CODE \
{ \
typename FieldType<A>::type val=0; \
Dune::dvverb << INDENT << "Infinity Norm of Expression\n"; \
++INDENT; \
for (size_t i=0; i<a.N(); ++i) { val += eval_one_norm_real(a[i]); } \
--INDENT; \
return val; \
}
#define VAL_CODE { return std::abs(a); }
#include "exprtmpl/norm.inc"
template<class K>
inline K eval_one_norm_real (const std::complex<K>& c)
{
return eval_one_norm_real(c.real()) + eval_one_norm_real(c.imag());
}
template <class A>
typename FieldType<A>::type
one_norm_real (const A & a)
{
return eval_one_norm_real(a);
}
/* two norm sqrt(sum over squared values of entries) */
template <class A>
typename FieldType<A>::type
two_norm (const A & a)
{
return sqrt(eval_two_norm2(a));
}
/* sqare of two norm (sum over squared values of entries), need for block recursion */
#define NORM eval_two_norm2
#define NORM_CODE \
{ \
typename FieldType<A>::type val=0; \
Dune::dvverb << INDENT << "Infinity Norm of Expression\n"; \
++INDENT; \
for (size_t i=0; i<a.N(); ++i) { val += eval_two_norm2(a[i]); } \
--INDENT; \
return val; \
}
#define VAL_CODE { return a*a; }
#include "exprtmpl/norm.inc"
template<class K>
inline K eval_two_norm2 (const std::complex<K>& c)
{
return c.real()*c.real() + c.imag()*c.imag();
}
template <class A>
typename FieldType<A>::type
two_norm2 (const A & a)
{
return eval_two_norm2(a);
}
/* infinity norm (maximum of absolute values of entries) */
#define NORM eval_infinity_norm
#define NORM_CODE { \
typename FieldType<A>::type val=0; \
Dune::dvverb << INDENT << "Infinity Norm of Expression\n"; \
++INDENT; \
for (size_t i=0; i<a.N(); ++i) { val = std::max(val,eval_infinity_norm(a[i])); } \
--INDENT; \
return val; \
}
#define VAL_CODE { return a; }
#include "exprtmpl/norm.inc"
template <class A>
typename FieldType<A>::type
infinity_norm (const A & a)
{
return eval_infinity_norm(a);
}
/* simplified infinity norm (uses Manhattan norm for complex values) */
#define NORM eval_infinity_norm_real
#define NORM_CODE { \
typename FieldType<A>::type val=0; \
Dune::dvverb << INDENT << "Infinity Norm of Expression\n"; \
++INDENT; \
for (size_t i=0; i<a.N(); ++i) { val = std::max(val,eval_infinity_norm(a[i])); } \
--INDENT; \
return val; \
}
#define VAL_CODE { return std::abs(a); }
#include "exprtmpl/norm.inc"
template<class K>
inline K eval_infinity_norm_real (const std::complex<K>& c)
{
return eval_one_norm_real(c.real()) + eval_one_norm_real(c.imag());
}
template <class A>
typename FieldType<A>::type
infinity_norm_real (const A & a)
{
return eval_infinity_norm(a);
}
/* vector * vector */
namespace ExprTmpl {
// Vector * Vector
template <class A>
typename FieldType<A>::type
operator * (const Vector<A> & a, const Vector<A> & b)
{
assert(a.N() == b.N());
typename FieldType<A>::type x = 0;
for (size_t i=0; i<a.N(); i++)
x = a[i] * b[i];
return x;
}
// Expression * Vector
template <class A, class B>
typename FieldType<A>::type
operator * (const Vector<A> & a, const Expression<B> & b)
{
dune_static_assert((is_same<FieldType<A>,FieldType<B> >::value),
"Field types of both operands must match!");
assert(a.N() == b.N());
typename FieldType<A>::type x = 0;
for (size_t i=0; i<a.N(); i++)
x = a[i] * b[i];
return x;
}
// Vector * Expression
template <class A, class B>
typename FieldType<A>::type
operator * (const Expression<A> & a, const Vector<B> & b)
{
dune_static_assert((is_same<FieldType<A>,FieldType<B> >::value),
"Field types of both operands must match!");
assert(a.N() == b.N());
typename FieldType<A>::type x = 0;
for (size_t i=0; i<a.N(); i++)
x = a[i] * b[i];
return x;
}
} // namespace ExprTmpl
} // namespace Dune
#endif // DUNE_EXPRTMPL_HH
dune/common/exprtmpl/Makefile.am deleted 100644 → 0
+ 0
− 7
View file @ b17f00a0
# $Id: $
SOURCES = exprexpr.inc norm.inc scalar.inc
exprtmpldir = $(includedir)/dune/common/exprtmpl
exprtmpl_HEADERS = $(SOURCES)
include $(top_srcdir)/am/global-rules
dune/common/exprtmpl/exprexpr.inc deleted 100644 → 0
+ 0
− 136
View file @ b17f00a0
// -*- 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);
}
// Expression op FieldVector<K,1>
template <class A, class K>
K operator OP (const Expression<A> & a,
const Expression< ConstRef< FieldVector<K,1> > >& b)
{
return a OP b[0];
}
// 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);
}
// FieldVector<K,1> op Expression
template <class K, class B>
K operator OP (const Expression< ConstRef< FieldVector<K,1> > > & a,
const Expression<B> & b)
{
return a[0] OP b;
}
// 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);
}
// FieldVector<K,1> op FieldVector<K,1>
template <class K>
K operator OP (const Expression< ConstRef< FieldVector<K,1> > > & a,
const Expression< ConstRef< FieldVector<K,1> > > & b)
{
return a[0] OP b[0];
}
// 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
dune/common/exprtmpl/norm.inc deleted 100644 → 0
+ 0
− 74
View file @ b17f00a0
// -*- C++ -*-
#ifndef NORM
#error NORM undefined
#endif
#ifndef NORM_CODE
#error NORM_CODE undefined
#endif
#ifndef VAL_CODE
#error VAL_CODE undefined
#endif
template <class A>
typename FieldType<A>::type
NORM (const ExprTmpl::Expression<A> & a)
{
NORM_CODE
}
template <class A>
typename FieldType<A>::type
NORM (const ExprTmpl::Vector<A> & a)
{
NORM_CODE
}
/*
template <class A>
double
ExprTmpl::Expression<A>::NORM () const
{
const ExprTmpl::Expression<A> & a = *this;
NORM_CODE
}
template <class A>
double
ExprTmpl::Vector<A>::NORM () const
{
const ExprTmpl::Vector<A> & a = *this;
NORM_CODE
}
template <class A>
double
ExprTmpl::ConstRef<A>::NORM () const
{
const ExprTmpl::ConstRef<A> & a = *this;
NORM_CODE
}
*/
inline double
NORM (const double & a)
{
VAL_CODE
}
inline float
NORM (const float & a)
{
VAL_CODE
}
inline int
NORM (const int & a)
{
VAL_CODE
}
#undef NORM
#undef NORM_CODE
#undef VAL_CODE
dune/common/exprtmpl/scalar.inc deleted 100644 → 0
+ 0
− 185
View file @ b17f00a0
// -*- 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;
};
// Scalar op FieldVector<K,1>
template <class K>
K operator OP (const Expression< ConstRef< FieldVector<K,1> > > & a,
const K & lambda)
{
return a[0] OP 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 FieldVector<K,1>
template <class K>
K operator OP (const K & lambda,
const Expression< ConstRef< FieldVector<K,1> > > & a)
{
return lambda OP a[0];
}
// 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 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
m4/Makefile.am
+ 0
− 1
View file @ 0378a976
......@@ -28,7 +28,6 @@ ALLM4S = \
dune_deprecated.m4 \
dune_deprecated_cppflags.m4 \
dune_docu.m4 \
dune_exprtmpl.m4 \
dune_fieldvector_size_is_method.m4 \
dune_linkcxx.m4 \
dune_mpi.m4 \
......
m4/dune_common.m4
+ 0
− 1
View file @ 0378a976
......@@ -29,7 +29,6 @@ AC_DEFUN([DUNE_COMMON_CHECKS],
AC_REQUIRE([DUNE_SET_MINIMAL_DEBUG_LEVEL])
AC_REQUIRE([DUNE_PATH_XDR])
AC_REQUIRE([DUNE_MPI])
AC_REQUIRE([DUNE_EXPRTMPL])
AC_REQUIRE([DUNE_TR1_HEADERS])
dnl Create the flag --enable-fieldvector-size-is-method
......
m4/dune_exprtmpl.m4 deleted 100644 → 0
+ 0
− 10
View file @ b17f00a0
AC_DEFUN([DUNE_EXPRTMPL],[
# enable experimental features
AC_ARG_ENABLE(expressiontemplates,
AS_HELP_STRING([--enable-expressiontemplates],[enable experimental expression templates in dune]))
AS_IF([test "x$enable_expressiontemplates" = "xyes"],[
AC_DEFINE([DUNE_EXPRESSIONTEMPLATES], [1],
[Define to 1 if the experimental expression templates should be used])])
AM_CONDITIONAL([EXPRESSIONTEMPLATES],
[test "x$enable_expressiontemplates" = "xyes"])
])
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment