Fixed size matrix implementation

Question

I implemented a fixed-size matrix class which supports basic matrix operations (cannot use C++11). What do you think, and how could it be improved?

#ifndef FIXED_SIZE_MATRIX_H__
#define FIXED_SIZE_MATRIX_H__

#include <iostream>
#include <cstdlib>

template <typename Ty, int N, int M = N>
struct FixedSizeMatrix {
  typedef Ty value_type;

  union {
    struct {
      Ty element[N][M];
    };
    struct {
      Ty flatten[N * M];
    };
  };

  // Access with bounds checking
  Ty& operator()(int r, int c) {
    assert(r >= 0 && r < N);
    assert(c >= 0 && c < M);
    return element[r][c];
  }

  // Return matrix transpose
  FixedSizeMatrix<Ty, M, N> T() const {
    FixedSizeMatrix<Ty, M, N> result;
    for (int r = 0; r < N; ++r) {
      for (int c = 0; c < M; ++c) {
        result.element[c][r] = element[r][c];
      }
    }
    return result;
  }

  // Return matrix inverse
  FixedSizeMatrix<Ty, N, M> I();
};

///
// Matrix Inverse
///
// A design choice was made to keep FixedSizeMatrix an aggregate class to enable
// FixedSizeMatrix<float, 2, 2> = {1.0f, 1.0f, 0.5f, 0.2f} initialization.
// With c++11 it would be possible to create a FixedSizeMatrix base class
// and derive all variations from it, while retaining the brace
// initialization.
//
// Matrix inverse helpers
namespace detail {
// General version not implemented
template <typename Ty, int N, int M>
struct inverse;

// Matrix inversion for 2x2 matrix
template <typename Ty>
struct inverse<Ty, 2, 2> {
  FixedSizeMatrix<Ty, 2, 2> operator()(FixedSizeMatrix<Ty, 2, 2> a) {
    FixedSizeMatrix<Ty, 2, 2> result;
    Ty det =
        a.element[0][0] * a.element[1][1] - a.element[0][1] * a.element[1][0];
    assert(det != 0);

    result.element[0][0] = a.element[1][1] / det;
    result.element[1][1] = a.element[0][0] / det;
    result.element[0][1] = -a.element[0][1] / det;
    result.element[1][0] = -a.element[1][0] / det;
    return result;
  }
};
}  // detail

// Define matrix inverse
template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> FixedSizeMatrix<Ty, N, M>::I() {
  return detail::inverse<Ty, N, M>()(*this);
}

///
//  Matrix operations
///
// Matrix product
template <typename Ty, int N, int M, int P>
FixedSizeMatrix<Ty, N, P> operator*(FixedSizeMatrix<Ty, N, M> a,
                                    FixedSizeMatrix<Ty, M, P> b) {
  FixedSizeMatrix<Ty, N, P> result;

  for (int r = 0; r < N; ++r) {
    for (int c = 0; c < P; ++c) {
      Ty accum = Ty(0);
      for (int i = 0; i < M; ++i) {
        accum += a.element[r][i] * b.element[i][c];
      }
      result.element[r][c] = accum;
    }
  }
  return result;
}

// Unary negation
template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> operator-(FixedSizeMatrix<Ty, N, M> a) {
  FixedSizeMatrix<Ty, N, M> result;
  for (int e = 0; e < N * M; ++e) result.flatten[e] = -a.flatten[e];
  return result;
}

#define MATRIX_WITH_MATRIX_OPERATOR(op_symbol, op)                            \
  template <typename Ty, int N, int M>                                        \
  FixedSizeMatrix<Ty, N, M> operator op_symbol(FixedSizeMatrix<Ty, N, M> a,   \
                                               FixedSizeMatrix<Ty, N, M> b) { \
    FixedSizeMatrix<Ty, N, M> result;                                         \
    for (int e = 0; e < N * M; ++e)                                           \
      result.flatten[e] = a.flatten[e] op b.flatten[e];                       \
    return result;                                                            \
  }

MATRIX_WITH_MATRIX_OPERATOR(+, +);  
MATRIX_WITH_MATRIX_OPERATOR(-, -); 
#undef MATRIX_WITH_MATRIX_OPERATOR

