Yet another matrix class in C++

Here is an atempt to provide a simple and correct matrix template class in C++. It is based on the great exemple from the C++ FAQ Lite.

Important note: this class is provided only for educational purpose, and is not meant to be used in a production environment. If ou're looking for production-ready classes, have a look on the Boost C++ libraries. E.g., this class is not optimized for speed by using expression templates.

It provides:

• arithmetic operators (+, -, /, *)
• matrix inversion method (be carefull, the matrix is then replaced by its inverse)
• stream operators (>> and <<)
• is_equal method, with special treatment for float and double
• brackets operator (), to allow things like
  Matrix<int> A(3,3); // create a 3x3 matrix of int
A(1,2) = -5; // put -5 in the 2nd row, 3rd column of A (since indexes starts at 0)

All of this in a template class that can be used with all nulber-like entities matching te following criteria:

• operators +, -, *, =, /, ==, !=, << and >> are provided
• autocast from double is available
• entity has a default constructor
• entity supports std::abs (provides an overload)

It also provides I/O operators with the following syntaxes:

• "[nrows x ncols]" and one line per row of tab-separated values for the output operator <<
• blank-separated (spaces, tabs or new lines) values, row-by-row for the input operator >> (it is convenient to use the same syntax as output, ie line-rows of tab-separated values, but it's not mandatory)

Example of use:
#include "matrix.h"

using namespace std;
using namespace mtx;

int main( int argc, char* argv[] )
{
    Matrix<double> A (3, 3);

    for ( unsigned i = 0; i < m.nrows(); ++ i )
        for ( unsigned j = 0; j < m.cols(); ++ j )
            A(i,j) = i + j;

    cout << A << endl;
}


The matrix inversion method is derived from the corresponding C method of Numerical Recipes in C, Second Edition.

Compiles with gcc4, gcc4-apple and MSVS 2005.

Since it's a template class, all the code fits in a signle header file, matrix.h:

#ifndef MTX_MATRIX_H
#define MTX_MATRIX_H 1

#include <vector>
#include <iostream>
#include <string>
#include <stdexcept>

namespace mtx {

// generic function for equals
template<typename T>
inline bool is_equal( const T& x, const T& y )
{
  return ( x == y );
}

// Its specialization for doubles
template<>
inline bool is_equal<double>( const double& x, const double& y )
{
  static const double epsilon = std::numeric_limits<double>::epsilon();
  const double diff = std::abs( x - y );
  return ( diff <= epsilon * std::abs(x) || diff <= epsilon * std::abs(y) );
}
// and for floats
template<>
static inline bool is_equal<float>( const float& x, const float& y )
{
  static const float epsilon = std::numeric_limits<float>::epsilon();
  const float diff = std::abs( x - y );
  return ( diff <= epsilon * std::abs(x) || diff <= epsilon * std::abs(y) );
}

// The Matrix<T> class is a generic 2-dimensional container for numeric elements.
// It can calculate basic operations +, -, *, /, << and >>
// if those elements provide :
//        operators: +, -, *, =, /, ==, !=, <<, >>
//        autocast from double
//        a default constructor
//        an overload of std::abs
template<typename T>
class Matrix {
  //typdefs
  typedef std::vector<T> vector_t;
  typedef std::vector<vector_t> matrix_t;
public:
// Constructor
  Matrix(const unsigned nrows, const unsigned ncols, const T& init_val = 0);

  // No need for destructor, copy constructor or assignment operator
  // as std::vector properly handle them

// Accessors
 // Access methods to get the (i,j) element:
  T&       operator() (unsigned i, unsigned j);        // not-const subscript operator
  const T& operator() (unsigned i, unsigned j) const;  // const subscript operator

  unsigned nrows() const;  // returns the number of rows in the matrix
  unsigned ncols() const;  // returns the number of columns in the matrix

// Exceptions
  // Thrown by constructor if either size is zero
  class BadSize : public std::invalid_argument
  { public: BadSize() : std::invalid_argument("Bad size for matrix.") {} };
  // Thrown by accessors if i or j is too big
  class BoundsViolation : public std::out_of_range
  { public: BoundsViolation() : std::out_of_range("Matrix index out of range.") {} };
  // Thrown by operations throw if the dimentions of one operand are not suitable for the operation
  class BadOpSize : public std::invalid_argument
  { public: BadOpSize() : std::invalid_argument("Bad size of matrix operand.") {} };
  // Thrown by LU decomposition if the matrix is singular
  class Singular : public std::invalid_argument
  { public: Singular() : std::invalid_argument("Invalid argument: singular matrix.") {} };

// Functions
  Matrix<T>& invert();

private:
// Private functions called by invert()
  void lu_decomposition(std::vector<unsigned>& index);
  void lu_substitution(const std::vector<unsigned>& index, vector_t& col) const;

// Mermber variable
  matrix_t matrix;
};

template<typename T>
inline unsigned Matrix<T>::nrows() const
{
  return matrix.size();
}

template<typename T>
inline unsigned Matrix<T>::ncols() const
{
  return matrix[0].size();
}

// in the function definition we use (row, col) instead of (i,j) for clarity
template<typename T>
inline T& Matrix<T>::operator() (unsigned row, unsigned col)
{
  if ( row >= nrows() || col >= ncols() )
    throw BoundsViolation();
  return matrix[row][col];
}

template<typename T>
inline const T& Matrix<T>::operator() (unsigned row, unsigned col) const
{
  if ( row >= nrows() || col >= ncols() )
    throw BoundsViolation();
  return matrix[row][col];
}

//////////////////////////////////////////////////


// Replace M by the LU decompusition of a rowwise-permutation of itself, and
// takesa vector of unsigned int that will become the row-permutations index
template<typename T>
void Matrix<T>::lu_decomposition(std::vector<unsigned>& index)
{
  // We only perform LU decomposition on squared matrices
  if( ncols() != nrows() )
    throw typename Matrix<T>::BadOpSize();
 
  // hence n is the size of the Matrix
  unsigned n = ncols();

  register unsigned i, j, k; // for clarity we won't re-declare them every loop
  // for speed we declare them as register
  unsigned i_pivot;   // a row index used for permutation
  // other temporary objects declared here to avoid multiple contructor calls
  T tmp, max, sum;

  // This will store the the maximum of each row (to be used for the pivot search)
  std::vector<T> max_inv( n );

  index.resize( n ); // The row-permutation index has one entry for each row

  using std::abs;

  // Find the maximum of each row and store its inverse in inv_max
  for( i = 0; i < n; ++i ) {
    max = 0;
    for( j = 0; j < n; ++j ) {
      if( abs( matrix[i][j] ) > abs( max ) )
        max = matrix[i][j];
    }

    if( is_equal<T>(max, 0) ) // Largest element == 0 implies singular matrix
      throw typename Matrix<T>::Singular();
    max_inv[i] = max;
  }

  // Crout's method to solve N≤+N equations with N≤ unknowns :
  // remplace each element of the row by the lower or upper matrix elements
  for( j = 0; j < n ; ++j ) {
    for( i = 0; i < j; ++i ) {  // first part: 0 <= i < j
      sum = matrix[i][j];
      for( k = 0; k < i; k++ )
        sum = sum - matrix[i][k] * matrix[k][j];
      matrix[i][j] = sum;
    }
    max = 0;
    i_pivot = j; // not mandatory just to prevent compiler warning
    for( i = j; i < n; ++i ) {  // seconde part: j <= i < n (with pivoting)
      sum = matrix[i][j];
      for( k = 0; k < j; k++ )
        sum = sum - matrix[i][k] * matrix[k][j];
      matrix[i][j] = sum;
        // search for the largest element
      tmp = sum * max_inv[i];
      if( abs(tmp) >= abs(max) ) {
        max = tmp;
        i_pivot = i;
      }
    }
    
    if(j != i_pivot) {                    // if we need to permute row
      std::swap(matrix[j], matrix[i_pivot]);  // swap vectors is fast
      max_inv[i_pivot] = max_inv[j];          // actualize the maximum inverse
    }
    index[j] = i_pivot;       // index the column number

    if( is_equal<T>(matrix[j][j], 0) )   // Largest element == 0 implies singular matrix
      throw typename Matrix<T>::Singular();
    
    // Divide the remaining elements in the columns by the pivot
    if( j < n-1 ) {
      tmp = 1.0 / matrix[j][j];   // only one division, several multiplications
      for(i = j+1; i < n; i++)
        matrix[i][j] = matrix[i][j] * tmp;
    }
  } // next column j
}

// One the matrix is LU-decomposed, applies forward & backward substitution
// Optimized for vectors with many zeroes
template<typename T>
void Matrix<T>::lu_substitution(const std::vector<unsigned>& index, vector_t& col) const {
  T sum;
  unsigned n = ncols();
  register unsigned i, j;
  bool has_non_zero = false;
  unsigned i_non_zero = 0;

  // forward substitution + permutations
  for( i = 0; i < n; ++i ) {
    // applies the permutation to solution vector
    sum = col[index[i]];
    col[index[i]] = col[i];
    // forward substitution for non-zero columns
    if( has_non_zero ) {
      for( j = i_non_zero; j < i; ++j )
        sum = sum - matrix[i][j] * col[j];
    }
    else if( ! is_equal<T>(sum, 0) ) {   // if we encounter a non-zero element
      has_non_zero = true;    // we'll apply the forward substitution for the next rows
      i_non_zero = i;         // starting from the current row number
    }
    col[i] = sum;     // place the solution in the vector;
  }
  // back substitution
  for( i = n; i >= 1; --i ) {
      sum = col[i-1];
    for( j = i; j < n; ++j )
      sum = sum - matrix[i-1][j] * col[j];
    col[i-1] = sum / matrix[i-1][i-1];
  }
}

// Replace the current matrix by its inverse and returns a reference to itself.
template<typename T>
Matrix<T>& Matrix<T>::invert() {
  unsigned n = ncols();
  std::vector<unsigned> index( n );
  matrix_t identity( n, vector_t( n, 0 ) ); // Note : identify is stored by COLUMN, and NOT by ROW like the current matrix
  // Perform LU decomposition only once
  try {
    lu_decomposition(index);
  }
  catch(...) {
    throw;
  }
  // Perform the substitution for each column of the identity matrix
  // ie identity becomes the inverse of the current matrix
  for( unsigned i = 0; i < n; ++i ) {
    identity[i][i] = 1;       // place the 1 in the diaginal
    lu_substitution(index, identity[i]);
  }

  // We can now exchange our values with the identity matrix which
  // has been replaced by our inverse, by doing a simple transpose
  for( unsigned i = 0; i < n; ++i )
    for( unsigned j = 0; j < n; ++j )
      matrix[i][j] = identity[j][i];

  return (*this);
}

//////////////////////////////////////////////////

template<typename T>
Matrix<T>::Matrix(const unsigned nrows, const unsigned ncols, const T& init_val)
: matrix (nrows, vector_t(ncols, init_val)) // initialize all vectors
{
  if( nrows == 0 || ncols == 0 )
    throw BadSize();
  for( unsigned i = 0; i < nrows; ++i )
    matrix[i].resize(ncols);
}

// Matrix operations.
// In comments we use the notation: Matrix[number of rows * number of columns]

// Addition : A[m * n] + B[m * n] = C[m * n]
template<typename T>
Matrix<T> operator+( const Matrix<T>& A, const Matrix<T>& B )
{
  if( A.nrows() != B.nrows() || A.ncols() != B.ncols() )
    throw typename Matrix<T>::BadOpSize();

  Matrix<T> C( A.nrows(), A.ncols() );

  for( unsigned i=0, nrows=C.nrows(); i < nrows; ++i )
    for( unsigned j=0, ncols=C.ncols(); j < ncols; ++j )
      C(i,j) = A(i,j) + B(i,j);

  return C;
}

// Substraction : A[m * n] - B[m * n] = C[m * n]
template<typename T>
Matrix<T> operator-( const Matrix<T>& A, const Matrix<T>& B )
{
  if( A.nrows() != B.nrows() || A.ncols() != B.ncols() )
    throw typename Matrix<T>::BadOpSize();

  Matrix<T> C( A.nrows(), A.ncols() );

  for( unsigned i=0, nrows=C.nrows(); i < nrows; ++i )
    for( unsigned j=0, ncols=C.ncols(); j < ncols; ++j )
      C(i,j) = A(i,j) - B(i,j);

  return C;
}

// Matrix multiplication:
// A[m * p] * B[p * n] = C[m * n]
template<typename T>
Matrix<T> operator*( const Matrix<T>& A, const Matrix<T>& B )
{
  if( A.ncols() != B.nrows() )
    throw typename Matrix<T>::BadOpSize();

  unsigned m = A.nrows();
  unsigned n = B.ncols();
  unsigned p = A.ncols();
  // create the matrix of size m*n and initialize all elements to 0
  Matrix<T> C( m, n, 0 );

  for( unsigned i = 0; i < m; ++i )
    for( unsigned j = 0; j < n; ++j )
      for( unsigned k = 0; k < p; ++k )
        C(i,j) = C(i,j) + ( A(i,k) * B(k,j) );

  return C;
}

// Matrix division
// C[n * n] = A[n * n] / B[n * n] is equivalent to
// C = A * B.inverse()        except that only a copy of B is inverted
// The multiplication code is reapeated to gain -1 Matrix instance
template<typename T>
Matrix<T> operator/( const Matrix<T>& A, const Matrix<T>& B )
{
  if( A.nrows() != B.nrows() || A.ncols() != B.ncols() || A.nrows() != A.ncols() )
    throw typename Matrix<T>::BadOpSize();

  unsigned n = A.nrows();
  Matrix<T> C( n, n ), INVB = B;
  try {
    INVB.invert();
  }
  catch(...) {
    throw;
  }

  // multiply A by the inverse of B
  for( unsigned i = 0; i < n; ++i )
    for( unsigned j = 0; j < n; ++j )
      for( unsigned k = 0; k < n; ++k )
        C(i,j) = C(i,j) + ( A(i,k) * INVB(k,j) );

  return C;
}

//////////////////////////////////////////////////

// Write the matrix m in the output stream s.
// The first line is: [number_of_rows x number_of_columns]
// then it writes 1 line per matrix row, separating the values with tabs.
//Ex:     [3 x 3]
//        -1    5       0
//        7     -2      1
//        4     1       6
// This function need operator << to be available for the type T,
// otherwise it'll cause a compile-time error.
template<typename T>
std::ostream& operator<<( std::ostream& s, const Matrix<T>& M )
{
  s << "[ " << M.nrows() << " x " << M.ncols() << " ]\n";

  for( unsigned i=0, nrows=M.nrows(); i < nrows; ++i )
    for( unsigned j=0, ncols=M.ncols(); j < ncols; ++j )
      s << M(i,j) << (j == ncols-1 ? '\n' : '\t');

  return s;
}

// Read each elements of the matrix m from the input stream s.
// As '\t' and '\n' are white spaces (hence input delimiters), it is
// convenient to input the matrix as 1 line = 1 row of tab-separated values.
// Note: the dimensions of the matrix are already entered in m.
// This function needs operator >> to be available for the type T,
// otherwise it'll cause a compile-time error.
// The error handling is leaved to the underlying type, so on error all
// the remaining elements will keep their actual values, and the ios fail
// flag will be maintained, allowing a check like: while(cin >> m) {...}
template<typename T>
std::istream& operator>>( std::istream& s, Matrix<T>& M )
{
  for( unsigned i=0, nrows=M.nrows(); ( i < nrows && s); ++i )
    for( unsigned j=0, ncols=M.ncols(); ( j < ncols && s); ++j )
      if( !(s >> M(i,j)) )
        std::cerr << "Error while reading matrix, at row " << i << " column " << j << "." << std::endl ;

  return s;
}

} // end namespace mtx

#endif

If you have any question or comment regarding this code, feel free to email me.