Boost logo

Boost :

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

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

would be superdoubleplusgood.

Incidentally you may be interested in the Tensor Contraction Engine,

or the tensor template library

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

Boost list run by bdawes at, gregod at, cpdaniel at, john at