Diggins PDP #5: KMatrix (fixed row and column matricieis)

// Diggins Public Domain Post #5
// Contributed to Public Domain by Christopher Diggins
// This is an efficient matrix class implementation for when rows and
columns are known at compile-time

#include <algorithm>
#include <iterator>
#include <cassert>

//#include "static_assert.hpp"

template<class Value_T, int Size_N, int Step_N>
class kslice_iter
{
public:
// public typedefs
typedef Value_T value_type;
typedef Value_T& reference;
typedef const Value_T& const_reference;
typedef std::size_t size_type;
typedef std::ptrdiff_t difference_type;
typedef kslice_iter self;
// constructors
explicit kslice_iter(Value_T* x) : m(x) {
}
// operators
self& operator++() {
m += step();
return *this;
}
self operator++(int) {
self tmp = *this;
m += step();
return tmp;
}
Value_T& operator[](size_type n) {
assert(n < size());
return m[n * Step_N];
}
const Value_T& operator[](size_type n) const {
assert(n < N);
return m[n * Step_N];
}
static int size() {
return Size_N;
}
static int step() {
return Step_N;
}
Value_T& operator*() {
return *m;
}
friend difference_type operator-(const self& x, const self& y) {
return (y.m - x.m) / step();
}
friend bool operator==(const self& x, const self& y) {
return x.m == y.m;
}
friend bool operator!=(const self& x, const self& y) {
return x.m != y.m;
}
friend bool operator<(const self& x, const self& y) {
return x.m < y.m;
}
private:
// default constructor hidden from view
kslice_iter() { };
Value_T* m;
};

template<int N, class T, class RowIter_T, class ColIter_T>
T dot_product(RowIter_T x, ColIter_T y)
{
T ret(0);
for (int i=N; i > 0; --i) {
ret += *x++ * *y++;
}
return ret;
}

template<class T, class M1, class M2, class M3>
kmatrix_multiply(const M1& m1, const M2& m2, M3& result)
{
//STATIC_ASSERT(M1::cols == M2::rows);
//STATIC_ASSERT(M3::rows == M1::rows);
//STATIC_ASSERT(M3::cols == M2::cols);
for (int i=0; i < M1::rows; ++i) {
for (int j=0; j < M2::cols; ++j) {
result[i][j] = dot_product<M1::cols, T>(m1.row(i), m2.column(j));
}
}
}

template<class Value_T, int Rows_N, int Cols_N>
class kmatrix
{
public:
// public typedefs
typedef Value_T value_type;
typedef kmatrix self;
typedef Value_T* iterator;
typedef const Value_T* const_iterator;
typedef kslice_iter<Value_T, Cols_N, 1> row_type;
typedef kslice_iter<Value_T, Rows_N, Cols_N> col_type;
// public constants
static const int rows = Rows_N;
static const int cols = Cols_N;
// constructors
kmatrix() {
std::fill(begin(), end(), Value_T(0));
}
kmatrix(Value_T& x) {
std::fill(begin(), end(), x);
}
kmatrix(const self& x) {
operator=(x);
}
// public functions
self& operator=(const self& x) {
std::copy(x.begin(), x.end(), begin());
return *this;
}
self& operator=(Value_T& x) {
std::fill(begin(), end(), x);
return *this;
}
row_type row(int n) const {
assert(n < rows);
return row_type(m + (n * Cols_N));
}
col_type column(int n) const {
assert(n < cols);
return col_type(m + n);
}
row_type operator[](int n) {
return row(n);
}
const row_type operator[](int n) const {
return row(n);
}
iterator begin() {
return m;
}
iterator end() {
return m + size();
}
const_iterator begin() const {
return m;
}
const_iterator end() const {
return m + size();
}
static int size() {
return Rows_N * Cols_N;
}
// operators
self& operator+=(const self& x) {
for (int i=0; i < size(); ++i) {
m[i] += x.m[i];
}
return *this;
}
self& operator-=(const self& x) {
for (int i=0; i < size(); ++i) {
m[i] -= x.m[i];
}
return *this;
}
self& operator+=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] += x;
}
return *this;
}
self& operator-=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] -= x;
}
return *this;
}
self& operator*=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] *= x;
}
return *this;
}
self& operator/=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] /= x;
}
return *this;
}
self& operator=(value_type x) {
std::fill(begin(), end(), x);
return *this;
}
self operator-() {
self x;
for (int i=0; i < size(); ++i) {
x.m[i] = -m[i];
}
return x;
}
// friends
friend self operator+(const self& x, const self& y) {
return self(x) += y;
}
friend self operator-(const self& x, const self& y) {
return self(x) -= y;
}
friend self operator+(const self& x, value_type y) {
return self(x) += y;
}
friend self operator-(const self& x, value_type y) {
return self(x) -= y;
}
friend self operator*(const self& x, value_type y) {
return self(x) *= y;
}
friend self operator/(const self& x, value_type y) {
return self(x) /= y;
}
template<int Cols2_N>
friend kmatrix<Value_T, Rows_N, Cols2_N>
operator*(const self& x, const kmatrix<Value_T, Cols_N, Cols2_N>& y) {
kmatrix<Value_T, Rows_N, Cols2_N> ret;
kmatrix_multiply(x, y, ret);
return ret;
}
private:
mutable Value_T m[Rows_N * Cols_N];
};

