Matrix Class, =, =*, * operators

[SLH]

I've got a matrix class, that manages an 1D float array (using indexing to make it act as a 2D array).

I'm having trouble implementing various operators. The problem i get is that it's trying to delete the array twice.

I've stepped through the code, and i still don't see why it's calling the deconstructor when it does.

int main(void)
{
	CMatrix testMat(2,2,1.0f);
	CMatrix testMat2(2,2,1.0f);
	testMat = testMat * testMat2; //Destructor called after this line
	testMat.Print();
	return 0;
}

/////////////////////////
// Function prototypes in header file
/////////////////////////
	CMatrix operator*(const float Value);
	CMatrix operator*(const CMatrix &RHS); //Matrix Multiply
	CMatrix &operator*=(const CMatrix &RHS); //Matrix Multiply
	CMatrix &operator*=(const float Value);
	CMatrix &operator=(const CMatrix &RHS); //To duplicate a matrix
/////////////////////////
// End Function prototypes
/////////////////////////
CMatrix::CMatrix()
{
	m_pMatrix = NULL;
	m_pInvMatrix = NULL;
	m_pTranMatrix = NULL;
	m_pAdjMatrix = NULL;
}

CMatrix::CMatrix(int Rows, int Columns)
{
	m_pMatrix = new float[Rows*Columns];
	m_pInvMatrix = new float[Rows*Columns];
	m_pTranMatrix = new float[Rows*Columns];
	m_pAdjMatrix = new float[Rows*Columns];
	
	m_bInvIsDirty = m_bTranIsDirty = m_bDetIsDirty = m_bAdjIsDirty = true;
	
	m_iColumns = Columns;
	m_iRows = Rows;
}

CMatrix::CMatrix(int Rows, int Columns, float Value)
{
	m_pMatrix = new float[Rows*Columns];
	m_pInvMatrix = new float[Rows*Columns];
	m_pTranMatrix = new float[Rows*Columns];
	m_pAdjMatrix = new float[Rows*Columns];
	for(int i=0;i<Rows*Columns;i++)
		m_pMatrix[i] = Value;
	
	m_bInvIsDirty = m_bTranIsDirty = m_bDetIsDirty = m_bAdjIsDirty = true;
	
	m_iColumns = Columns;
	m_iRows = Rows;
}

CMatrix::CMatrix(int Rows,int Columns,float *Matrix)
{
	m_pMatrix = new float[Rows*Columns];
	m_pInvMatrix = new float[Rows*Columns];
	m_pTranMatrix = new float[Rows*Columns];
	m_pAdjMatrix = new float[Rows*Columns];
	for(int i=0;i<Rows*Columns;i++)
		m_pMatrix[i] = Matrix[i];
	
	m_bInvIsDirty = m_bTranIsDirty = m_bDetIsDirty = m_bAdjIsDirty = true;
	
	m_iColumns = Columns;
	m_iRows = Rows;
}

CMatrix::~CMatrix()
{
	delete[] m_pMatrix;
	delete[] m_pInvMatrix;
	delete[] m_pTranMatrix;
	delete[] m_pAdjMatrix;
}

CMatrix CMatrix::operator*(const float Value)
{
	CMatrix Result(m_iColumns,m_iRows,m_pMatrix);
	
	for(int i=0;i<m_iRows*m_iColumns;i++)
		m_pMatrix[i] *= Value;
	
	return Result;
}

CMatrix CMatrix::operator*(const CMatrix &RHS)
{
	CMatrix Result(m_iRows,RHS.m_iColumns);
	
	for(int i=0;i<m_iRows;i++)
		for(int j=0;j<RHS.m_iColumns;j++)
			for(int k=0;k<m_iColumns;k++)
				Result.m_pMatrix[i+j*RHS.m_iColumns] += m_pMatrix[i+k*m_iColumns] * RHS.m_pMatrix[k+j*m_iColumns];
	
	return Result;
}

CMatrix &CMatrix::operator*=(const CMatrix &RHS)
{
	assert(RHS.m_iColumns == m_iColumns);
	
	CMatrix tmp(m_iRows,RHS.m_iColumns);
	
	tmp = *this * RHS;
	
	memcpy(m_pMatrix,tmp.m_pMatrix,sizeof(float)*m_iColumns*m_iRows);

	m_bInvIsDirty = m_bTranIsDirty = m_bDetIsDirty = m_bAdjIsDirty = true;

	return *this;
}

