## Concept-based matrix library

[21st October 2007]

A couple of years ago, I started to put together a computer game but I haven't had enough time to work on it, recently. I've had a little more spare time as of late, so I think it's about time I get going again. I managed to dig up a video of my old game, NunsWithGuns^{[1]}. Click on the image below to see a DivX video (1.7Mb):

I got rather nostaligic watching it, but I'm not sure whether I'll continue with this game or start a new one. I've learnt a lot over the last year or two and I have a bunch of ideas that at least warrant plenty of experimentation, if not a complete re-write.

One such idea involves upgrading my matrix and vector maths library^{[2]}. The existing library does the job but I can't help feeling that I'm missing a trick or two. So I've decided to experiment with an upgrade so that the facilities I define can be used in conjunction with existing libraries, such as Blitz++.

### The aims

The library needs to very flexible. I'm particularly interested in it having the following characteristics:

- generic
- very efficient, or at least implemented so that there is ample room for future low-level optimisation
- compatible with common graphics APIs, especially OpenGL
- natural syntax
- should be able to interoperate with existing libraries
- portable

That's a pretty tall order. Some of these items may at first appear to conflict with one another, but if a concept-based approach is adopted, we should be able to meet all of these goals.

### What does a concept-based approach

actually entail?

Concepts will be formally introduced in to the next C++ standard as a first class citizen of the language, but people have been simulating them in an adhoc fashion for a long time, now. Even the existing C++ standard effectively defines a number of concepts related to iterators. You might have heard of the distinction between InputIterator, ForwardIterator, RandomAccessIteraror and so on? These can all be thought of as concepts, in the sense that I'll be using in this post.

A concept is a way of bundling together a bunch of properties that hold true about a set of types. If a particular C++ type has all the given properties, then it is said to be a *model of the concept*.

Let's take the `std::sort()`

algorithm as an example. This function takes a pair of iterators and a function-object to be used to compare two elements of the sequence.

The iterators however, must be of a type that models the RandomAccessIterator concept. You can't sort the elements of an `std::list`

using `std::sort`

because its iterators do not model the RandomAccessIterator concept^{[3]} because its iterators can't be incremented by an arbitrary amount in constant time.

`std::list`

iterators instead model the BidirectionalIterator concept. The set of properties needed of a BidirectionalIterator model are a subset of those needed of a RandomAccessIterator model. In this situation, we say that RandomAccessIterator is a refinement of the BidirectionalIterator concept. Any type that models RandomAccessIterator also models BidirectionalIterator.

So a number of the C++ standard library algorithms are only callable with types that model a particular refinement of the InputIterator and/or OutputIterator concepts. The requirements on the types of arguments passed to generic function, such as a standard library algorithm can only be documented externally to the code. Currently, there is nothing in the language to explicitly check that the type of the iterators given to `std::sort()`

do indeed model the RandomAccessIterator concept. The first time you know that you have a concept violation

is when you get a horrendously long compiler error.

However, as I mentioned earlier, concepts are being added to the next revision of the C++ standard as part of the core language. You will be able to add additional decoration to your generic code in order to catch concept violations early, so that better error messages can be given.

GCC has already implemented the feature, if you're interested in playing around with it. In another post, I linked to a video that summaries the concepts feature and the benefits it will bring.

But for this library, I'm just going to rely on documentation to specify my concepts. Perhaps when the feature is widely implemented, I can add C++0x concepts support.

### Why a concept-based library?

In short, a concept-based approach is needed to have any chance of fulfilling my aims for the library. In particular, the ability to operate with existing matrix types coupled with the need for an efficient and generic implementation pretty much seals the deal. Let's have a look at why this is.

First of all I should mention that this library is only going to be dealing with matrices of fixed size and those sizes are available at compile time. So we might expect to be able to use a 2D array of doubles as a matrix. Indeed we will be able to do just that with this approach.

All of the code I'll present in this post will live inside the `mafs`

namespace. I'll probably change it to something else in time. It'll do for now, though. I mention it now because I'll be omitting it for brevity in this post.

The basic interface of the default matrix type provided by the library will look something like this sketch:

```
template<typename Element, typename Height, typename Width>
class matrix
{
public:
typedef Element element_type;
typedef unsigned size_type;
static const size_type height = Height;
static const size_type width = Width;
// ...
element_type &operator() (size_type i, size_type j);
const element_type &operator() (size_type i, size_type j) const;
// ...
};
```

So we can access the element at (i, j) in a `matrix m`

, with the syntax `m(i, j)`

. We can also find out various pieces of information about the `matrix`

's size and types by looking at the class-scope `typedefs`