// public domain matrix benchmark code

#include "kmatrix.hpp"

#include <cassert>
#include <iostream>
#include <cstdlib>
#include <boost/numeric/ublas/matrix.hpp>

class timer
{
public:
timer() : start(std::clock()) {
}
double elapsed() const {
return (double(std::clock() - start) / CLOCKS_PER_SEC) * 1000;
}
private:
std::clock_t start;
};

void output_time_elapsed(const timer& t, std::ostream& o = std::cout) {
o << static_cast<int>(t.elapsed()) << " msec elapsed" << std::endl;
}

struct reporting_timer {
~reporting_timer() {
output_time_elapsed(m);
}
timer m;
};


template<int Rows_N, int Cols_N, class Matrix_T>
void init_kmatrix(Matrix_T& m)
{
for (int i=0; i < Rows_N; ++i) {
for (int j=0; j < Cols_N; ++j) {
m[i][j] = i * Cols_N + j;
}
}
}

template<int Rows_N, int Cols_N, class Matrix_T>
void init_ublas_matrix(Matrix_T& m)
{
for (int i=0; i < Rows_N; ++i) {
for (int j=0; j < Cols_N; ++j) {
m(i, j) = i * Cols_N + j;
}
}
}

template<class T, int Iters_N, int M, int R, int N>
void kmatrix_performance_test()
{
std::cout << "kmatrix ";
kmatrix<T, M, R> km1;
kmatrix<T, R, N> km2;
kmatrix<T, M, N> km3;
init_kmatrix<M, R>(km1);
init_kmatrix<R, M>(km2);
{
reporting_timer t;
for (int i=0; i < Iters_N; ++i)
{
kmatrix_multiply<T>(km1, km2, km3);
}
}

std::cout << "ublas ";
boost::numeric::ublas::matrix<T> m1(M,R);
boost::numeric::ublas::matrix<T> m2(R,N);
boost::numeric::ublas::matrix<T> m3(M,N);
init_ublas_matrix<M, R>(m1);
init_ublas_matrix<R, M>(m2);
{
reporting_timer t;
for (int i=0; i < Iters_N; ++i)
{
boost::numeric::ublas::noalias(m3) =
boost::numeric::ublas::prod(m1,m2);
}
}

for (int i=0; i<M; ++i) {
for (int j=0; j<N; ++j) {
T x = m3(i,j);
T y = km3[i][j];
ASSURE(x == y);
}
}
}

void benchmarks()
{
std::cout << "comparing integer matricies" << std::endl;
std::cout << "2,2 X 2,2" << std::endl;
kmatrix_performance_test<int, 100000, 2, 2, 2>();
std::cout << "3,3 X 3,3" << std::endl;
kmatrix_performance_test<int, 100000, 3, 3, 3>();
std::cout << "100,100 X 100,100" << std::endl;
kmatrix_performance_test<int, 10, 100, 100, 100>();
std::cout << "100,1 X 1,100" << std::endl;
kmatrix_performance_test<int, 100, 100, 1, 100>();
std::cout << "1,100 X 100, 1" << std::endl;
kmatrix_performance_test<int, 100000, 1, 100, 1>();

std::cout << "comparing double matricies" << std::endl;
std::cout << "2,2 X 2,2" << std::endl;
kmatrix_performance_test<double, 100000, 2, 2, 2>();
std::cout << "3,3 X 3,3" << std::endl;
kmatrix_performance_test<double, 100000, 3, 3, 3>();
std::cout << "100,100 X 100,100" << std::endl;
kmatrix_performance_test<double, 10, 100, 100, 100>();
std::cout << "100,1 X 1,100" << std::endl;
kmatrix_performance_test<double, 100, 100, 1, 100>();
std::cout << "1,100 X 100, 1" << std::endl;
kmatrix_performance_test<double, 100000, 1, 100, 1>();
}

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

/*
my results:

comparing integer matricies
2,2 X 2,2
kmatrix 78 msec elapsed
ublas 172 msec elapsed
3,3 X 3,3
kmatrix 219 msec elapsed
ublas 437 msec elapsed
100,100 X 100,100
kmatrix 516 msec elapsed
ublas 1156 msec elapsed
100,1 X 1,100
kmatrix 125 msec elapsed
ublas 203 msec elapsed
1,100 X 100, 1
kmatrix 500 msec elapsed
ublas 1172 msec elapsed
comparing double matricies
2,2 X 2,2
kmatrix 79 msec elapsed
ublas 171 msec elapsed
3,3 X 3,3
kmatrix 250 msec elapsed
ublas 422 msec elapsed
100,100 X 100,100
kmatrix 656 msec elapsed
ublas 1344 msec elapsed
100,1 X 1,100
kmatrix 156 msec elapsed
ublas 219 msec elapsed
1,100 X 100, 1
kmatrix 703 msec elapsed
ublas 1172 msec elapsed
Press any key to continue . . .
*/