Generic "expression template" library?

D

Daniel Pitts

Are there any popular expression template libraries that support
all/most arithmetic operators, and that allow "compound" expressions to
be created?

for example...
template<typename A, typename B>
struct Multiply;
template<typename A, typename B>
struct Add;


x = Multiply<double, double>(1,3) + 4;

where x might be of type Add<Multiply<double, double>, int>.

I was thinking of creating my own, but if there is an existing one that
handles all the odd cases (such as avoiding infinite recursion), all the
better.
 
J

Juha Nieminen

Daniel said:
x = Multiply<double, double>(1,3) + 4;

Out of curiosity, exactly how would that be different from:

x = 1.0 * 3.0 + 4:

which a compiler will also be able to calculate at compile time (and
which is much easier to read)?
 
D

Daniel Pitts

Juha said:
Out of curiosity, exactly how would that be different from:

x = 1.0 * 3.0 + 4:

which a compiler will also be able to calculate at compile time (and
which is much easier to read)?
Because in the end I won't be dealing with doubles. It may be Vectors
(not std::vector), Matrices, Colors, WhatEverYouThinkOf, etc...

It may also be that A a; B b; C c(a*b), where A, B, and C are not the
same type. I need at least on expression template to handle that, but it
would be nice if the expression template could be arbitrary depth.
 
J

James Kanze

Are there any popular expression template libraries that
support all/most arithmetic operators, and that allow
"compound" expressions to be created?
for example...
template<typename A, typename B>
struct Multiply;
template<typename A, typename B>
struct Add;
x = Multiply<double, double>(1,3) + 4;
where x might be of type Add<Multiply<double, double>, int>.
I was thinking of creating my own, but if there is an existing
one that handles all the odd cases (such as avoiding infinite
recursion), all the better.

http://www.oonumerics.org/blitz/? I'm not too familiar with
this library, as I don't work in numerical processing, but from
what I understand, it is the reference for this sort of thing.
(The original author, Todd Veldhuizen, is one of the inventors
of template metaprogramming.)

I do think, however, that it expects overloaded operators,
rather than named obejcts (like Multiply).
 
J

Juha Nieminen

Daniel said:
Because in the end I won't be dealing with doubles. It may be Vectors
(not std::vector), Matrices, Colors, WhatEverYouThinkOf, etc...

I still fail to see the need for a templated version of the basic
operators. For example, what's the relevant difference between:

x = Multiply<Color, Color>(c1, c2);

and:

x = c1 * c2;

?

If I understood correctly, the template would effectively simply do
the latter. I completely fail to see any advantage in doing it in such
an obfuscated way.
It may also be that A a; B b; C c(a*b), where A, B, and C are not the
same type.

Exactly how does template arithmetic help that? Inside, the template
will just do the a*b anyways.
 
N

Noah Roberts

Juha said:
I still fail to see the need for a templated version of the basic
operators. For example, what's the relevant difference between:

x = Multiply<Color, Color>(c1, c2);

and:

x = c1 * c2;

?

Consider the problem of multiplying 3 3x3 matrices. A straight forward
method of doing so might be to implement operator * like so:

matrix_type operator * (matrix_type const& a, matrix_type const& b)
{
matrix_type ret;
for all i,j do:
ret(i,j) = a(i,j) * b(j,i);
return ret;
}

or something like that. Important thing to note here is that you've a
copying operator here. Any time * appears in an expression there's a
temporary.

The expression template version will instead do something like this:

matrix_multiplication operator * (mt const& a, mt const& b)
{
return matrix_multiplication(a,b);
}

matrix_multiplication then has a conversion operator to convert to
matrix_type that does each of the necessary operations for each element
accessed.

See
http://ubiety.uwaterloo.ca/~tveldhui/papers/Expression-Templates/exprtmpl.html
for a more complete example and rationale.
 
J

James Kanze

Consider the problem of multiplying 3 3x3 matrices. A
straight forward method of doing so might be to implement
operator * like so:
matrix_type operator * (matrix_type const& a, matrix_type const& b)
{
matrix_type ret;
for all i,j do:
ret(i,j) = a(i,j) * b(j,i);
return ret;
}
or something like that.
Important thing to note here is that you've a copying operator
here. Any time * appears in an expression there's a
temporary.
The expression template version will instead do something like this:
matrix_multiplication operator * (mt const& a, mt const& b)
{
return matrix_multiplication(a,b);
}

Of course, but you're still using operator *, and not
Multiple<double,double>().

