The Artima Developer Community
Sponsored Link

Weblogs Forum
Storing and Tracking Matrix Rows and Columns at Compile-Time

1 reply on 1 page. Most recent reply: May 14, 2005 9:44 PM by Christopher Diggins

Welcome Guest
  Sign In

Go back to the topic listing  Back to Topic List Click to reply to this topic  Reply to this Topic Click to search messages in this forum  Search Forum Click for a threaded view of the topic  Threaded View   
Previous Topic   Next Topic
Flat View: This topic has 1 reply on 1 page
Christopher Diggins

Posts: 1215
Nickname: cdiggins
Registered: Feb, 2004

Storing and Tracking Matrix Rows and Columns at Compile-Time (View in Weblogs)
Posted: May 14, 2005 9:44 PM
Reply to this message Reply
Summary
Here is a method for for representing and computing rows and columns at compile-time, and to do matrix multiplication accordingly. (trickier than it might sound)
Advertisement
Disclaimer: I had horrible food poisoning recently so I can't be held responsible for the quality of this work ... good or bad ;-)

#include <iostream>
#include <numeric>
#include <valarray>
#include <cassert>

template<class Value_T, unsigned int Size_N, unsigned int Stride_N>
class Slice
{
public:
  // constructor
  Slice(Value_T* x) : p(x) { }
  // typedef
  typedef Value_T value_type;
  // only a forward iterator is provided,
  // others are left as an exercise for the reader
  struct Slice_forward_iter {
    Slice_forward_iter(Value_T* x) : p(x) { }
    Value_T& operator*() { return *p; }
    const Value_T& operator*() const { return *p; }
    Value_T& operator++() { Value_T* tmp = p; p += Stride_N; return *tmp; }
    Value_T& operator++(int) { return *(p+=Stride_N); }
    bool operator==(const Slice_forward_iter& x) { return p == x.p; }
    bool operator!=(const Slice_forward_iter& x) { return p != x.p; }
    Value_T* p;
  };
  typedef Slice_forward_iter iterator;
  typedef const Slice_forward_iter const_iterator;
  // public functions
  iterator begin() { return iterator(p); }
  iterator end() { return iterator(p + (Size_N * Stride_N)); }
  value_type& operator[](size_t n) { return *(p + (n * Stride_N)); }
  const value_type& operator[](size_t n) const { return *(p + (n * Stride_N)); }
  static size_t size() { return Size_N; }
  static size_t stride() { return Stride_N; }
private:
  // prevent default construction
  Slice() { };
  Value_T* p;
};

template<class Value_T, unsigned int Rows_N, unsigned int Cols_N>
class Matrix
{
public:
  typedef Slice<Value_T, Cols_N, 1> row_type;
  typedef Slice<Value_T, Rows_N, Cols_N> col_type;
  const row_type row(unsigned int n) const {
    assert(n < rows);
    return row_type(data + (n * Cols_N));
  }
  const col_type column(unsigned int n) const {
    assert(n < cols);
    return col_type(data + n);
  }
  row_type operator[](unsigned int n) { return row(n); }
  const row_type operator[](unsigned int n) const { return row(n); }
  const static unsigned int rows = Rows_N;
  const static unsigned int cols = Cols_N;
  typedef Value_T value_type;
private:
  mutable Value_T data[Rows_N * Cols_N];
};

template<class Matrix1, class Matrix2>
struct MatrixMultiplicationType {
  typedef Matrix<typename Matrix1::value_type, Matrix1::rows, Matrix2::cols> type;
};

template<class Slice1_T, class Slice2_T>
typename Slice1_T::value_type dot_product(Slice1_T x, Slice2_T y) {
  assert(x.size() == y.size());
  return std::inner_product(x.begin(), x.end(), y.begin(), typename Slice1_T::value_type(0));
}

