Re: [lmi] overview of C++ expression template libraries

[Greg Chicares]

On 2005-8-31 12:33 UTC, Vadim Zeitlin wrote:
> On Tue, 30 Aug 2005 20:17:42 +0000 Greg Chicares <address@hidden> wrote:
> 
> GC> Furthermore, we should consider what error-detection facilities each final
> GC> candidate provides. Trying to add two vectors of different lengths is an
> GC> error that I'd like to be detected automatically, even if only in a
> GC> nondefault 'debug' mode.
> 
>  I have some doubts that many libraries do this by default but I'll check.
> 
> GC> Let me sketch out a test suite.

I grepped the lmi sources for 'std::transform' and added inline comments,
marked with "ET !!", that suggest translations to the sort of notation
we'd like to be able to use. In a couple of places, I could hardly
understand the C++ I had written, but saw that I had described a simple
APL equivalent in a comment, and then everything became clear.

Then I sifted through all those new comments and narrowed them down to a
reasonably brief test suite:

  data:
    vectors of double: vd0, vd1, vd2
    vectors of int: vi0
    scalar double: s0
    unary and binary functions:
      double uf(double);
      double bf(double, double);

  expressions to test:
    vd0 = -1.0 + 1.0 / vd1;
    vd0 -= vd0 - vd1 * vd2 * vd3;
    vd0 = uf(greater_of(s0, vd1 - vd1[0]) - vi0);
    vd0 = bf((greater_of(0.0, lesser_of(1.0, vd1 * vd2)), s0);

I think it would be sufficient to measure how long those expressions take
to run, and to compile, with various libraries. It's good enough if all
vectors have length 50: it's almost always in [20, 100] for lmi.

> GC> I uploaded 'vector_test.cpp' to cvs the
> GC> other day. It's just the beginning of an investigation I started years ago
> GC> but never finished.

In cleaning it up, I found it to be of little use for the present purpose...

> I think I can adapt it to help test typical lmi
> operations with various techniques.

...so I committed a new 'expression_template_0_test.cpp' in case you want
to use it to investigate other libraries, though of course you might find
a better way to set that up.

>  Ok, I think it's tbe best we can do. To recapitulate, I'll then
> reimplement the tests using:
> 
> 1. PETE if I can make it work as easily as you say
> 2. uBLAS
> 3. boost::lambda
> 4. std::valarray
> 5. maybe something hand-coded if I can do it quickly (< day?)
> 
> I'll then show you the code and we'll also compare the performances of
> different approaches. Does this sound good?

Yes, thanks. Notes on two of those items:

> 3. boost::lambda

Please take a look at 'expression_template_0_test.cpp' first. The comments
preceding mete_stl_smart() explain an obstacle I ran into. Very often we
need to write expressions that I don't see how to write with the combination
of std::transform and boost::lambda. That combination supports
  vector0 = [arbitrarily complicated expression]
well because that's boost::lambda's strength, but doesn't seem to support
  vector0 += [rather simple expression]
  vector0 -= [rather simple expression]
  vector0 *= [rather simple expression]
  ...
or
  vector0 = vector1 + vector2 + vector3 + vector4;
because those are limitations of std::transform.

> 5. maybe something hand-coded if I can do it quickly (< day?)

I'd suggest postponing that until we see what uBLAS, PETE, and std::valarray
can do. We might find that we really need to create something new, and I got
some interesting ideas about that from rereading lots of lmi code, e.g.

 - Maybe we could use boost::lambda as a deferred-evaluation expression
   engine, and just write an array-transversal wrapper on top of it.
   That would make it easy to do bounds checking, BTW.

 - Because APL really is the natural notation for much of this work,
   we might consider implementing a subset of it (with ASCII symbols),
   using boost::spirit as a parser (so we'd have RPN) and boost::phoenix
   instead of boost::lambda (because the same author wrote spirit and
   phoenix).

If we're going to develop a custom solution, then let's make sure it's tuned
to lmi's needs, which don't really match what the scientific-computing
community has done. They ask "how can I make this as fast as FORTRAN?", but
I'd ask "how can I make this as clear as APL?". APL has a reputation as a
read-only language, but that's because it lacks good control constructs
(we can use C's), and because many programmers don't know how to use it--as
a mathematical notation, it's very clear.