Matrix computations on the C++11 templates.

matrix.cc

#include "matrix.h"

template class Matrix<4, 1>;
template class Matrix<4, 4>;

matrix.h

#ifndef CAD_MATRIX_H_
#define CAD_MATRIX_H_

// Empty type to identify return type of functions.
struct None {};

// templates prototypes to work with linked types.
template <class T> class MatrixDelegate;
template <unsigned int N,
          unsigned int M,
          class value_t = double,
          template <class> class Delegate = MatrixDelegate,
          class this_t = None> class Matrix;

// All work with entire matrix data are made by MatrixDelegate classes.
// You cannot construct an object of this class but you can redefine some
// methods to make it faster. Than you can use Matrix template argument.
template <class T>
class MatrixDelegate {
 public:
  // Entire matrix value type.
  typedef typename T::value_type value_t;
  // Matrix size.
  static const unsigned int N = T::HEIGHT;
  static const unsigned int M = T::WIDTH;
  // Addition.
  template <class other_value_t, template <class> class OtherDelegate>
  static void add(T& dest,
                  const Matrix<N, M, other_value_t, OtherDelegate>& arg) {
    other_value_t *from = arg.data();
    for (value_t *it = dest.data(), *end = it + N * M; it != end; ++it, ++from)
      *it += *from;
  }
  static void add(T& dest, const value_t& arg) {
    for (value_t *it = dest.data(), *end = it + N * M; it != end; ++it)
      *it += arg;
  }
  // Subtraction.
  template <class other_value_t, template <class> class OtherDelegate>
  static void sub(T& dest,
                  const Matrix<N, M, other_value_t, OtherDelegate>& arg) {
    other_value_t *from = arg.data();
    for (value_t *it = dest.data(), *end = it + N * M; it != end; ++it, ++from)
      *it -= *from;
  }
  static void sub(T& dest, const value_t& arg) {
    for (value_t *it = dest.data(), *end = it + N * M; it != end; ++it)
      *it -= arg;
  }
  // Multiplication
  template <unsigned int P,
            class other_value_t,
            template <class> class OtherDelegate>
  static Matrix<N, P, value_t, MatrixDelegate> mul(const T& left,
      const Matrix<M, P, other_value_t, OtherDelegate>& right) {
    Matrix<N, P, value_t, MatrixDelegate> res;
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < P; ++j) {
        value_t v = 0.0f;
        for (int k = 0; k < M; ++k)
          v += left.get(i, k) * right.get(k, j);
        res.set(i, j, v);
      }
    return res;
  }
  static void mul(T& dest, const value_t& arg) {
    for (value_t *it = dest.data(), *end = it + N * M; it != end; ++it)
      *it *= arg;
  }
  // Division by number.
  static void div(T& dest, const value_t& arg) {
    for (value_t *it = dest.data(), *end = it + N * M; it != end; ++it)
      *it /= arg;
  }
  // Dot product.
  template <class other_value_t, template <class> class OtherDelegate>
  static value_t dot(const T& left,
                     const Matrix<N, M, other_value_t, OtherDelegate>& right) {
    value_t res(0.0f);
    other_value_t *rit = right.data();
    for (value_t *it = left.data(), *end = it + N * M; it != end; ++it, ++rit)
      res += *it * *rit;
  }

 private:
  // Constructor is private to prevent object creation.
  MatrixDelegate();
};

// MatrixType structure used to detect the return value type of operations.
template <unsigned int N,
          unsigned int M,
          class value_t,
          template <class> class Delegate,
          class this_t>
struct MatrixType {
  typedef this_t this_type;
};
// Specialization for None as this_t.
template <unsigned int N,
          unsigned int M,
          class value_t,
          template <class> class Delegate>
struct MatrixType<N, M, value_t, Delegate, None> {
  typedef Matrix<N, M, value_t, Delegate, None> this_type;
};

// Matrix class. You can inherit it in your class but must provide correct
// this_t argument to make sure that matrix operations works with your type.
template <unsigned int N,
          unsigned int M,
          class value_t,
          template <class> class Delegate,
          class this_t>
