Heron-Centric: Ruminations of a Language Designer
Storing and Tracking Matrix Rows and Columns at Compile-Time
by Christopher Diggins
May 15, 2005
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;
}
```

# Talk Back!

If you'd like to be notified whenever Christopher Diggins adds a new entry to his weblog, subscribe to his RSS feed.