CMatrix &CMatrix::operator*=(const float Value)
{
	for(int i=0;i<m_iRows*m_iColumns;i++)
	{
		m_pMatrix[i] *= Value;
	}
	return *this;
}

CMatrix &CMatrix::operator=(const CMatrix &RHS)
{
	//Copy the Columns and Rows of the matrix
	m_iColumns = RHS.m_iColumns;
	m_iRows = RHS.m_iRows;
	delete[] m_pMatrix;
	delete[] m_pInvMatrix;
	delete[] m_pTranMatrix;
	delete[] m_pAdjMatrix;
	
	m_pMatrix = new float[m_iRows*m_iColumns];
	m_pInvMatrix = new float[m_iRows*m_iColumns];
	m_pTranMatrix = new float[m_iRows*m_iColumns];
	m_pAdjMatrix = new float[m_iRows*m_iColumns];

	//Copy the matrix array
	memcpy(m_pMatrix,RHS.m_pMatrix,sizeof(float)*m_iColumns*m_iRows);
	//Copy the cached matrices dirty flag values, and the cached matrices themselves if appropriate
	m_bInvIsDirty = RHS.m_bInvIsDirty;
	if(!m_bInvIsDirty)
		memcpy(m_pInvMatrix,RHS.m_pInvMatrix,sizeof(float)*m_iColumns*m_iRows);
	
	m_bTranIsDirty = RHS.m_bTranIsDirty;
	if(!m_bTranIsDirty)
		memcpy(m_pTranMatrix,RHS.m_pTranMatrix,sizeof(float)*m_iColumns*m_iRows);
	
	m_bAdjIsDirty = RHS.m_bAdjIsDirty;
	if(!m_bAdjIsDirty)
		memcpy(m_pAdjMatrix,RHS.m_pAdjMatrix,sizeof(float)*m_iColumns*m_iRows);
	
	//Copy the cached determinent flag and value if appropriate
	m_bDetIsDirty = RHS.m_bDetIsDirty;
	if(!m_bDetIsDirty)
		m_fDet = RHS.m_fDet;
	
	return *this;
}

Any help would be greatly appreciated.

[CornedBee]

I can't tell you, without debugging (and for that I lack the complete source), why it's deleting twice. It might be something about a missing copy constructor, but I doubt it.

I can tell you a few other things, however.
1) Your implementation of *(Matrix) is broken. You add to uninitialized and therefore undefined values, yielding undefined values, of course. It does not test for multipliability either.
2) Your implementation of *=(Matrix) is broken. The debug check at the start of the function does not make mathematical sense. (You must test RHS.m_iRows against m_iColumns, and you should make it a check that exists in release builds as well.) However, curiously enough the test ensures that the second problem of the function does not occur (in debug builds): incomplete copying of the result matrix if the RHS matrix has more columns than the LHS. The third problem is one of inefficiency: you construct tmp with a row and column count, thus allocating memory. Then you assign to that matrix, which frees that memory and allocates its own.
3) Your implementation of *=(float) is broken: it does not invalidate the dependent matrices and values.
4) The way you're implementing the operators is the wrong way round, according to the usual techniques. You should implement * in terms of *=, not the other way round. * should be a free function, then.
5) You shouldn't immediately allocate memory for the dependent matrices. Many matrices cannot even have inverts, so what's the use in having a place to store it? This way, you can get rid of the dirty flags, by simply setting the pointers to NULL. Also, is the transposed matrix useful to store? It is so easy to get.
6) Your indexing is broken. You use j as the row index, which makes the multiplication forums wrong.

#include <limits>
#include <stdexcept>

// Use this if no determinant was calculated.
const float NaN = std::numeric_limits<float>::quiet_NaN();

Matrix & Matrix::operator *=(float f)
{
  for(int i = 0; i < m_iRows*m_iColumns; ++i) {
    m_pMatrix[i] *= f;
  }
  delete[] m_pInvMatrix; m_pInvMatrix = NULL;
  delete[] m_pTransMatrix; m_pTransMatrix = NULL;
  delete[] m_pAdjMatrix; m_pAdjMatrix = NULL;
  if(m_fDet != NaN) {
    m_fDet *= pow(f, m_iColumns); // Should use an optimized version for integers powers here.
  }
  return *this;
}

const Matrix operator *(const Matrix &lhs, float rhs)
{
  return Matrix(lhs) *= rhs;
}

const Matrix operator *(float lhs, const Matrix &rhs)
{
  // Multiplication by scalar is commutative.
  return Matrix(rhs) *= lhs;
}

