#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;
}

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;
}