and constants. But this is not the case with a 2D array of doubles, for example. So in order to treat these two entities uniformly, we'll say that any type to be used as a matrix will have to specialize the a handful of templates provided by the library. Something like:

```
template<typename T>
struct element_type { typedef typename T::element_type type; };
template<typename T>
struct size_type { typedef typename T::size_type type; };
template<typename Matrix>
struct matrix_access
{
typedef typename size_type<Matrix>::type size_type;
static typename element_type<Matrix>::type get(const Matrix&m, size_type i, size_type j)
{
return m(i, j);
}
};
```

You'll notice that the default definition will work with the matrix template class outlined above. To make it work with 2D arrays, however, we'll need to add a specialisations:

```
template<typename T, std::size_t H, std::size_t W>
struct element_type<T[H][W]> { typedef T type; };
template<typename T, std::size_t H, std::size_t W>
struct size_type<T[H][W]> { typedef std::size_t type; };
template<typename T, std::size_t H, std::size_t W>
struct matrix_access<T[H][W]>
{
typedef typename size_type<T[H][W]>::type size_type;
static const typename element_type<T[H][W]>::type &
get(const T m[H][W], size_type i, size_type j)
{
return m[i][j];
}
};
```

Now matrix elements can be accessed uniformly by writing `matrix_access::get(m, i, j)`

. This is quite verbose, but it doesn't matter as we'll only be using this syntax in the guts of the library. The client won't have to use it.

Now, if you have your own matrix type, of fixed size, you can define your own specialisations as we've done above for 2D arrays and it will then work seamlessly with the library.

This allows us to write efficient, generic code that can interoperate with numerous other libraries, without having to butcher their existing code.

### The key concepts within the library

We can use the above as a basis for one of the library's concepts, of which there will be three, each being a refinement of the one before it. Here are the names of the concepts we'll define. The details will be fleshed out shortly afterwards.

*FixedMatrix*— models of this concept can be thought of as immutable 2D arrays*MutableFixedMatrix*— models of this concept can be thought of as mutable 2D arrays. Clearly, mutable 2D arrays can be treated as immutable 2D arrays, so this really is a refinement of FixedMatrix*NaturalMatrix*— models of this concept can be thought of as mutable 2D arrays that support element access through the`m(i, j)`

syntax.

The FixedMatrix concept is what we'll use internally for the majority of arithmetic operations, such as multiplication. But when we have operations that change the elements of a matrix in-situ, we will require that a model of the MutableFixedMatrix concept is given.

The NatualMatrix concept is what we'll want to use most in the outside world

, because of its convenient and natural element-access syntax.

So I'll quickly run through the concepts now in a bit more detail, fleshing out what it means to be a model of each.

### FixedMatrix

Given a type `M`

and an object `m`

of that type, we say that `M`

is a model of the FixedMatrix concept if the following conditions are met:

- A
`typedef M::element_type`

is exposed that is the type of each matrix element in an object of type`M`

, else`element_type<>`

is specialised such that`element_type<M>::type`

yields such a type - A
`typedef M::size_type`

for an integral type is exposed, instances of which should be used as indices when accessing the elements of an object of type`M`

, else`size_type<>`

is specialised such that`size_type<M>::type`

yields such a type - A
`typedef M::operand_category`

is exposed which is the same type as`matrix_tag`

, else`operand_category<>`

is specialised such that`operand_category<M>::type`

yields such a type (we'll use this to distinguish between matrices, vectors and scalars in mixed operations such as multiplication.`matrix_tag`

is just an empty struct) `M::height`

is a value available at compile time which represents the height of objects of type`M`

and is of a type that is convertible to`size_type<M>::type`

, else`matrix_height<>`

is specialised such that`matrix_height<M>::value`

yields such a value`M::width`

is a value available at compile time which represents the width of objects of type`M`

and is of a type that is convertible to`size_type<M>::type`

, else`matrix_width<>`

is specialised such that`matrix_width<M>::value`

yields such a value- Either
`matrix_access<M>::get_const(m, i, j)`

else`m(i, j)`

returns the element at position (i, j) in`m`

, where`i`

and`j`

are objects of type convertible to`size_type<M>::type`

. Any type of object may be returned so long as it is convertible to`element_type<M>::type`

e.g.

```
extern M m;
extern element_type<M>::type x;
extern size_type<M>::type i, j;
x = m(i, j); // If this is not valid, the next line must be
x = matrix_access<M>::get_const(m, i, j);
```

Note that when we define a new matrix type to work with the library, we can add class-scope `typedefs`