Matrix & Matrix::operator *=(const Matrix &rhs)
{
  if(m_iColumns != rhs.m_iRows) {
    throw std::invalid_argument("Cannot multiply matrices.");
  }
  float *target = new float[m_iRows * rhs.m_iColumns];
  memset(target, 0, sizeof(float) * m_iRows * rhs.m_iColumns);
  for(int i = 0; i < m_iRows; ++i) {
    for(int j = 0; j < rhs.m_iColumns; ++j) {
      for(int k = 0; k < m_iColumns; ++k) {
        target[i*rhs.m_iColumns + j] += m_pMatrix[i*m_iColumns + k] *
            rhs.m_pMatrix[k*rhs.m_iColumns + j];
      }
    }
  }
  delete[] m_pInvMatrix; m_pInvMatrix = NULL;
  delete[] m_pTransMatrix; m_pTransMatrix = NULL;
  delete[] m_pAdjMatrix; m_pAdjMatrix = NULL;
  if(m_fDet != NaN && rhs.m_fDet != NaN) {
    m_fDet *= rhs.m_fDet;
  } else {
    m_fDet = NaN;
  }
  delete[] m_pMatrix;
  m_pMatrix = target;
}

const Matrix operator *(const Matrix &lhs, const Matrix &rhs)
{
  return Matrix(lhs) *= rhs;
}

Matrix::Matrix(const Matrix &rhs)
  : m_iRows(rhs.m_iRows),
  m_iColumns(rhs.m_iColumns),
  m_fDet(rhs.m_fDet)
{
  m_pMatrix = new float[m_iRows * m_iColumns];
  memcpy(m_pMatrix, rhs.m_pMatrix, sizeof(float) * m_iRows * m_iColumns);

  if(rhs.m_pInvMatrix) {
    m_pInvMatrix = new float[m_iRows * m_iColumns];
    memcpy(m_pInvMatrix, rhs.m_pInvMatrix, sizeof(float) * m_iRows * m_iColumns);
  } else {
    m_pInvMatrix = NULL;
  }

  if(rhs.m_pTranMatrix) {
    m_pTranMatrix = new float[m_iRows * m_iColumns];
    memcpy(m_pTranMatrix, rhs.m_pTranMatrix, sizeof(float) * m_iRows * m_iColumns);
  } else {
    m_pTranMatrix = NULL;
  }

  if(rhs.m_pAdjMatrix) {
    m_pAdjMatrix = new float[m_iRows * m_iColumns];
    memcpy(m_pAdjMatrix, rhs.m_pAdjMatrix, sizeof(float) * m_iRows * m_iColumns);
  } else {
    m_pAdjMatrix = NULL;
  }
}

Matrix & Matrix::operator =(const Matrix &rhs)
{
  bool needRealloc = m_iRows * m_iColumns != rhs.m_iRows * rhs.m_iColumns;

  m_iRows = rhs.m_iRows;
  m_iColumns = rhs.m_iColumns;

  if(needRealloc) {
    delete[] m_pMatrix;
    m_pMatrix = new float[m_iRows * m_iColumns];
  }
  memcpy(m_pMatrix, rhs.m_pMatrix, sizeof(float) * m_iRows * m_iColumns);

  if(rhs.m_pInvMatrix) {
    if(needRealloc) {
      delete[] m_pInvMatrix;
      m_pInvMatrix = new float[m_iRows * m_iColumns];
    }
    memcpy(m_pInvMatrix, rhs.m_pInvMatrix, sizeof(float) * m_iRows * m_iColumns);
  } else {
    delete[] m_pInvMatrix; m_pInvMatrix = NULL;
  }

  if(rhs.m_pTranMatrix) {
    if(needRealloc) {
      delete[] m_pTranMatrix;
      m_pTranMatrix = new float[m_iRows * m_iColumns];
    }
    memcpy(m_pTranMatrix, rhs.m_pTranMatrix, sizeof(float) * m_iRows * m_iColumns);
  } else {
    delete[] m_pTranMatrix; m_pTranMatrix = NULL;
  }

  if(rhs.m_pAdjMatrix) {
    if(needRealloc) {
      delete[] m_pAdjMatrix;
      m_pAdjMatrix = new float[m_iRows * m_iColumns];
    }
    memcpy(m_pAdjMatrix, rhs.m_pAdjMatrix, sizeof(float) * m_iRows * m_iColumns);
  } else {
    delete[] m_pAdjMatrix; m_pAdjMatrix = NULL;
  }

  m_fDet = rhs.m_fDet;
}