#define MATRIX_WITH_SCALAR_OPERATOR(op_symbol, op)                          \
  template <typename Ty, int N, int M>                                      \
  FixedSizeMatrix<Ty, N, M> operator op_symbol(FixedSizeMatrix<Ty, N, M> a, \
                                               Ty scalar) {                 \
    FixedSizeMatrix<Ty, N, M> result;                                       \
    for (int e = 0; e < N * M; ++e)                                         \
        result.flatten[e] = a.flatten[e] op scalar;                         \
    return result;                                                          \
  }

MATRIX_WITH_SCALAR_OPERATOR(+, +);
MATRIX_WITH_SCALAR_OPERATOR(-, -);
MATRIX_WITH_SCALAR_OPERATOR(*, *);
MATRIX_WITH_SCALAR_OPERATOR(/, / );
#undef MATRIX_WITH_SCALAR_OPERATOR

template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> operator+(Ty scalar, FixedSizeMatrix<Ty, N, M> a) {
  return a + scalar;
}

template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> operator*(Ty scalar, FixedSizeMatrix<Ty, N, M> a) {
  return a * scalar;
}

template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> operator-(Ty scalar, FixedSizeMatrix<Ty, N, M> a) {
  return -a + scalar;
}

template <typename Ty, int N>
FixedSizeMatrix<Ty, N, N> identity_matrix() {
  FixedSizeMatrix<Ty, N, N> result = FixedSizeMatrix<Ty, N, N>();  
  for (int i = 0; i < N; ++i) result.element[i][i] = Ty(1);
  return result;
}

template <typename Ty, int N, int M>
std::ostream& operator<<(std::ostream& out, FixedSizeMatrix<Ty, N, M> a) {
  for (int r = 0; r < N; ++r) {
    for (int c = 0; c < M; ++c) {
      out << a.element[r][c] << " ";
    }
    out << std::endl;
  }
  return out;
}
#endif

Corbin · Accepted Answer · 2014-03-05 21:35:23Z

Looks good over all, but I have a few thoughts. A few are stylistic, some are technical, and some are just trade offs I want to make sure you're aware of.

I would use the value_type typedef for method return types, temporary types, etc. It should make for a bit prettier signature in IDEs and any kind of auto generated documentation.

I would try to be consistent with variable names with regards to N/M/r/c. I either use N/M/n/m or R/C/r/c. Some people like to avoid using different capitalizations though, so that might be a valid concern.

Template type parameters default to the current ones inside of a templated class, so all of the FixedSizeMatrix<Ty, M, N> types can be simplified to FixedSizeMatrix.

It's good that you're consistent, but 2 spaces is rather unusual for C++. I would use either 4 spaces or tabs.

You should also have a const version of operator().

If you ever decide to make the class more sophisticated, you might regret exposing element. Perhaps hide the struct as a private member, have operator() be unchecked, and have a checked version called at. You could even have a 1 parameter version of each one that would provide the flattened functionality.

This unfortunately would break aggregate initialization though, which I suspect is why you've made the decision to have it public.

Still though, even if you make keep it public, you can get in a habit of using operator() and at() and then if you ever use C++11, you can properly encapsulate it with minimal code breakage.

I would throw a std::out_of_range instead of using asserts in operator(). Note that this and having the checked version named at while using the operator as a pass through is consistent with the behavior of std::vector.

Same thing for your other asserts. Exceptions can be caught and handled if desired. Asserts not so much. if you're using asserts for performance reasons, please consider if the tiny difference really matters. Modern C++ compilers have a zero cost exception model where exceptions are only costly if they're actually thrown.

Don't use double underscores. You technically can in certain situations, but the rules are very difficult to remember. Your FIXED_SIZE_MATRIX_H__ is technically invalid.

Is I common notation for a matrix inverse (makes me think identity matrices)? Even if it is, I would name the method invert.

Same thing with the transpose.

Any time you're passing the matrix as an argument to a function and you're not modifying it, pass it by constant reference. There's no need to make a potentially very expensive copy just to read from it. There's also the option of having the argument copy but then using the copy to store the result (though that requires that you can operate on the matrix without a dependency cycle between elements). This would allow move semantics in C++11 on temporaries passed as the argument.

