[rfc] A light weight statically typed matrix template class (again?)
Hi @, shouldn't there be at least an elementary matrix class available in the standard? I notice there have been some discussions in the past, but there doesn't seem to be anything in the pipeline. So, is there any interest in a (really!) light weight statically typed matrix template class? I'm thinking of something like: #include <algorithm> #include <functional> template<typename T, unsigned int M, unsigned int N> class matrix { protected: T m[M][N]; public: typedef T value_type; enum { rows = M, columns = N }; T* data() { return &**m; } T const* data() const { return &**m; } matrix() {} matrix(T const& c) { std::fill(data(), data()+M*N, c); } matrix(matrix const& m) { std::copy(m.data(), m.data()+M*N, data()); } template<typename U> matrix(matrix<U,M,N> const& m) { std::copy(m.data(), m.data()+M*N, data()); } template<typename U> matrix& operator=(matrix<U,M,N> const& m) { std::copy(m.data(),m.data()+M*N,data()); return *this; } template<class UnaryFunction> matrix& map(UnaryFunction); template<class BinaryFunction> matrix& join(matrix const&, BinaryFunction); matrix& operator+=(matrix const& m) { return join(m, std::plus<T>() ); } matrix& operator-=(matrix const& m) { return join(m, std::minus<T>() ); } matrix& operator*=(T const& c) { return map(std::bind2nd(std::multiplies<T>(), c)); } matrix& operator/=(T const& c) { return map(std::bind2nd(std::divides<T>(), c)); } matrix& operator%=(T const& c) { return map(std::bind2nd(std::modulus<T>(), c)); } T& operator()(unsigned int i, unsigned int j) { return m[i][j]; } T const& operator()(unsigned int i, unsigned int j) const { return m[i][j]; } }; template<typename T, unsigned int M, unsigned int N> template<class UnaryFunction> matrix<T,M,N>& matrix<T,M,N>::map( UnaryFunction op) { std::transform(data(), data()+M*N, data(), op); return *this; } template<typename T, unsigned int M, unsigned int N> template<class BinaryFunction> matrix<T,M,N>& matrix<T,M,N>::join(matrix const& m, BinaryFunction op) { std::transform(data(), data()+M*N, m.data(), data(), op); return *this; } // global converter; allows type deduction template<typename T, unsigned int M, unsigned int N> inline matrix<T,M,N> mtrx(T (&a)[M][N]) { matrix<T,M,N> tmp; std::copy(&**a, &**a+M*N, tmp.data()); return tmp; } // unary and binary operators template<typename T, unsigned int M, unsigned int N> inline bool operator==(matrix<T,M,N> const& m1, matrix<T,M,N> const& m2) { return std::equal( m1.data(), m1.data()+M*N, m2.data()); } template<typename T, unsigned int M, unsigned int N> inline bool operator!=(matrix<T,M,N> const& m1, matrix<T,M,N> const& m2) { return ! std::equal( m1.data(), m1.data()+M*N, m2.data()); } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator-(matrix<T,M,N> const& m) { matrix<T,M,N> tmp = T(); return tmp-=m; } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator+(matrix<T,M,N> tmp, matrix<T,M,N> const& m) { return tmp+=m; } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator-(matrix<T,M,N> tmp, matrix<T,M,N> const& m) { return tmp-=m; } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator*(matrix<T,M,N> tmp, T const& c) { return tmp*=c; } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator*(T const& c, matrix<T,M,N> tmp) { return tmp*=c; } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator/(matrix<T,M,N> tmp, T const& c) { return tmp/=c; } template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator%(matrix<T,M,N> tmp, T const& c) { return tmp%=c; } The advantage of this approach is its simplicity, in fact, it's so straight forward I'm surprised it's not already in the standard. There's nothing in there (yet), that's not already in the STL. It's just like bundling a few things together and making them available in a more convenient notation and/or from a different point of view. There are still a few essential but fairly basic things missing, like matrix multiplication and transposing a matrix. Elementary row operations might be useful as well for those, who want to provide complex algorithms, but that's it, basically. There shouldn't be anything in there that's more complicated than that. (Matrices are such a general concept that it can be very tricky, if not impossible, to provide any guaranties for anything above the basic stuff.) Such a matrix class would be useful in all situations, where matrix notation can lead to shorter and better readable code, while matrices remain reasonably small and/or performance is not a primary concern. I'm thinking of domains like graphical programming. Matrices could be passed around by value and almost be treated as build in types. They could also be used as building block for more advanced libraries. I came across the problem of a missing matrix class when refining my Euclid class. There are some application problems, that are best solved using matrices (like vector space transformation). The current release therefore contains a (diagonal) matrix class, but I feel it doesn't really belong there. Matrices are by far general enough, to live in there own library. The available Linear Algebra libraries, however, are by far too complex for such simple tasks. I've uploaded a first version in the vault (simple_matrix.zip), to have something to start with. It already contains some missing features in a very naive implementation and comes with a test suit and an example program. Any interest? Comments most welcome! Andreas
On 12/06/07 12:37, Andreas Harnack wrote:
Hi @,
shouldn't there be at least an elementary matrix class available in the standard? I notice there have been some discussions in [snip] Any interest? Comments most welcome! Andreas What's the advantage of this over the specialization of multi_array:
http://www.boost.org/libs/multi_array/doc/reference.html to just 2 dimensions?
Larry Evans schrieb:
What's the advantage of this over the specialization of multi_array:
http://www.boost.org/libs/multi_array/doc/reference.html
to just 2 dimensions?
I might have missed something, but I couldn't find any matrix operations in the multi_array class?! Further, a MultiArray 'defines an interface to hierarchically nested containers'. A matrix is neither really nested nor am I too happy with the way of viewing a matrix primarily as a container. Technically, of course, it is, but for the applications I have in mind it's more important what you can do with a matrix as a whole (like adding, multiplying with a scalar etc.) rather than their internal structure. That's one reason there are no iterators defined yet. They will be important at some point for people who want to implement algorithms over matrices (or operators for combining matrices with other types like vectors), but for now I think the focus should be on operations over matrices as a whole. Andreas
You might be interested in looking at blitz++ or tvmet. They both use expression templates to statically unroll various operators (for small sized vectors/matricies anyway). If I ever find a spare minute, I've always wanted to boostify one of those libs. Boost also has a UBLAS impl. which has a matrix and vector class (although not static). Cheers, Chris On Dec 6, 2007 1:37 PM, Andreas Harnack <ah.boost.04@justmail.de> wrote:
Hi @,
shouldn't there be at least an elementary matrix class available in the standard? I notice there have been some discussions in the past, but there doesn't seem to be anything in the pipeline. So, is there any interest in a (really!) light weight statically typed matrix template class? I'm thinking of something like:
#include <algorithm> #include <functional>
template<typename T, unsigned int M, unsigned int N> class matrix { protected: T m[M][N]; public: typedef T value_type; enum { rows = M, columns = N }; T* data() { return &**m; } T const* data() const { return &**m; }
matrix() {} matrix(T const& c) { std::fill(data(), data()+M*N, c); } matrix(matrix const& m) { std::copy(m.data(), m.data()+M*N, data()); }
template<typename U> matrix(matrix<U,M,N> const& m) { std::copy(m.data(), m.data()+M*N, data()); }
template<typename U> matrix& operator=(matrix<U,M,N> const& m) { std::copy(m.data(),m.data()+M*N,data()); return *this; }
template<class UnaryFunction> matrix& map(UnaryFunction); template<class BinaryFunction> matrix& join(matrix const&, BinaryFunction);
matrix& operator+=(matrix const& m) { return join(m, std::plus<T>() ); } matrix& operator-=(matrix const& m) { return join(m, std::minus<T>() ); } matrix& operator*=(T const& c) { return map(std::bind2nd(std::multiplies<T>(), c)); } matrix& operator/=(T const& c) { return map(std::bind2nd(std::divides<T>(), c)); } matrix& operator%=(T const& c) { return map(std::bind2nd(std::modulus<T>(), c)); }
T& operator()(unsigned int i, unsigned int j) { return m[i][j]; } T const& operator()(unsigned int i, unsigned int j) const { return m[i][j]; } };
template<typename T, unsigned int M, unsigned int N> template<class UnaryFunction> matrix<T,M,N>& matrix<T,M,N>::map( UnaryFunction op) { std::transform(data(), data()+M*N, data(), op); return *this; }
template<typename T, unsigned int M, unsigned int N> template<class BinaryFunction> matrix<T,M,N>& matrix<T,M,N>::join(matrix const& m, BinaryFunction op) { std::transform(data(), data()+M*N, m.data(), data(), op); return *this; }
// global converter; allows type deduction template<typename T, unsigned int M, unsigned int N> inline matrix<T,M,N> mtrx(T (&a)[M][N]) { matrix<T,M,N> tmp; std::copy(&**a, &**a+M*N, tmp.data()); return tmp; }
// unary and binary operators template<typename T, unsigned int M, unsigned int N> inline bool operator==(matrix<T,M,N> const& m1, matrix<T,M,N> const& m2) { return std::equal( m1.data(), m1.data()+M*N, m2.data()); }
template<typename T, unsigned int M, unsigned int N> inline bool operator!=(matrix<T,M,N> const& m1, matrix<T,M,N> const& m2) { return ! std::equal( m1.data(), m1.data()+M*N, m2.data()); }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator-(matrix<T,M,N> const& m) { matrix<T,M,N> tmp = T(); return tmp-=m; }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator+(matrix<T,M,N> tmp, matrix<T,M,N> const& m) { return tmp+=m; }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator-(matrix<T,M,N> tmp, matrix<T,M,N> const& m) { return tmp-=m; }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator*(matrix<T,M,N> tmp, T const& c) { return tmp*=c; }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator*(T const& c, matrix<T,M,N> tmp) { return tmp*=c; }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator/(matrix<T,M,N> tmp, T const& c) { return tmp/=c; }
template<typename T, unsigned int M, unsigned int N> matrix<T,M,N> inline operator%(matrix<T,M,N> tmp, T const& c) { return tmp%=c; }
The advantage of this approach is its simplicity, in fact, it's so straight forward I'm surprised it's not already in the standard. There's nothing in there (yet), that's not already in the STL. It's just like bundling a few things together and making them available in a more convenient notation and/or from a different point of view.
There are still a few essential but fairly basic things missing, like matrix multiplication and transposing a matrix. Elementary row operations might be useful as well for those, who want to provide complex algorithms, but that's it, basically. There shouldn't be anything in there that's more complicated than that. (Matrices are such a general concept that it can be very tricky, if not impossible, to provide any guaranties for anything above the basic stuff.)
Such a matrix class would be useful in all situations, where matrix notation can lead to shorter and better readable code, while matrices remain reasonably small and/or performance is not a primary concern. I'm thinking of domains like graphical programming. Matrices could be passed around by value and almost be treated as build in types. They could also be used as building block for more advanced libraries.
I came across the problem of a missing matrix class when refining my Euclid class. There are some application problems, that are best solved using matrices (like vector space transformation). The current release therefore contains a (diagonal) matrix class, but I feel it doesn't really belong there. Matrices are by far general enough, to live in there own library. The available Linear Algebra libraries, however, are by far too complex for such simple tasks.
I've uploaded a first version in the vault (simple_matrix.zip), to have something to start with. It already contains some missing features in a very naive implementation and comes with a test suit and an example program.
Any interest? Comments most welcome! Andreas
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
Chris Fairles schrieb:
You might be interested in looking at blitz++ or tvmet. They both use expression templates to statically unroll various operators (for small sized vectors/matricies anyway). If I ever find a spare minute, I've always wanted to boostify one of those libs.
Boost also has a UBLAS impl. which has a matrix and vector class (although not static).
Hi Chris, I checked all of them and I found them way to complicated for what I have in mind. Let's assume we have an application with a graphical user interface, an object described by coordinates in three-dimensional space and we need to project that object onto a two-dimensional array of pixels; mathematically, that's a simple 3x2 matrix multiplication, most likely in integer arithmetic. Personally, the last thing I want to do is to include a full scale Algebraic Library into the user face just for such a simple task. What I want in such a case is a matrix class, that is as easy to use and as widely available as, let's say, the a complex numbers class. By the way, the test suite actually uses uBLAS to verify the results of the matrix operations. Not exactly a fool prove way, but for the time being I think it's OK. Andreas
Hi, On Thu, Dec 06, 2007 at 07:37:51PM +0100, Andreas Harnack wrote:
shouldn't there be at least an elementary matrix class available in the standard? I notice there have been some discussions in the past, but there doesn't seem to be anything in the pipeline. So, is there any interest in a (really!) light weight statically typed matrix template class? I'm thinking of something like:
template<typename T, unsigned int M, unsigned int N> class matrix
a few days ago I started exactly as you described (except that I used old style for loops instead of STL stuff :-) a TinyDenseMatrix class to implement the Moore Penrose pseudo inverse for real tiny matrices (4 x n, n<20). A very big disadvantage is using fixed dimensions. There is no need for memory allocation in this case but it restricts the usage dramatically! Nearly always the dimensions are variables, there are only a few rare cases where you know these before. So how about making N and M variables and providing one or two good memory allocators?
T* data() { return &**m; }
Why not &m[0][0]?
The advantage of this approach is its simplicity, in fact, it's
Right, I wrote such stuff already multiple times even if I care normally about sparse matrices with dimensions > 1 000 000. Jens
Jens Seidel schrieb:
Hi,
a few days ago I started exactly as you described (except that I used old style for loops instead of STL stuff :-) a TinyDenseMatrix class to implement the Moore Penrose pseudo inverse for real tiny matrices (4 x n, n<20).
A very big disadvantage is using fixed dimensions. There is no need for memory allocation in this case but it restricts the usage dramatically!
Nearly always the dimensions are variables, there are only a few rare cases where you know these before. So how about making N and M variables and providing one or two good memory allocators?
T* data() { return &**m; }
Why not &m[0][0]?
Just a matter of taste, I guess. Actually I'm not even sure yet that it's wise to use a two dimensional array at all. T m[M*N] would work equally well and since there's no interface yet that relies on that structure, the latter might even be the simpler approach. The two-dimensional array just felt to be more 'natural' for a matrix. (whatever that means :-) ) But thanks for your support, Andreas
The advantage of this approach is its simplicity, in fact, it's
Right, I wrote such stuff already multiple times even if I care normally about sparse matrices with dimensions > 1 000 000.
Jens _______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
On 12/06/07 16:07, Andreas Harnack wrote: [snip]
Jens Seidel schrieb: [snip]
Nearly always the dimensions are variables, there are only a few rare cases where you know these before. So how about making N and M variables and providing one or two good memory allocators?
T* data() { return &**m; }
Why not &m[0][0]?
Just a matter of taste, I guess. Actually I'm not even sure yet that it's wise to use a two dimensional array at all. T m[M*N] would work equally well and since there's no interface yet that relies on that structure, the latter might even be the simpler approach. The two-dimensional array just felt to be more 'natural' for a matrix. (whatever that means :-) ) I think n-dimension would be pretty easy. Instead of:
template<typename T, unsigned int M, unsigned int N> class matrix { ... }; you could provide: template<typename Value, typename Shape> class array_fix_shape { ... }; where, using concepts maybe, you could constrain Shape to be some instance of mpl::vector_c<unsigned,S0,S1,...,Sn>. Then, the rank of array_fix_size would be Shape::size. You could calculate of strides using some mpl metafunction on Shape: typedef some-mpl-function < Shape >::type strides ; and then allocate storage as: Value values [ mpl::at_c<strides,Strides::size>::value ] ; Then operator[] could be defined as: Value& operator[](boost::array<unsigned,Shape::size>const& indices) { return values[some_offset_function_<strides>(indices)]; } where some_offset_function simply performs a dot product of the compile time vector, strides, with the runtime vector, indices.
On 12/06/07 18:40, Larry Evans wrote: [snip]
Then, the rank of array_fix_size would be Shape::size. You could calculate of strides using some mpl metafunction on Shape:
typedef some-mpl-function < Shape >::type strides ;
and then allocate storage as:
Value values [ mpl::at_c<strides,Strides::size>::value ] ;
Then operator[] could be defined as:
Value& operator[](boost::array<unsigned,Shape::size>const& indices) { return values[some_offset_function_<strides>(indices)]; }
where some_offset_function simply performs a dot product of the compile time vector, strides, with the runtime vector, indices. Andreas,
There actually some code to illustrate this in the fold_seq_test.zip file at: http://www.boost-consulting.com/vault/index.php?&directory=Template%20Metaprogramming It's in the test_array_index_offset function. The fold_seq.hpp file includes several other files on my local disk and I could upload those if you're interested, but I think the test_array_index_offset code illustrates the general idea. Hm... No actually it differs from what's suggested in my previous post. Instead of that, there would be something like: typedef unsigned index_type ; struct inner_product_static_dynamic_ftor { explicit inner_product_static_dynamic_ftor ( index_type const* indices ) : cur_index_ptr(indices) , offset(0) {} index_type const* cur_index_ptr ; index_type offset ; template<class PartialProductElement> void operator()(PartialProductElement)const { offset+=PartialProductElement::value*(*cur_index_ptr); ++cur_index_ptr; } }; template < class Strides
index_type offset_at_indices ( array<index_type,mpl::size<Strides>::value>const& a ) { index_type const*indices=a.data(); inner_product_static_dynamic_ftor ip(indices); typedef mpl::pop_back<Strides>::type strides;//rm total array size. mpl::for_each<strides>(ip); return ip.offset; } *WARNING*: not compiled.
On 12/07/07 11:05, Larry Evans wrote: [snip]
There actually some code to illustrate this in the fold_seq_test.zip file at:
http://www.boost-consulting.com/vault/index.php?&directory=Template%20Metaprogramming
It's in the test_array_index_offset function. The fold_seq.hpp file includes several other files on my local disk and I could upload those if you're interested, but I think the test_array_index_offset code illustrates the general idea. Hm... No actually it differs from what's suggested in my previous post. Instead of that, there would be something like: fold_seq_test.zip has been replaced with:
array_dyn_indices_static_shape_offset_demo.zip which requires nothing on my local disk. It will compile with just mpl and it's much simpler and is essentially what was suggested in my previous post (i.e. there's something like inner_product_static_dynamic_ftor ).
On 12/07/07 19:14, Larry Evans wrote:
On 12/07/07 11:05, Larry Evans wrote: [snip]
There actually some code to illustrate this in the fold_seq_test.zip file at:
http://www.boost-consulting.com/vault/index.php?&directory=Template%20Metaprogramming
[snip] fold_seq_test.zip has been replaced with:
array_dyn_indices_static_shape_offset_demo.zip
which requires nothing on my local disk. It will [snip] In same vault directory, there's now:
array_multi_dyn_indices_static_shape_demo.cpp which doesn't bother with offset calcs and array<int,size> indices. It just uses boost::array and a new template: array_multi < ValueType , mpl::vector_c<extent_type,size0,size1,...,sizeN>
for any N to implement any rank array. It inherits from boost::array; hence, has all the it's methods. The .cpp file compiles and runs. That's all the testing that's been done.
Larry Evans schrieb:
On 12/07/07 19:14, Larry Evans wrote:
On 12/07/07 11:05, Larry Evans wrote: [snip]
Hi Larry, I'm sorry it took some time for me to answer, but I'm sure you know what a deadline is. :-) Anyway, thanks a lot for all your efforts; I'm afraid, however, I don't quite know what you're aiming at. As I understand you want to extend the concept of the matrix as a two dimensional array to that of a multi dimensional array, assuming that they follow the same principles? That sounds interesting, though I'm not sure that's sensible. I'm aiming at simplicity. A matrix is a very simple mathematical concept that can be useful in many situations, yet there's no data type supporting it adequately (whatever that means) in C++; All that's available (well, actually it's quite a lot) are fully-fledged libraries that apparently all are aiming to high level scientific computing. Tensors seem even more advanced to me. I do plan indeed to add a more template parameters, first and foremost to allow to specify special matrix forms, like diagonal matrices, which have some special properties, that might be interesting to the user. Different storage models might interesting too. I'm not yet sure how deal with variable sized matrices (or rather whether to deal with them at all). They certainly will be important for some people, but would come at the price to give up compile time size checking. Andreas
On 12/10/07 15:31, Andreas Harnack wrote:
Larry Evans schrieb:
On 12/07/07 19:14, Larry Evans wrote:
On 12/07/07 11:05, Larry Evans wrote: [snip] Hi, Andreas.
[snip]
efforts; I'm afraid, however, I don't quite know what you're aiming at. As I understand you want to extend the concept of the matrix as a two dimensional array to that of a multi dimensional array, assuming that they follow the same principles?
Yes.
That sounds interesting, though I'm not sure that's sensible. I'm aiming at simplicity.
I understand, and I agree that this: template < class Value , class Shape
struct array_multi : public mpl::eval_if_c < bool(mpl::size<Shape>::type::value-1)//if rank>1 , mpl::identity //rank>1 < array < array_multi < Value , typename mpl::pop_front < Shape >::type > , mpl::front < Shape >::type::value > > , mpl::identity //rank==1 < array < Value , mpl::front < Shape >::type::value > >
::type { };
is certainly more complex that this: template < typename T , unsigned int M , unsigned int N
class matrix { ... }; however, the array_multi "reuses" array; hence, is less complex *and* more general since it can handle any rank array. Hence, I thought the extra complexity of, let's call it the inheritance spec, would be worth the extra generality. After all, if just simplicity is what you want, then just use boost::array. OTOH, if your particular needs are just for rank 2 arrays, then your matrix is "made to order". OTOH, as soon as someone comes up with a use-case for rank>2, you have to reinvent what's essentially already provided by array_multi. I guess it just depends on whether one thinks the "bit" extra complexity of array_multi is worth the increased generality (I know, that's obvious, but that seems to be the essential reason for the matrix vs. array_multi argument). [snip]
I'm not yet sure how deal with variable sized matrices (or rather whether to deal with them at all). They certainly will be important for some people, but would come at the price to give up compile time size checking.
Doesn't multi_array implement variable size arrays. I see there's: template <typename ValueType, std::size_t NumDims, typename Allocator = std::allocator<ValueType> > class multi_array { ... template <typename ExtentList> explicit multi_array(const ExtentList& sizes, const storage_order_type& store = c_storage_order(), const Allocator& alloc = Allocator()); ... }; on page: http://www.boost.org/libs/multi_array/doc/reference.html#multi_array so I guess it just fixes the rank (the NumDims template parameter). And, BTW, apparently the multi_array authors seem to think that the cost of extra template complexity is worth the gain in generality; otherwise, they wouldn't have bothered with their library. -regards, Larry
Andreas Harnack wrote:
Hi @,
shouldn't there be at least an elementary matrix class available in the standard? I notice there have been some discussions in the past, but there doesn't seem to be anything in the pipeline.
You're likely to find that such a thing would be a huge disappointment (cf valarray :-) Doing numerics "right" is REALLY hard, particularly when it comes to matrix operations. Anything that manages to pass the hurdles of standardization is likely to be: 1) underspecified to the point of worthlessness. This is not a shot at the Committee ... it's just a really hard problem to do right, 2) missing key functionality for most subsets of the C++ user community. Unfortunately, it won't be the SAME missing functionality for different those different subgroups :-) For instance, you may only need the ability to multiply small matrices together, but I need to do Gauss-Jordan elimination, my neighbor only needs to find eigenvalues, and his cousin only needs to find matrix permutations. Where do you draw the line on what is "simple" enough to go into your "simple" matrix class? 3) so slow and numerically unstable as to be worthless for most serious uses. Naive implementations are really BAD news on these fronts. This is an enormous QOI issue, and reasonable implementations need to be implemented by numerics experts, not your garden variety programmer (I know just enough to know how much I don't know). To see just how difficult a problem this is, you might want to survey the attempts to implement matrix functionality in C++: ublas, Blitz, MTL, POOMA, and a host of others. All of them are enormous efforts, they all make different trade offs, all of them are missing key functionality, and all of them provide different performance and correctness guarantees. -- ------------------------------------------------------------------------------- Kevin Lynch voice: (617) 353-6025 Physics Department Fax: (617) 353-9393 Boston University office: PRB-361 590 Commonwealth Ave. e-mail: krlynch@bu.edu Boston, MA 02215 USA http://budoe.bu.edu/~krlynch -------------------------------------------------------------------------------
participants (5)
-
Andreas Harnack -
Chris Fairles -
Jens Seidel -
Kevin Lynch -
Larry Evans