to represent the dimensions and various types associated with the matrix. However, we can use existing types with the library without changing them by defining a handful of template specialisations. This traits-based technique is quite closely related to the notion of a *concept map*, part of the concepts proposal for C++0x.

If you look back at the primary template for `element_type<>`

, you'll see that it looks for an class-scope `typedef`

by default. Thus any type that exposes this `typedef`

will automatically work with the more generic `element_type<M>::type`

. The same goes for the other properties needed by the concept. This way, we can use a uniform syntax inside the code of the library to access matrix properties.

There are some extra details that we should strictly define, such as the properties needed of `element_type<M>::type`

(i.e. it should act like

a floating point numeric type), but we'll gloss over these for now.

### MutableFixedMatrix

Given a type `M`

and an object `m`

of that type, we say that `M`

is a model of the MutableFixedMatrix concept if the following conditions are met:

`M`

is a model of the FixedMatrix concept- For non-
`const m`

, either`matrix_access<M>::get_mutable(m, i, j)`

else`m(i, j)`

returns a reference to the element at position (i, j) in`m`

, where`i`

and`j`

are objects of type convertible to`size_type<M>::type`

e.g.

```
extern M m;
extern element_type<M>::type x;
extern size_type<M>::type i, j;
m(i, j) = x; // If this is not valid, the next line must be
matrix_access<M>::get_mutable(m, i, j) = x;
```

### NaturalMatrix

Given a type `M`

and an object `m`

of that type, we say that `M`

is a model of the NaturalMatrix concept if the following conditions are met:

`M`

is a model of the MutableFixedMatrix concept- If
`i`

and`j`

are objects convertible to`size_type<M>::type`

,`m(i, j)`

must return the same object or reference as`matrix_access<M>::get_mutable(m, i, j)`

if`m`

is non-`const`

or`matrix_access<M>::get_const(m, i, j)`

if`m`

is`const`

Now we have the concepts defined, let's have a look at how we can put them in to action.

### Example 1: matrix transpose

```
template<typename FixedMatrix>
matrix
<
typename element_type<FixedMatrix>::type,
matrix_width<FixedMatrix>::value,
matrix_height<FixedMatrix>::value
>
transpose(const FixedMatrix &m)
{
typedef FixedMatrix M; // for brevity!
typedef typename size_type<M>::type sz_t;
matrix
<
typename element_type<M>::type,
matrix_width<M>::value,
matrix_height<M>::value
>
ret;
for (sz_t i = 0; i != matrix_height<M>::value; ++i)
for (sz_t j = 0; j != matrix_width<M>::value; ++j)
ret(j, i) = at(m, i, j);
return ret;
}
```

You'll no doubt have noticed the `at()`

function, here, which is simply a short-hand for `matrix_access<>::get_const()`

or `matrix_access<>::get_mutable()`

depending on whether its first argument is `const`

or not.

The return type is somewhat long-winded, but that's the price you pay for generality. The client code will still be concise. Here's the body of my test function, for example:

```
typedef mafs::matrix<double, 3, 3> mat3x3;
void test_transpose()
{
using mafs::transpose;
using mafs::fill;
mat3x3 m1;
fill(m1) = 1, 2, 3,
4, 5, 6,
7, 8, 9;
float m2[3][3];
fill(m2) = 1, 2, 3,
4, 5, 6,
7, 8, 9;
CHECK(m1 == m2);
mat3x3 t1(transpose(m1));
mat3x3 t2(transpose(m2));
CHECK(t1 == t2);
CHECK(t1(0, 0) == 1);
CHECK(t1(0, 1) == 4);
CHECK(t1(0, 2) == 7);
CHECK(t1(1, 0) == 2);
CHECK(t1(2, 0) == 3);
CHECK(t2(0, 0) == 1);
CHECK(t2(0, 1) == 4);
CHECK(t2(0, 2) == 7);
CHECK(t2(1, 0) == 2);
CHECK(t2(2, 0) == 3);
}
```

I use a couple of little utility functions, here, most notably `fill()`

. I won't go in details, but all it does is allow the use of the concise comma-based syntax to fill a matrix. Perhaps I'll describe this in another post. It's a pretty standard technique for this kind of library, though.

But you can see that we're able to use both the matrix type supplied with the library and a 2D array uniformly. We can even compare two different models using `operator==`

. This works because the `matrix<>`

template class lives in the `mafs namespace`

, as does the `operator==`

I've defined, allowing Koenig lookup to kick in.

### Example 2: matrix multiplication

This example is a little more complex, because matrices can be multiplied by scalars and vectors, as well as other matrices. Allowances need to be made for this.

The `operator*`

provided by the library looks like this:

