Here is some code I wrote for Matrix multiplication for arbitrary

dimensionality known at compile-time. I am curious how practical it is. For

instance, is it common to know the dimensionality of matricies at

compile-time? Any help would be appreciated. Hopefully this code comes in

useful for someone, let me know if you find it useful, or if you have

suggestions on how to improve it.

// Public Domain by Christopher Diggins, 2005

#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

Object Oriented Template Library (OOTL)

http://www.ootl.org