Note that this is not exception safe, and will not be until you use smart pointers for the arrays. (On a good compiler, a boost::scoped_array has no overhead at all.)

[SLH]

Thanks a lot for your help, i think i understand all your points.

I've since decided that it might be better to use templates, but i've tried to apply many of your points to my new template class.
If you can spot any errors, or can see where i've done things in a really bad way i'd appreciate any pointers.

template <int m_iRows, int m_iColumns, int m_iLength>
class CMatrix
{
public:
	//Class functions
	CMatrix()
	{
		m_bInvIsDirty = m_bTranIsDirty = m_bAdjIsDirty = true;
	}
	CMatrix(float Value)
	{
		for(int i=0;i<m_iLength;i++)
			m_pMatrix[i] = Value;
		m_fDet = NaN;
	}
	~CMatrix(){}
	//Operator functions
	CMatrix operator*(const float Value)
	{
		CMatrix<m_iRows,m_iColumns,m_iRows*m_iColumns> Result(m_pMatrix);
		for(int i=0;i<m_iRows*m_iColumns;i++)
			m_pMatrix[i] *= Value;
		
		return Result;
	}
	CMatrix operator*(const CMatrix &RHS)
	{
		return CMatrix(RHS) *= *this;
	}
	CMatrix &operator*=(const CMatrix &RHS)
	{
		float target[m_iLength];
		memset(target, 0, sizeof(float) * m_iLength);

		for(int i = 0; i < m_iRows; ++i)
		{
			for(int j = 0; j < m_iColumns; ++j)
			{
				for(int k = 0; k < m_iColumns; ++k)
				{
					target[i*m_iColumns + j] += m_pMatrix[i*m_iColumns + k] *
					RHS.m_pMatrix[k*m_iColumns + j];
				}
			}
		}
		m_bInvIsDirty = m_bTranIsDirty = m_bAdjIsDirty = true;
		memcpy(m_pMatrix,target,sizeof(float)*m_iLength);
		return *this;
	}
	CMatrix &operator*=(const float Value)
	{
		for(int i=0;i<m_iLength;i++)
		{
			m_pMatrix[i] *= Value;
		}
		return *this;
	}
	CMatrix &operator=(const CMatrix &RHS)
	{
		//Copy the matrix array
		memcpy(m_pMatrix,RHS.m_pMatrix,sizeof(float)*m_iLength);
		//Copy the cached matrices dirty flag values, and the cached matrices themselves if appropriate
		m_bInvIsDirty = RHS.m_bInvIsDirty;
		if(!m_bInvIsDirty)
			memcpy(m_pInvMatrix,RHS.m_pInvMatrix,sizeof(float)*m_iLength);
		
		m_bTranIsDirty = RHS.m_bTranIsDirty;
		if(!m_bTranIsDirty)
			memcpy(m_pTranMatrix,RHS.m_pTranMatrix,sizeof(float)*m_iLength);
		
		m_bAdjIsDirty = RHS.m_bAdjIsDirty;
		if(!m_bAdjIsDirty)
			memcpy(m_pAdjMatrix,RHS.m_pAdjMatrix,sizeof(float)*m_iLength);
		
		m_fDet = RHS.m_fDet;
		
		return *this;
	}
	__forceinline float &operator()(const int Column,const int Row)
	{
		return m_pMatrix[Row+Column*m_iColumns];
	}
	//Matrix operations
	//CMatrix Inverse(void);
	//CMatrix Transpose(void);
	//CMatrix Adjoint(void);
	//float Determinant(void);

	int GetColumns(void){return m_iColumns;}
	int GetRows(void){return m_iRows;}

	void Print(void)
	{
		for(int i=0;i<m_iRows;i++)
		{
			for(int j=0;j<m_iColumns;j++)
			{
				printf("%f ",m_pMatrix[i+j*m_iColumns]);
			}
			cout << endl;
		}
	}
protected:
	//An array representing the matrix
	float m_pMatrix[m_iLength];
	//Cached matrix calculations to avoid re-calculation
	float m_pInvMatrix[m_iLength];
	float m_pTranMatrix[m_iLength];
	float m_pAdjMatrix[m_iLength];
	float m_fDet;
	//Booleans to determine whether a matrix calculation needs updating
	bool m_bInvIsDirty;
	bool m_bTranIsDirty;
	bool m_bAdjIsDirty;
};

Note: My class makes the assumption that all matrix manipulations will be between matracies of the same dimensions (as they will be in my prog).

