[boost] [rfc] A light weight statically typed matrix template class (again?)

From: Andreas Harnack (ah.boost.04_at_[hidden])
Date: 2007-12-06 13:37:51

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