Skip to main content
GameDev.net gamedev.net
🔒 Locked

my matrix class

Started by tendifo Jul 9, 2006 at 10:06 PM 20 replies 2.5k views
Original Post
tendifo
tendifo
I'm working on a physics engine and I'm writing up a nice Matrix library to tag along with it. I'd like some feedback on it though. Before writing the bulk of my physics code, I'd like this to be a solid as possible because the more efficient this is, the easier it will be down the road. I've tested it out more or less and it seems to be working math wise, but if anyone sees an error, please point it out. Thanks! matrix.h 
#ifndef matrix_h
#define matrix_h

class Matrix {
public:
	// constructors
	Matrix(int r, int c);
	Matrix();

	// create the element array
	void initElements(int r, int c);

	// manipulate the matrix
	double getElement(int r, int c);
	double get(); // gets the first element in the matrix
	void setElement(int r, int c, double e);

	// retrieve information about the matrix
	int getRows();
	int getColumns();
	bool fetchError(); // returns if there's an error from last calculation or whatever
	
	// clearing functions
	void blank(); // create a blank matrix (all 0's)
	void identity(); // create an identity matrix

	// matrix math functions
	void add(Matrix a, Matrix b);
	void subt(Matrix a, Matrix b);
	void mult(Matrix a, Matrix b);
	void scale(double s);

private:
	int rows, columns; // number of rows and columns in the matrix
	double *elements; // the element array
	bool error; // contains information about errors
};

#endif
matrix.cpp 
#include "matrix.h"

Matrix::Matrix(int r, int c) {
	initElements(r, c);
}

Matrix::Matrix() {
	Matrix(1, 1);
}

void Matrix::initElements(int r, int c) {
	rows = r;
	columns = c;
	elements = (double*)malloc(sizeof(double) * rows * columns);
	blank(); // set everything to 0
}

double Matrix::getElement(int r, int c) {
	int elementId = r * columns + c;
	return elements[elementId];
}

double Matrix::get() {
	return getElement(0, 0):
}

void Matrix::setElement(int r, int c, double e) {
	int elementId = r * columns + c;
	elements[elementId] = e;
}

int Matrix::getRows() {
	return rows;
}

int Matrix::getColumns() {
	return columns;
}

bool Matrix::fetchError() {
	return error;
}

void Matrix::blank() {
	for (int i = 0; i < getRows(); i += 1)
		for (int j = 0; j < getColumns(); j += 1)
			setElement(i, j, 0);
}

void Matrix::identity() {
	if (getRows() == getColumns()) {
		blank();
		for (int i = 0; i < getRows(); i += 1)
			setElement(i, i, 1);
		error = false;
	} else {
		error = true;
	}
}

void Matrix::add(Matrix a, Matrix b) {
	// is subtraction legal?
	if (a.getRows() == b.getRows() && a.getColumns() == b.getRows()) {
		// recreate the matrix
		initElements(a.getRows(), a.getColumns());
		for (int i = 0; i < a.getRows(); i += 1)
			for (int j = 0; j < a.getColumns(); j += 1)
				setElement(i, j, a.getElement(i, j) + b.getElement(i, j));
		error = false;
	} else {
		error = true;
	}
}

void Matrix::subt(Matrix a, Matrix b) {
	// is subtraction legal?
	if (a.getRows() == b.getRows() && a.getColumns() == b.getRows()) {
		// recreate the matrix
		initElements(a.getRows(), a.getColumns());
		for (int i = 0; i < a.getRows(); i += 1)
			for (int j = 0; j < a.getColumns(); j += 1)
				setElement(i, j, a.getElement(i, j) - b.getElement(i, j));
		error = false;
	} else {
		error = true;
	}
}

void Matrix::mult(Matrix a, Matrix b) {
	// first check to see if multiplication is even legal
	if (a.getColumns() == b.getRows()) {
		// recreate the matrix
		initElements(a.getRows(), b.getColumns());
		for (int i = 0; i < a.getRows(); i += 1)
			for (int j = 0; j < b.getColumns(); j += 1) {
				double sum = 0;
				for (int k = 0; k < a.getColumns(); k += 1)
					sum += a.getElement(i, k) * b.getElement(k, j);
				setElement(i, j, sum);
			}
		error = false;
	} else {
		// however, if we're dealing with vector matrices, then we can multiply
		if (a.getColumns() == 1 && b.getColumns() == 1 && a.getRows() == b.getRows()) {
			// create a 1x1 matrix to store the result
			initElements(1, 1);
			double sum = 0;
			for (int i = 0; i < a.getRows(); i += 1)
				sum += a.getElement(i, 0) * b.getElement(i, 0);
			setElement(0, 0, sum);
			error = false;
		} else {
			error = true;
		}
	}
}

