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_i ter {

Slice_forward_i ter(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==(cons t Slice_forward_i ter& x) { return p == x.p; }

bool operator!=(cons t Slice_forward_i ter& x) { return p != x.p; }

Value_T* p;

};

typedef Slice_forward_i ter iterator;

typedef const Slice_forward_i ter 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 MatrixMultiplic ationType {

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(Sli ce1_T x, Slice2_T y) {

assert(x.size() == y.size());

return std::inner_prod uct(x.begin(), x.end(), y.begin(), typename

Slice1_T::value _type(0));

}

template<class Matrix1, class Matrix2>

typename MatrixMultiplic ationType<Matri x1, Matrix2>::type

matrix_multiply (const Matrix1& m1, const Matrix2& m2)

{

typename MatrixMultiplic ationType<Matri x1, 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<typena me Matrix_T>

void output_matrix(c onst 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(m 2);

Matrix<int, 1, 2> m3 = matrix_multiply (m1, m2);

output_matrix(m 3);

}

int main() {

test_matrix();

system("pause") ;

return 0;

}

--

Christopher Diggins

Object Oriented Template Library (OOTL)

http://www.ootl.org