[CornedBee]

Not got much time right now, just a few things:
1) I dislike the iLength template parameter. The length is rows*cols, the additional parameter is only a source of errors.
2) Template parameters have the status of real constants. This means you can use real 2-d arrays now. float data[rows][cols];
3) If all your matrices in your app are the same size and you ever hope to multiply them, then a logical consequence is that they're all square. Make use of this - there's a lot of assumptions you can make with square matrices. Not least, you can get rid of another template parameter.
4) Now that all your arrays are direct members, they directly increase the size of the class. Thus, a 3x3 matrix now requires no less than 151 bytes, and perhaps more. A 4x4 matrix (the most common in 3d graphics) requires 263 bytes. That's quite large, especially if the cached values are not always used and thus wasted.
5) I still don't believe caching the transposed matrix brings you any performance advantage. It would be better to write a matrix-like proxy object that simply flips the indices in its accesses to a real matrix.

[SLH]

Thanks for taking the time to give me lots of constructive criticism, it's greatly appreciated.

1 & 2) I included the Length parameter so that i wouldn't have to keep calculating the row * column when iterating through my array. Since i can now use a 2D array anyway, i won't need that parameter any more. Perhaps i should have made it a member variable, that was initialised in the constructor?
3) Thanks for that, should have seen that myself really!
4) For now i'm gonna leave it as it is, but i'll certainly consider changing it if memory usage becomes an issue.
5) I've seen proxy object's mentioned in other threads (relating to allowing the programmer to do this: My2DArrayClass[i][j]). I'm not too sure how i'd go about implementing one though, any help would be appreciated.

[CornedBee]

At the beginning of this discussion I dug up my old matrix class and was rather dismayed to find just how bad it really was. (Relic from my early days.) I think I'll reimplement it, and then post a thoroughly documented version here.

[SLH]

Cool, thanks!

[CornedBee]

Different idea. I'll be documenting my progress here, inviting comments at every step.

First, some basic considerations.

We want to create a class for square matrices. It should be as generic as possible, and yet speed is a very important consideration. And we want to write as little code as possible, yet stay portable across compilers.

For the first, we will make it a template. It will have two parameters: the underlying numeric type and the size of the matrix.

For the second, matrix operations lend themselves quite excellently to SIMD instructions. How do we use those without getting bogged down in low-level assembly, though? Can we trust the compiler to unroll and parallellize our loops?
In part, we can, but C++ has something even better: std::valarray is a class that is intended for the very purpose of high-speed mathematical calculations on many numbers at once. It is a 1-dimensional container that offers various operations on all its numbers. We will trust the library implementors to optimize this class. As such, our matrix will at its core consist of

So, first some boilerplate code.

////////////////////////////////////////////////////////////////////////////
// (c) 2005 Sebastian Redl
// Available under the GNU General Public License. For details, go here:
// http://www.gnu.org/licenses/gpl.html

// A generic square matrix class.

#ifndef _MATRIX_HPP_SR_H_
#define _MATRIX_HPP_SR_H_

template <unsigned int SIDE, typename T = float>
class matrix
{
};

#endif

[CornedBee]

The next step is to define the interface that our class should have. We'll start with a very basic class and extend it later.

We want:
Default-construction to all-0 (because that's what valarray does).
Copy-construction.
Fill-construction from a single value.
Fill-construction from an input iterator. (Sufficient elements are assumed.)
A static method to return the identity matrix.

Destruction.

Matrix addition.
Matrix multiplication.
Multiplication with a scalar.

Indexed and bounds-checked element access using [][].

That should be enough for the start.

The interface looks like this. We will think later on what unspecified looks like.

template <unsigned int SIDE, typename T = float>
class matrix
{
public:
	typedef matrix<SIDE, T>		my_type;
	typedef T					value_type;
	typedef T &					reference;
	typedef const T &			const_reference;

	typedef unspecified			index_proxy;

	matrix();

	matrix(const my_type &rhs);

	explicit matrix(const_reference rhs);

	template <class InputIterator>
	explicit matrix(InputIterator in);

	static my_type identity();

	~matrix();


	my_type & operator +=(const my_type &rhs);

	my_type & operator *=(const my_type &rhs);

	my_type & operator *=(const T &rhs);

	index_proxy operator[](size_type row);

};

[SLH]

How's this little project going? I'm looking foward to how it turns out!
No probs if you don't have enough free time or anything.

[CornedBee]

It's going slowly - I have very little time for it indeed, especially as I don't really need it right now.