Expression templates

[Don Clugston]

Christian Kamm wrote:
>> I'm making some headway on an nd-array class that uses BLAS/LAPACK 
>> under the hood. It would be nifty if I could somehow get this:
>>
>>       R = A * x + y
>>
>> to be groked and translated to a single call to of BLAS's 
>> saxpy/daxpy/caxpy/zaxpy routines.
> 
> This would be fabulous! (and shouldn't it work similarly to the way 
> Don's code detects constant folding? If there's a matrix-vector product 
> and one uses opAdd with another vector, modify the call to xGEMV 
> accordingly)
 >
> I would definitely use something that makes accessing the BLAS routines 
> as easy as that. It should also be possible to automatically use the 
> calls for symmetric, hermitian, etc. matrices, if the matrix type 
> reflects these properties...

That's all feasible. It's not yet feasible to turn
  A = t * B + A;
into a DAXPY call, because there's no way of knowing that the two 'A's 
are the same. But
A += t * B;
could be identified as DAXPY, without any trouble. Ditto for the 
matrix/vector operations.
(Now if my proposed enhancement to operator overloading gets implemented 
(which seems likely to happen eventually), the compiler would turn
A = t*B + A;
into A.opAddAssign_r(t*B);
and then it could be identified as DAXPY.
)

>> The theory with BLAS anyway, is that you'll have access to a 
>> hand-tuned version for the platform you're on, so it'll be hard to 
>> beat with any kind of register-level compile-time micromanagement.
> 
> I agree. Writing near-optimal assembler code might be doable for vector 
> operations (and a pretty cool thing too), but for matrix-vector and 
> matrix-matrix operations, I'd go with the functions from BLAS level 2 
> and 3.

Absolutely agree. Cache effects are dominant in BLAS 2 and 3, and a lot 
of very clever people have put a lot of work put into getting that right.
It's only in BLAS1 that there's much chance of benefit from on-the-fly 
code creation.

> Please keep us updated on your progress or make it a project on dsource 
> - I'd want to check it out and contribute, if time allows.

It's been a bit of a diversion for me. I'll put the code into my 
'mathextra' project on dsource (name: mathextra.Blade), but won't do 
anything with it till after  D 1.0. (I've done a lot of proof-of-concept 
things to try to find the weak points of the D template system; but it 
generally doesn't make sense to put much effort into them until we get 
the first feature freeze; the language has improved so much over the 
past few months, this stuff tends to get obselete very quickly).

-Don.