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