void Matrix::scale(double s) {
	for (int i = 0; i < getRows(); i += 1)
		for (int j = 0; j < getColumns(); j += 1)
			setElement(i, j, getElement(i, j) * s);
}
[Edited by - tendifo on July 9, 2006 10:29:42 PM]
yadango
yadango
i wouldn't have the default constructor initialize a 1x1 matrix. it seems like a wasteful allocation.
Saruman
Saruman
I totally agree with that statement as it is not only wasteful but hurts performance for absolutely no reason.
tendifo
tendifo
Yeah, good point.

I'll just keep it a null pointer then. And that actually gives me an idea of having boundry checks on the get and set functions.
Scythen
Scythen
Two comments.

First, having the default constructor initialize the matrix to identity is not entirely wasteful. Not initializing the matrix can result in hard to track down errors due to invalid floats (NAN and such). A better idea is to have a constructor that explicitly does not initialize the matrix. Say, a bool defaulted to true or something. This way someone won’t accidentally create and use an uninitialized or partially uninitialized matrix.

Second, I don’t know what platform this is but you may want to keep an internal representation that is aligned in some way, particularly for PC work. If the matrix is always aligned then you can use SSE math operations.
tendifo
tendifo
I don't fully understand what you mean by having an internal representation of the matrix and having it aligned and such.

I'm thinking it might be a better idea instead of having the add, subt, and mult functions alter the instance itself to create a bunch of static functions that would go something like this:

static bool add(Matrix *result, Matrix a, Matrix b);