class Matrix {
 public:
  // Entire value_type. Defined for the delegate class.
  typedef value_t value_type;
  // Return type. Matrix<N, M, value_t, Delegate, None> if this_t is None.
  // Else - this_t.
  typedef typename MatrixType<N, M, value_t, Delegate, this_t>::this_type
      this_type;
  // Delegate instance.
  typedef Delegate<Matrix> delegate_type;
  // Size constants for delegate.
  static const unsigned int HEIGHT = N;
  static const unsigned int WIDTH = M;

  // Default constructor.
  Matrix()
    : data_(new value_t[N * M]) {
    clear();
  }

  // Copy constructor.
  Matrix(const this_type& arg)
    : data_(new value_t[N * M]) {
    for (value_t *it = data_, *end = it + N * M, *from = arg.data_; it != end;
        ++it, ++from)
      *it = *from;
  }

  // Convertion constructor.
  template <class other_value_t, template <class> class OtherDelegate>
  explicit Matrix(const Matrix<N, M, other_value_t, OtherDelegate>& arg)
    : data_(new value_t[N * M]) {
    other_value_t *from = arg.data_;
    for (value_t *it = data_, *end = data_ + N * M; it != end; ++it, ++from)
      *it = *from;
  }
  virtual ~Matrix() {};

  // Addition.
  template <class T>
  this_type& operator+= (const T& arg) {
    Delegate<this_type>::add(*this, arg);
    return *this;
  }
  template <class T>
  this_type operator+ (const T& arg) const {
    return this_type(*this) += arg;
  }
  // Subtraction.
  template <class T>
  this_type& operator-= (const T& arg) {
    Delegate<this_type>::sub(*this, arg);
    return *this;
  }
  template <class T>
  this_type operator- (const T& arg) const {
    return this_type(*this) -= arg;
  }
  // Multiplication.
  template <class T>
  this_type& operator*= (const T& arg) {
    Delegate<this_type>::mul(*this, arg);
    return *this;
  }
  template <class T>
  this_type operator* (const T& arg) const {
    return this_type(*this) *= arg;
  }
  // Division by number.
  this_type& operator/= (const value_t& arg) {
    Delegate<this_type>::div(*this, arg);
    return *this;
  }
  this_type operator/ (const value_t& arg) const {
    return this_type(*this) /= arg;
  }
  // Dot product.
  template <class T>
  this_type& operator%= (const T& arg) {
    Delegate<this_type>::dot(*this, arg);
    return *this;
  }
  template <class T>
  this_type operator% (const T& arg) const {
    return this_type(*this) %= arg;
  }
  
  // Assignment operator.
  template <class... Args>
  this_type& operator= (const Matrix<N, M, Args...>& arg) {
    this_type(arg).Swap(*this);
    return *this;
  }
  // Swap entire contents of matrixes.
  void Swap(this_type& other) {
    value_t *ptr = data_;
    data_ = other.data_;
    other.data_ = ptr;
  }

  // Returns value of cell [i+1, j+1]
  value_t get(unsigned int i, unsigned int j) const {
    if (i < N && j < M)
      return data_[i * M + j];
    return 0.0f;
  }

  // Sets value of cell [i+1, j+1]
  void set(unsigned int i, unsigned int j, value_t v) {
    if (i < N && j < M)
      data_[i * M + j] = v;
  }

  // Sets contents of matrix to 0.0f.
  void clear() {
    for (value_t *it = data_, *end = data_ + N * M; it != end; ++it)
      *it = 0.0f;
  }

  // Returns pointer to the entire data.
  value_t * data() { return data_; }
  const value_t * data() const { return data_; }

 private:
  // Entire data array. It made as pointer to make possible fast swap operation.
  value_t *data_;
};

// Make compilation faster.
// Instance this template specialisations in the matrix.cc
extern template class Matrix<4, 1>;
extern template class Matrix<4, 4>;

#endif  // CAD_MATRIX_H_