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

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

 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 6:44 PM
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)
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 14, 2005 11:20 PM
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
 Previous Topic Next Topic