```
template<typename Lhs, typename Rhs>
typename detail::product
<
Lhs,
Rhs,
typename operand_category<Lhs>::type,
typename operand_category<Rhs>::type
>
::result_type
operator* (const Lhs &lhs, const Rhs &rhs)
{
typedef typename operand_category<Lhs>::type lhs_cat;
typedef typename operand_category<Rhs>::type rhs_cat;
return detail::product<Lhs, Rhs, lhs_cat, rhs_cat>::multiply(lhs, rhs);
}
```

This simply delegates its work to `detail::product<>::multiply()`

(`detail`

is a nested namespace of `mafs`

). Note that the `operand_category`

s of the left and right operands are given as template parameters. This will allow the implementation of mixed scalar/vector/matrix multiplication in future through template specialisation. For now, we'll just consider matrix vs matrix.

Here's the rest of the code needed:

```
namespace detail
{
template<typename Lhs, typename Rhs, typename LhsCategory, typename RhsCategory>
struct product;
template<typename Lhs, typename Rhs>
struct product<Lhs, Rhs, matrix_tag, matrix_tag>
{
typedef typename element_type<Lhs>::type lhs_element;
typedef typename element_type<Rhs>::type rhs_element;
typedef typename detail::widest_type<lhs_element, rhs_element>::type element_type;
typedef matrix<element_type, matrix_height<Lhs>::value, matrix_width<Rhs>::value> result_type;
static result_type multiply(const Lhs &lhs, const Rhs &rhs)
{
BOOST_STATIC_ASSERT(matrix_width<Lhs>::value == matrix_height<Rhs>::value);
typedef result_type::size_type sz_t;
result_type ret;
for (sz_t i = 0; i != result_type::height; ++i)
for (sz_t j = 0; j != result_type::width; ++j)
for (sz_t k = 0; k != matrix_width<Lhs>::value; ++k)
ret(i, j) += at(lhs, i, k) * at(rhs, k, j);
return ret;
}
};
} // close namespace detail
```

As it happens, this is nothing like the actual multiplication code I've implemented, but it illustrates how things work. We can now multiply 2D arrays and `mafs::matrix<>`

objects, or any other objects whose type models the FixedMatrix concept.

I'll be discussing the actual multiplication code in my next post, as it's highly experimental and warrants it's own discussion. If you're desperate to know, have a look in the code archive available for download at the end of this post.

I've also cheated a little by using an as-yet-undefined `widest_type<>`

meta-function. This simply yields the largest

of the two types it is given, so that if you're multiplying a matrix of `float`

s with a matrix of `double`

s, you'll end up with a matrix of `double`

s to maintain precision.

### Example 3: matrix views

The concepts-based approach allows us to make models of FixedMatrix that reference existing matrices in order to provide a restricted or modified view of those matrices. Specific examples include:

- A view that looks at a smaller area of an existing matrix
- A view that looks at an existing matrix, but with a row and column removed (useful for calculating determinants)
- A view that presents the transpose of an existing matrix

In fact, our transpose function could instead be defined as follows:

```
template<typename FixedMatrix>
transposed_view<FixedMatrix> transpose(const FixedMatrix &m)
{
return transposed_view<FixedMatrix>(m);
}
```

where transposed_view<> would look something like this:

```
template<typename FixedMatrix>
class transposed_view
{
public:
typedef typename element_type<FixedMatrix>::type element_type;
typedef typename size_type<FixedMatrix>::type size_type;
typedef operand_category matrix_tag;
static const size_type width = matrix_height<FixedMatrix>::value;
static const size_type height = matrix_width<FixedMatrix>::value;
transposed_view(const FixedMatrix &m) : m_(&m) { }
element_type operator() (size_type i, size_type j) const { return at(*m_, j, i); }
private:
const FixedMatrix *m_;
};
```

And this might not be as silly as it sounds. In the previous version of our transpose function, we return a new `matrix`

object by value. A temporary gets created containing however many elements are necessary. The `transposed_view<>`

version, however, merely references an existing matrix and while it still generates a temporary, it is much smaller. So it's possible we'd see a performance improvement with this version. Of course this is entirely conjecture at this point and I'd need to benchmark to find out which really performs better. But the point is that the library simply doesn't care which method is used as they each return an object whose type is a model of the FixedMatrix concept.

### Interoperability with common graphics APIs

If I'm going to be using this library for games development, it had better work nicely with existing graphics, physics and collision APIs. My favourite graphics library, OpenGL, expects the storage of its matrices to be column-major by default. But wouldn't you know it, the Direct3D fixed function pipeline requires row-major matrices!

The `matrix<>`