template<class Matrix1, class Matrix2>
typename MatrixMultiplicationType<Matrix1, Matrix2>::type
matrix_multiply(const Matrix1& m1, const Matrix2& m2)
{
  typename MatrixMultiplicationType<Matrix1, Matrix2>::type result;
  assert(Matrix1::cols == Matrix2::rows);
  for (int i=0; i < Matrix1::rows; ++i)
  for (int j=0; j < Matrix2::cols; ++j)
    result[i][j] = dot_product(m1.row(i), m2.column(j));
  return result;
}

template<typename Matrix_T>
void output_matrix(const Matrix_T& x) {
  for (int i=0; i < x.rows; ++i) {
    for (int j=0; j < x.cols; ++j) {
      std::cout << x[i][j] << ", ";
    }
    std::cout << std::endl;
  }
  std::cout << std::endl;
}

void test_matrix() {
  Matrix<int, 1, 2> m1;
  m1[0][0] = 2;
  m1[0][1] = 3;
  Matrix<int, 2, 2> m2;
  m2[0][0] = 2;
  m2[0][1] = 0;
  m2[1][0] = 0;
  m2[1][1] = 2;
  output_matrix(m2);
  Matrix<int, 1, 2> m3 = matrix_multiply(m1, m2);
  output_matrix(m3);
}

int main() {
  test_matrix();
  system("pause");
  return 0;
}


Christopher Diggins

Posts: 1215
Nickname: cdiggins
Registered: Feb, 2004

Re: Storing and Tracking Matrix Dimensionality at Compile-Time Posted: May 15, 2005 2:20 AM
Reply to this message Reply
I just compared benchmarks to boost::ublas and preliminary results show that it runs at over twice as fast for integer matricies. Here is my benchmark code:

using namespace boost::numeric::ublas;
 
template<class Number, unsigned int M, unsigned int R, unsigned int N, unsigned int I>
void compare_naive_and_ublas()
{
  Matrix<Number, M, N> result1;
  matrix<Number> result2;
  {
    std::cout << "about to run Naive test " << std::endl;
    Matrix<Number, M, R> m1;
    Matrix<Number, R, N> m2;
    int n = 0;
    for (unsigned i = 0; i < M; ++i)
    for (unsigned j = 0; j < R; ++j)
      m1[i][j] = ++n;
    n = 0;
    for (unsigned i = 0; i < R; ++i)
    for (unsigned j = 0; j < N; ++j)
      m2[i][j] = ++n;
    boost::progress_timer t;
    for (int i=0; i < I; i++) {
      result1 = matrix_multiply(m1, m2);
    }
    std::cout << "naive time elapsed " << t.elapsed() << std::endl;
  }
  {
    std::cout << "about to run uBlas test " << std::endl;
    matrix<Number> m1(M, R);
    matrix<Number> m2(R, N);
    int n = 0;
    for (unsigned i = 0; i < M; ++i)
    for (unsigned j = 0; j < R; ++j)
      m1(i, j) = ++n;
    n = 0;
    for (unsigned i = 0; i < R; ++i)
    for (unsigned j = 0; j < N; ++j)
      m2(i, j) = ++n;
    boost::progress_timer t;
    for (int i=0; i < I; i++) {
      result2 = prod(m1, m2);
    }
    std::cout << "boost time elapsed " << t.elapsed() << std::endl;
  }
  for (unsigned i = 0; i < M; ++i)
  for (unsigned j = 0; j < N; ++j)
    assert(result1[i][j] == result2(i, j));
}
 
int main() {
  compare_naive_and_ublas<int, 17,23,31,1000>();
  system("pause");
  return 0;
}

Flat View: This topic has 1 reply on 1 page
Topic: More on languages and objects Previous Topic   Next Topic Topic: Working Around Non-Virtual Destructors


Sponsored Links



Google
  Web Artima.com   

Copyright © 1996-2014 Artima, Inc. All Rights Reserved. - Privacy Policy - Terms of Use - Advertise with Us