(FWIW: this technique was used even before the language had
templates, using virtual functions. Since the compiler knew the
actual types involved in the expression, it could still inline
the functions.)
matrix_multiplication then has a conversion operator to
convert to matrix_type that does each of the necessary
operations for each element accessed.

Or there is an overload of the assignment operator and a
constructor which takes matrix_multiplication. Of course,
normally, the overloads will be templates, so that it can also
handle matrix_add, matrix_subtract, etc. (Or in earlier days, it
took a const reference to the base class of all of these
classes).
 
J

Juha Nieminen

Noah said:
Important thing to note here is that you've a
copying operator here. Any time * appears in an expression there's a
temporary.

The expression template version will instead do something like this:

matrix_multiplication operator * (mt const& a, mt const& b)
{
return matrix_multiplication(a,b);
}

matrix_multiplication then has a conversion operator to convert to
matrix_type that does each of the necessary operations for each element
accessed.

I'm not sure if I understood correctly, but are you saying that with
the templates it's basically possible to perform lazy evaluation of the
operations involved which, with complex expressions, will minimize the
amount of temporaries created?

(Btw, the move semantics of the upcoming standard will also offer a
tool to lessen the impact of temporaries.)
 
K

Kai-Uwe Bux

Juha said:
I'm not sure if I understood correctly, but are you saying that with
the templates it's basically possible to perform lazy evaluation of the
operations involved which, with complex expressions, will minimize the
amount of temporaries created?

To some degree. The matrix multiplication example of square matrices is
interesting. If you had the expression

trace( A * B )

then an expression template would indeed only evaluate the diagonal entries
of the product A*B, resulting in quadratic runtime. If you had

( A * B ) * C

on the other hand, a simple minded expression template would need to
re-evaluate each entry of A*B anytime it is needed in evaluating the triple
product (the entry needs to be re-evaluated since it is not stored in a
temporary). This would result in bi-quadratic runtime whereas introducing a
temporary would result in cubic run-time. One common trick to speed up
things is therefore to force the creation of temporaries like

Matrix( A*B ) * C

In principle, expression templates could detect whether they appear inside
an expression of the first kind or an expression of the second kind.
However, that is somewhat tricky.


Best

Kai-Uwe Bux
 
J

James Kanze

I'm not sure if I understood correctly, but are you saying
that with the templates it's basically possible to perform
lazy evaluation of the operations involved which, with complex
expressions, will minimize the amount of temporaries created?

with templates or polymorphic classes, sort of. The number of
temporaries stays the same, but the temporaries are of very
small types (e.g. two pointers to already existing objects),
rather than enormous matrixes. The maxtrix_add class might look
something like:

template< typename T, typename LHS, typename RSH >
class maxtrix_add
{
public:
maxtrix_add( LHS const& lhs, RHS const& rhs )
: lhs( lhs )
, rhs( rhs )
{
assert( lhs.getX() == rhs.getX() ) ;
assert( lhs.getY() == rhs.getY() ) ;
}

T get( size_t i, size_t j ) const
{
return lhs.get( i, j ) + rhs.get( i, j ) ;
}
private:
LHS const& lhs ;
RHS const& rhs ;
} ;

template< typename T, typename LHS, typename RHS >
matrix_add< T, LHS, RHS >
operator+(
LHS const& lhs,
RHS const& rhs )
{
return matrix_add< T, LHS, RHS >( lhs, rhs ) ;
}

and matrix_multiply:

template< typename T, typename LHS, typename RHS >
class matrix_multiply
{
public:
matrix_multiply( LHS& lhs, RHS& rhs )
: lhs( lhs )
, rhs( rhs )
{
assert( lhs.getX() == rhs.getY() ) ;
}

T get( size_t i, size_t j ) const
{
T result = 0 ;
for ( size_t k = 0 ; k != lhs.getX() ; ++ k ) {
result += lhs.get( i, k ) * rhs.get( k, j ) ;
}
return result ;
}

private:
LHS const& lhs ;
RHS const& rhs ;
} ;

template< typename T, typename LHS, typename RHS >
matrix_multipy< T, LHS, RHS >
operator*(
LHS const& lhs,
RHS const& rhs )
{
return matrix_multipy< T, LHS, RHS >( lhs, rhs ) ;
}

The intermediate matrixes in a complicated expression are never
constructed. (This is obviously just an example. You'ld have
to do something to ensure that the template operators are only
found when you're actually dealing with matrices---when I
learned this idiom, we used virtual functions, with the
operators taking references to the base class, and there was no
problem.)
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

Forum statistics

Threads
473,773
Messages
2,569,594
Members
45,123
Latest member
Layne6498
Top