template class is the lynch-pin of the library, the default model of the NaturalMatrix concept. This type at least should allow either storage convention. Direct access to the underlying array of elements also needs to be allowed, else the storage convention is entirely irrelevant! Fortunately, allowing either storage convention is a very simple task, as we'll soon see. The first thing we'll need to do is to add a fourth template parameter to the `matrix<>`

template class in order to specify the layout of the contained elements:

```
enum matrix_layout
{
row_major = 0,
column_major = 1
};
template<typename Element, unsigned Height, unsigned Width, matrix_layout Layout = column_major>
class matrix
{
public:
typedef Element element_type;
typedef unsigned size_type;
static const size_type height = Height;
static const size_type width = Width;
// ...
private:
element_type a_[height * width];
};
```

Since we know the layout of the elements array at compile time, we can dispatch element look-up to an appropriate function using template specialisation. For example, the non-const `operator()`

member operator can be implemented like this:

```
namespace detail
{
template<unsigned Height, unsigned Width, unsigned Layout>
struct index;
template<unsigned Height, unsigned Width>
struct index<Height, Width, column_major>
{
static unsigned calc(unsigned row, unsigned col)
{
assert(row < Height);
assert(col < Width);
return col * Height + row;
}
};
template<unsigned Height, unsigned Width>
struct index<Height, Width, row_major>
{
static unsigned calc(unsigned row, unsigned col)
{
assert(row < Height);
assert(col < Width);
return row * Width + col;
}
};
}
template<typename Element, unsigned Height, unsigned Width, matrix_layout Layout = column_major>
class matrix
{
private:
typedef detail::index<Height, Width, Layout> index;
public:
// ...
element_type &operator() (size_type i, size_type j)
{
return a_[index::calc(i, j)];
}
// ...
};
```

So whenever an element of the array is accessed using a 2D index, we make sure that the implementation goes through the `index::calc()`

function. Now all that needs to be done is to add a way to access this internal array. This is a simple matter of adding a member function. I won't bore you that, here.

In fact, I may generalise this by introducing an additional concept in future for matrices that allow access to their internal array of elements. Something like ArrayFixedMatrix or AccessibleFixedMatrix, perhaps?

### Room for optimisation

I haven't done anything particular in this post to optimise any of the operations we've defined. But there's ample room for template specialisation in order to vectorise operations and so forth, or maybe even use something like OpenMP for larger matrices.

In my next post I'll be presenting some techniques to optimize matrix multiplication chains and provide hooks to enable user-defined optimizations for one's own matrix types. This is already implemented in the version of code downloadable at the end. It's rather still experimental in so far as it's only only makes the code faster under some circumstances according to my benchmarks and I may roll-back the implementation to a more straight-forward one. Client code is none the wiser, though, as the interface is the same. But this library is a good example of templates being able to provide both generality and performance. Portability

This code should work with any standards conforming compiler, but I'll only be supporting g++ and MSVC 8. We'll see in the next post that I'll be employing some really rather heavy meta-programming techniques to enable these multiplication optimisations I keep mentioning, so I won't be able to support the current borland and digital mars compilers. This is portable enough for me, anyhow.

### Still todo

The library is nowhere near finished, or ready for use in the real world. In particular, I still need to implement the following things:

- The remaining matrix arithmetic operations e.g. addition, subtraction, etc
- The remaining matrix functions e.g.
`inverse()`

,`determinant()`

, etc - Integration of scalar and vector concepts
- Various matrix and vector views e.g. sub-matrix view, row/column view, etc
- Research in to vectorisation and optimisations
- It would be nice to provide a header that enables Blitz++ compatibility once
`#include`

d

But once the library is in a good state, I'll make it available on this web site.

### Conclusion

Even though the library is nowhere near finished, I'm confident that the concept-based approach is the right one. I certainly think all my goals for the library can be met by heading down this road. In fact I'm thinking about applying the concept-based approach in a bunch of other places within my game, particularly the geometry library. There might by concepts such as Polygon, MutablePolygon, FixedPolygon, MutableFixedPolygon and so on.

Of course, the concept-based approach is complicated and does require a high quality compiler, so it might not be to everyone's liking.

### Downloads

##### Footnotes

### Comments

All original content copyright© Edd Dawson.

Any opinions expressed by Edd are his own and are not necessarily shared by his employer. Or by anyone else, in fact.

All source code appearing on this website that was written by Edd Dawson is made available under the terms of the Boost software license version 1.0 unless otherwise stated or implied by the license associated with the work from which the code is derived.

Brad

[27/10/2007 at 05:29:00]

Wow, and I thought my math library was cool...