Everywhere you do Ty(0), you force an interface constraint on Ty in that it must have a constructor that is compatible with a single integer parameter. If you default construct the Ty, you avoid this interface constraint, though you of course add a new one. Consider the two cases this can affect: is a type used more likely to have a default constructor that represents 0 or a single integer constructor? I would lean towards the first one, oddly enough. For example, a Complex might not have an integer constructor.

It might be nice to put this in a namespace if you plan on the usage of this to grow to more than just 1 or 2 simple programs. That helps mentally organize your library as one coherent part, and it provides a few technical advantages. Namely, your detail namespace could conflict with someone else's ::detail namespace, and it would be nice to give your macros, though short lived, a prefix of some kind. It's relatively unlikely that anything in this would actually conflict, but if someone were using this with another matrix library, you never know. Namespaces don't cost anything to you as the developer (or to consumers of the class), but they offer a few advantages.

@BlazBratanic I don't know if I'd accept my answer just yet. There's a few other C++ regulars that might have some stuff to add. Might want to give it a day or two just in case :). — Corbin, Mar 5 '14 at 22:20

Yuushi · Answer 2 · 2014-03-06 03:12:13Z

There are better ways of defining operators than using macros that are only (slightly) more verbose:

#define MATRIX_WITH_MATRIX_OPERATOR(op_symbol, op)                            \
template <typename Ty, int N, int M>                                        \
FixedSizeMatrix<Ty, N, M> operator op_symbol(FixedSizeMatrix<Ty, N, M> a,   \
                                           FixedSizeMatrix<Ty, N, M> b) { \
  FixedSizeMatrix<Ty, N, M> result;                                         \
  for (int e = 0; e < N * M; ++e)                                           \
    result.flatten[e] = a.flatten[e] op b.flatten[e];                       \
  return result;                                                            \
}

Firstly, these parameters should be passed by const & (you're creating full copies of them currently). Secondly, you can easily do this with templates:

template <typename Ty, int N, int M, typename BinaryOp>
FixedSizeMatrix<Ty, N, M> base_op(const FixedSizeMatrix<Ty, N, M>& a,
                                  const FixedSizeMatrix<Ty, N, M>& b,
                                  BinaryOp f)
{
    FixedSizeMatrix<Ty, N, M> result;
    for(int e = 0; e < N * M; ++e) {
        result.flatten[e] = f(a.flatten[e], b.flatten[e]);
    }
    return result;
}

Your operators can then be one liners for the implementation:

#include <functional>

template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> operator+(const FixedSizeMatrix<Ty, N, M>& a,
                                    const FixedSizeMatrix<Ty, N, M>& b)
{
    return base_op(a, b, std::plus<Ty>());
}

template <typename Ty, int N, int M>
FixedSizeMatrix<Ty, N, M> operator-(const FixedSizeMatrix<Ty, N, M>& a,
                                    const FixedSizeMatrix<Ty, N, M>& b)
{
    return base_op(a, b, std::minus<Ty>());
}

And so on. This is slightly more verbose, but avoids macro ugliness (and will be just as efficient, once the compiler is done inlining everything).

Great point on the operators. Wish that had occurred to me :). — Corbin, Mar 6 '14 at 3:14
I did this modification and changed scalar operators the same. However, looking at the code with all the modifications it seems like I added large sections of boilerplate code. Are there any advantages of doing things this way, besides avoiding macros? — Blaz Bratanic, Mar 6 '14 at 8:03
@BlazBratanic It's mainly for avoiding macros. You might decide it's a worthwhile tradeoff to stick with macros here, but most of the boilerplate code is just function signatures. I think it's worthwhile to get rid of macros, personally. — Yuushi, Mar 6 '14 at 8:12
@Yuushi Thanks for your helpful suggestion. I would accept both answers if I could. — Blaz Bratanic, Mar 7 '14 at 8:56

asked	1 year ago
viewed	766 times
active	1 year ago

current community

your communities

more stack exchange communities

Fixed size matrix implementation

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged c++ matrix c++03 or ask your own question.

Hot Network Questions

current community

your communities

more stack exchange communities

Fixed size matrix implementation

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged c++ matrix c++03 or ask your own question.

Related

Hot Network Questions