Re: [boost] different matrix library?

Subject: Re: [boost] different matrix library?
From: Edward Grace (ej.grace_at_[hidden])
Date: 2009-08-11 10:16:14

>> In principle yes - however the efficient implementation of the above
>> is likely to be a lot tricker than it might seem.
> will you as a user concern about how tricky an implementation is?
> or will you rather care about how convinient and clear public  
> interface
> is?

The latter - naturally! If you can make it extensible so that the  
back end could (for example) be thread, CUDA, and MPI aware so much  
the better. ;-)

>> The 'advantage' of BLAS and other routines as I see it is many many
>> man-years of optimisation and tweaking on various architectures as
>> well as good generic implementations (ATLAS).  I'd venture that it's
>> still  going to be hard to beat!
> the advantage of BLAS as far as i can see is if you take original
> fortran implementation and use it as it is (wow, how many 'as's and
> 'is's)
> as i understand that's the implementation you are talking about when
> mentioning 'man-years'
> i hope to make an implementation 'as good as' but which will exploit
> all of c++ advantages

I look forward to it!  Competing with top-end platform tuned  
FORTRAN90 compilers is not a challenge for the feint hearted!

>> I am slightly unsure, is your proposal a rewrite of the linear
>> algebra routines or a wrapper that conceptually maps calls in the
>> following manner?
>>   yourlib::operator*(foo,bar) -> ublas::prod(foo,bar)
> definitely it's not a wrapper
> so i guess it's a rewrite in c++ style

Well.  You're a braver man than I!  Building a templated linear  
algebra library is a lot of work!

>> Are you aware of Blitz++ and POOMA?
>> Blitz offers a very (for the mathematical physicist) intuitive  
>> tensor-
>> like approach, an example:
>>      // This expression will set
>>      //
>>      // c    = a   * b
>>      //  ijk    ik    kj
>>      C = A(i,k) * B(k,j);
> i'm aware of blitz++ (i peeped at the implementation a little when i
> wrote my own lib)
> i don't use it because i don't like the concept

How so?  I don't use Blitz++ as it appears to have stagnated - very  
unfortunate. The concept seems quite attractive however.  For example  
Blitz++ style expressiveness need not constrain you to scalars  
(tensors of order zero), vectors (tensors of order one) or matricies  
(tensors of order 2) and can allow far more flexibility than simple  
matrix calculus.

For example, viewing _1 as placeholders in the a similar manner to  
boost::bind, being able to write things (I'm ignoring contra- vs. co- 
variance) like:

   typedef std::complex complex;
   boost::yourlib::tensor<complex,0> epsilon; // Same as std::complex  
epsilon since it's a scalar;
   boost::yourlib::tensor<complex,1> P,E;
   boost::yourlib::tensor<complex,2> Chi1; // Other template  
arguments could indicate upper or lower indices.
   boost::yourlib::tensor<complex,3> Chi2;
   boost::yourlib::tensor<complex,4> Chi3;

   //... stuff
   using boost::yourlib::einstein_placeholders;

   P(_1) = Chi1(_1,_2)*E(_2) + Chi2(_1,_2,_3)*E(_2)*E(_3) + Chi3 
(_1,_2,_3,_4)*E(_2)*E(_3)*E(_4);

so that the Einstein summation convention was observed (summation  
over repeated indices),

   http://en.wikipedia.org/wiki/Einstein_notation
   http://en.wikipedia.org/wiki/Tensor

would be superdoubleplusgood.

Incidentally you may be interested in the Tensor Contraction Engine,

   http://www.csc.lsu.edu/~gb/TCE/

or the tensor template library

   http://ttl.yanenko.com/

-ed
------------------------------------------------
"No more boom and bust." -- Dr. J. G. Brown, 1997