Something like that. That way, the function would return if it went ok or not instead of changing an error bool in the instance.
Scythen
Scythen
If you like methods rather than operators (I don't) that fine but there is really no use for return values to indicate error conditions. Think about it, what error could there possibly be adding two matrices? Most matrix operations will never fail unless the matrix class fails to maintain a valid matrix. Obviously passing bad floats to math operations on the matrix can corrupt it but that’s not something you can check for and get any type of reasonable performance. You could have some ifdefs in there for a check build or something but that’s up to you. Also, you definitely don’t want a bunch of if statements in you math routines checking all the return values (if the programmer even bothers to).

My advice is have member methods or operators (I greatly prefer operators) and don’t bother with the error returns.

About the alignment, if your really out for performance you could have a member float *m_matrix that points to an aligned matrix. The reason for this is to allow use of SSE functions that require the data to be aligned. For example you could use _mm_stream_ps in your assignment operator if the memory were aligned. If you do decided to align the matrix don’t use any heap allocations just use an array.
yadango
yadango
another tip: add a copy constructor and assignment operator.
The Parrot
The Parrot
Something I noticed is that you compare the rows and the column count; you could have one member int function that stores both numbers, saving one check. If you say you never go beyond a 4×4 matrix you could just store the matrix number with (rows + (cols << 3)) and compare that.

But that's seriously minor :)
--Jeroen Stout - WebsiteCreator of-Divided
Zahlman
Zahlman
- malloc()? In MY C++?

- When you resize a matrix with initElements(), you leak the old memory. When a matrix falls out of scope, you leak its memory. When you assign or copy-construct a matrix, bad things happen (the memory is aliased).

- Notwithstanding all the things you need to do to make the dynamic allocation work properly, it's not even really called for. When does a matrix ever *change* size once instantiated? All we want is different size matrices; we can do that with templating, and as a free bonus, get some compile-time checking on our mathematical operations.

- Speaking of mathematical operations: In C++, we can overload operators, and it is normal to do so when there is an obvious and unambiguous meaning for the operator. This includes subscripting (operator[]). Anyway, having a special get() for the (0,0) element is highly unlikely to be useful.

- The error handling is seriously annoying. Using assertions or throwing exceptions would make more sense. I use assertions because it's assumed to be programmer error rather than user error if, say, an array subscript is wrong. That might not be true depending on your application and how much control is given to the user.

- The special handling for "vector matrices" is suspect: just transpose one of them and multiply the row and column normally.

template <typename T, size_t size>class Array {  T data[size];  public:  Array() { clear(); } // not sure this is needed  T& operator[](size_t idx) {    std::assert(idx >= 0 && idx < size);    return data[idx];  }  const T& operator[] const (size_t idx) {    std::assert(idx >= 0 && idx < size);    return data[idx];  }  void clear() {    std::fill(data, data + size, T());  }};template <size_t R, size_t C>class Matrix : public Array<Array<double, C>, R> {  public:  // clear() and subscripting operations are inherited, and  // subscripting twice yields a float&, as expected.  // retrieve information about the matrix  // These are probably useless, because dimensions are going to get checked at  // compile-time, but...  size_t rows() { return R; }  size_t columns() { return C; }  // We'll use a free function to create identity matrices.  // We can reset a matrix by assigning to it, with the default operator=.  // matrix math functions  Matrix& operator+= (const Matrix<R, C>& other) {    for (int i = 0; i < R; ++i) {      for (int j = 0; j < C; ++j) {        (*this)[j] += other[j];      }    }    return *this;  }  Matrix& operator-= (const Matrix<R, C>& other) {    for (int i = 0; i < R; ++i) {      for (int j = 0; j < C; ++j) {        (*this)[j] -= other[j];      }    }    return *this;  }  Matrix& operator*= (double scale) {    for (int i = 0; i < R; ++i) {      for (int j = 0; j < C; ++j) {        (*this)[j] *= scale;      }    }    return *this;  }  Matrix<C, R> transpose() {    Matrix<C, R> result;    for (int i = 0; i < R; ++i) {      for (int j = 0; j < C; ++j) {        result[j] = (*this)[j];      }    }    return result;  }};template <size_t R, size_t C>Matrix<R, C> operator+ (const Matrix<R, C>& lhs, const Matrix<R, C>& rhs) {  return Matrix(lhs) += rhs;}template <size_t R, size_t C>Matrix<R, C> operator- (const Matrix<R, C>& lhs, const Matrix<R, C>& rhs) {  return Matrix(lhs) -= rhs;}template <size_t R, size_t C>Matrix<R, C> operator* (const Matrix<R, C>& m, double scale) {  return Matrix(m) *= scale;}template <size_t R1, size_t C1R2, size_t C2>Matrix<R1, C2> operator* (const Matrix<R1, C1R2>& lhs,                          const Matrix<C1R2, C2>& rhs) {  Matrix<R1, C2> result;  for (int i = 0; i < R1; ++i) {    for (int j = 0; j < C2; ++j) {      for (int k = 0; k < C1R2; ++k) {        // Um, I think I got this right o_o        result[j] += lhs[k][j] * rhs[k];      }    }  }  return result;}template <size_t X>Matrix<X, X> identity() {  Matrix<X, X> result;  for (int i = 0; i < X; ++i) {    result = 1; // other values were set via constructors.  }  return result;}


[Edited by - Zahlman on July 10, 2006 2:17:28 AM]
SiCrane
SiCrane
Everything else aside, I don't think this function does what you think it does:
Matrix::Matrix() {	Matrix(1, 1);}
Lode
Lode
Quote:
Original post by SiCrane
Everything else aside, I don't think this function does what you think it does:
Matrix::Matrix() {	Matrix(1, 1);}


He's right, I've found this out the hard way: if you call the constructor of this class, it'll make a new one, it won't do anything to the current one. You cannot call one constructor from another one, you have to write every contructor on it's own or use an init function to avoid code duplication.
Homepage | LodePNG
Sharlin
Sharlin
Quote:
Original post by Lode
You cannot call one constructor from another one, you have to write every contructor on it's own or use an init function to avoid code duplication.


You can call a base class's constructor from an initializer list, though (this is the only way to do it, too). C++0x, when ready, will support delegating constructors.

tendifo
tendifo
Yea, I definitely see some problems with my code. I like the code posted by Zahlman. I might implement something like that.

Thanks for the input.
tendifo
tendifo
Ok so I reworked my code. I'm sticking with a static 3x3 matrix. That's all I'm ever going to need to be honest and it can handle both matrices and vectors so it's good.

Tell me what you think.

math.h

class Matrix {public:	Matrix();	double mget(int r, int c);	void mset(int r, int c, double v);	Matrix operator+(Matrix &rhs);	Matrix operator-(Matrix &rhs);	Matrix operator*(Matrix &rhs);	Matrix operator*(double s);	void clear();	void identity();	Matrix transpose();private:	double elements[3][3];};


math.cpp

#include "math.h"Matrix::Matrix() {	clear();}double Matrix::mget(int r, int c) {	return elements[r][c];}void Matrix::mset(int r, int c, double v) {	elements[r][c] = v;}Matrix Matrix::operator+(Matrix &rhs) {	Matrix temp;	for (int i = 0; i < 3; i += 1) {		for (int j = 0; j < 3; j += 1) {			temp.mset(i, j, mget(i, j) + rhs.mget(i, j));		}	}	return temp;}Matrix Matrix::operator-(Matrix &rhs) {	Matrix temp;	for (int i = 0; i < 3; i += 1) {		for (int j = 0; j < 3; j += 1) {			temp.mset(i, j, mget(i, j) - rhs.mget(i, j));		}	}	return temp;}Matrix Matrix::operator*(Matrix &rhs) {	Matrix temp;	double sum = 0;	for (int i = 0; i < 3; i += 1) {		for (int j = 0; j < 3; j += 1) {			sum = 0;			for (int k = 0; k < 3; k += 1) {				sum += mget(i, k) * rhs.mget(k, j);			}			temp.mset(i, j, sum);		}	}	return temp;}Matrix Matrix::operator*(double s) {	Matrix temp;	for (int i = 0; i < 3; i += 1) {		for (int j = 0; j < 3; j += 1) {			temp.mset(i, j, mget(i, j) * s);		}	}	return temp;}void Matrix::clear() {	for (int i = 0; i < 3; i += 1) {		for (int j = 0; j < 3; j += 1) {			elements[j] = 0;		}	}}void Matrix::identity() {	clear();	mset(0, 0, 1);	mset(1, 1, 1);	mset(2, 2, 1);}Matrix Matrix::transpose() {	Matrix temp;	for (int i = 0; i < 3; i += 1) {		for (int j = 0; j < 3; j += 1) {			temp.mset(i, j, mget(j, i));		}	}	return temp;}
deathkrush
deathkrush
How about overloading the [] operator?

class Matrix {public:    double* operator[](unsigned int i)    {        assert(i<rows*columns);        return elements + i*columns;    }};


Then, you can access your matrix like this:
Matrix m;
m[0][1] = 0;

But that kind of blows away getters and setters.
deathkrushPS3/Xbox360 Graphics Programmer, Mass Media.Completed Projects: Stuntman Ignition (PS3), Saints Row 2 (PS3), Darksiders(PS3, 360)
Red Ant
Red Ant
Also think about const-correctness. Member functions which do not modify the instance they're invoked on should be const. For instance your

double getElement(int r, int c);

method could easily be

double getElement(int r, int c) const;


What's more, avoid passing matrixes by value. Pass them by const reference instead. This saves you some unnecessary copying.
angry
angry
I think you should make the entire class inline. And instead of a 2-dimensional array I think you should have just a 1-dimensional array, and then skip the for's and use memcpy and memset instead.
Zahlman
Zahlman
Quote:
Original post by angry
I think you should make the entire class inline. And instead of a 2-dimensional array I think you should have just a 1-dimensional array, and then skip the for's and use memcpy and memset instead.


Please, use std::copy and std::fill, rather than memset and memcpy. :s Also, a 2D array of the sort in use here (both dimensions static) doesn't incur any overhead, so there's really no reason to switch to a 1-d array.

Note that implementing the matrix as an array of arrays is much nicer when it comes to making subscripting work. It also lets you replace a row at a time of the matrix. And it shouldn't be any slower, as long as everything is inline, except when you index with constants for both indices AND the compile-time complexity (nothing is really being added in the run-time data) confuses the compiler enough that it can't identify and do the constant folding.
tendifo
tendifo
Alright so just to make sure I'm getting the code right.. to implement const I should use:

double mget(int r, int c) const;


And to return a reference instead of value I should do:

Matrix &operator+(const Matrix &rhs);


Is that the correct coding?

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.