C++ Template Class for Linear Algebra Matrix

It seems that many projects come upon a need to perform some linear algebra maths. However, the a large number of the available libraries usually are huge with many unnecessary or system dependent functionalities or …they are not free.

The following template class provides some basic linear algebra matrix operations such as +,-,* as well as find the determinant, the covariance matrix etc…

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

enum Orientation
{
	ROWS,
	COLUMNS
};

template <typename T> class Matrix
{
	private:

	T* data;

	unsigned int dimensions[ 2 ];

	double getDeterminant ( Matrix &X );

	double getCovariance ( T *, T *, size_t );

	protected:

	public:

	// Destructor

	~Matrix ( );

	// Constructors

	Matrix ( );

	Matrix ( int );

	Matrix ( int, int );

	template<typename... Args> Matrix ( Args... args );

	// Operators overloading

	Matrix<T>& operator=( Matrix<T>& );

	Matrix<T>& operator=( std::initializer_list<T> );

	Matrix<T>& operator=( std::initializer_list<std::initializer_list<T>> );

	Matrix<T> operator+( const Matrix<T>& );

	Matrix<T> operator-( const Matrix<T>& );

	Matrix<T> operator*( const Matrix<T>& );

	// Other Methods

	Matrix<T> getCovariance ( );

	Matrix<T> getTransposed ( void );

	Matrix<T> getMean ( Orientation );

	void transpose ( void );

	double getDeterminant ( );

	void clear ( void );

	void print ( void );
};

Constructors

template<typename T> inline Matrix<T>::Matrix ( )
{
	clear ( );
}
template<typename T> inline Matrix<T>::Matrix ( int size1 )
{
	clear ( );

	dimensions[ 0 ] = size1;

	dimensions[ 1 ] = 1;

	data = ( T* ) calloc ( dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ), sizeof ( T ) );	
}

template<typename T> inline Matrix<T>::Matrix ( int size1 , int size2 )
{
	clear ( );

	dimensions[ 0 ] = size1;

	dimensions[ 1 ] = size2;

	data = ( T* ) calloc ( dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ), sizeof ( T ) );
}

template<typename T> template<typename ...Args> inline Matrix<T>::Matrix ( Args ...args )
{
	static_assert( sizeof...( args ) <= 2 && sizeof...( args ) >= 1 || sizeof...( args ) == 0, "Error: Supported Matrix Dimensions : 1, 2" );

	unsigned int res[ sizeof...( args ) ] = { args... };

	dimensions[ 0 ] = res[ 0 ];

	dimensions[ 1 ] = sizeof...( args ) == 2 ? res[ 1 ] : 1;

	data = ( T* ) calloc ( dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ), sizeof ( T ) );
}

template<typename T> inline Matrix<T>::~Matrix()
{
	//	clear();
}

Operators overloading

template<typename T> inline Matrix<T>& Matrix<T>::operator=( Matrix<T>& matrix )
{
	clear ( );

	dimensions[ 0 ] = matrix.dimensions[ 0 ];

	dimensions[ 1 ] = matrix.dimensions[ 1 ];

	data = ( T* ) calloc ( dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ), sizeof ( T ) );

	memcpy ( data, matrix.data, dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ) * sizeof ( T ) );

	return *this;
}

template<typename T> inline Matrix<T>& Matrix<T>::operator=( std::initializer_list<T> other )
{
	clear ( );

	dimensions[ 0 ] = 1;

	dimensions[ 1 ] = other.size ( );

	data = ( T* ) calloc ( dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ), sizeof ( T ) );

	memcpy ( data, other.begin ( ), other.size ( ) * sizeof ( T ) );

	return *this;
}

template<typename T> inline Matrix<T>& Matrix<T>::operator=( std::initializer_list<std::initializer_list<T>> other )
{
	clear ( );

	dimensions[ 0 ] = other.size ( );

	dimensions[ 1 ] = other.begin ( )->size ( );

	data = ( T* ) calloc ( dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ), sizeof ( T ) );

	T * t_data = data - other.begin ( )->size ( );

	for ( auto x : other )

		memcpy ( t_data += x.size ( ), x.begin ( ), x.size ( ) * sizeof ( T ) );

	return *this;
}

template<typename T> inline Matrix<T> Matrix<T>::operator+( const Matrix<T>& matrix )
{
	if ( this->dimensions[ 0 ] != matrix.dimensions[ 0 ] || this->dimensions[ 1 ] != matrix.dimensions[ 1 ] )

		exit ( -1 );

	Matrix<T> returned_matrix ( dimensions[ 0 ], dimensions[ 1 ] );

	for ( int i = 0; i < dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ); i++ )

		returned_matrix.data[ i ] = data[ i ] + matrix.data[ i ];

	return returned_matrix;
}

template<typename T> inline Matrix<T> Matrix<T>::operator-( const Matrix<T>& matrix )
{
	if ( this->dimensions[ 0 ] != matrix.dimensions[ 0 ] || this->dimensions[ 1 ] != matrix.dimensions[ 1 ] )

		exit ( -1 );

	Matrix<T> returned_matrix ( dimensions[ 0 ], dimensions[ 1 ] );

	for ( int i = 0; i < dimensions[ 0 ] * ( dimensions[ 1 ] == 0 ? 1 : dimensions[ 1 ] ); i++ )

		returned_matrix.data[ i ] = data[ i ] - matrix.data[ i ];

	return returned_matrix;
}

template<typename T> inline Matrix<T> Matrix<T>::operator*( const Matrix<T>& matrix )
{
	if ( this->dimensions[ 1 ] != matrix.dimensions[ 0 ] )

		exit ( -1 );

	Matrix<T> returned_matrix ( dimensions[ 0 ], matrix.dimensions[ 1 ] );

	for ( int i = 0; i < dimensions[ 0 ]; i++ ) // First Matrix Rows

		for ( int j = 0; j < matrix.dimensions[ 1 ]; j++ ) // Second Matrix Cols

			for ( int k = 0; k < dimensions[ 1 ]; k++ ) // First Matrix Cols

				returned_matrix.data[ i * matrix.dimensions[ 1 ] + j ]

				+=

				data[ i * ( dimensions[ 1 ] ) + k ]

				*

				matrix.data[ k * ( matrix.dimensions[ 1 ] ) + j ];

	return returned_matrix;
}

Other Methods

template<typename T> inline double Matrix<T>::getDeterminant ( )
{
	return getDeterminant ( *this ); 
}

template<typename T> inline double Matrix<T>::getDeterminant ( Matrix & X )
{
	if ( dimensions[ 0 ] != dimensions[ 1 ] || dimensions[ 0 ] < 2 )

		exit ( -1 );

	if ( dimensions[ 0 ] > 2 )
	{
		double result = 0;

		Matrix<T> Y ( dimensions[ 0 ] - 1, dimensions[ 0 ] - 1 );

		for ( int i = 0; i < dimensions[ 0 ]; i++ )
		{
			for ( int s_i = 0; s_i < Y.dimensions[ 0 ]; s_i++ )

				for ( int s_j = 0; s_j < Y.dimensions[ 0 ]; s_j++ )

					Y.data[ s_i * Y.dimensions[ 1 ] + s_j ] = X.data[ ( ( s_i + 1 ) * dimensions[ 1 ] ) + s_j >= i ? s_j + 1 : s_j ];

					result += ( i % 2 == 0 ? 1 : -1 ) * X.data[i] * getDeterminant ( Y );
		}

		return result;
	}

	return data[ 0 ] * data[ 3 ] - data[ 1 ] * data[ 2 ];
}

template<typename T> inline double Matrix<T>::getCovariance ( T * vector1, T * vector2, size_t size )
{
	double x_mean = 0, y_mean = 0, z_cov = 0;

	for ( int i = 0; i < size; x_mean += vector1[ i++ ] );

	x_mean /= size;

	for ( int i = 0; i < size; y_mean += vector2[ i++ ] );

	y_mean /= size;

	for ( int i = 0; i < size; i++ )

		z_cov += ( vector1[ i ] - x_mean ) * ( vector2[ i ] - y_mean );

	return z_cov / ( size - 1 );
}

template<typename T> inline Matrix<T> Matrix<T>::getCovariance ( )
{
	Matrix<T> covariance_matrix ( dimensions[ 0 ], dimensions[ 0 ] );

	this->transpose ( );

	for ( int i = 0; i < covariance_matrix.dimensions[ 0 ]; i++ )

		for ( int j = 0; j < covariance_matrix.dimensions[ 0 ]; j++ )

			covariance_matrix.data[ j * covariance_matrix.dimensions[ 0 ] + i ] = getCovariance ( &data[ i*dimensions[ 0 ] ], &data[ j*dimensions[ 0 ] ], dimensions[ 1 ] );

	this->transpose ( );

	return covariance_matrix;
}

template<typename T> inline Matrix<T> Matrix<T>::getTransposed ( )
{
	Matrix<T> returned_matrix ( dimensions[ 1 ], dimensions[ 0 ] );

	for ( int i = 0; i < dimensions[ 0 ]; i++ )

		for ( int j = 0; j < dimensions[ 1 ]; j++ )

			returned_matrix.data[ j*dimensions[ 0 ] + i ] = data[ i*dimensions[ 1 ] + j ];

	return returned_matrix;
}

template<typename T> inline void Matrix<T>::transpose()
{
	Matrix<T> returned_matrix(dimensions[1], dimensions[0]);

	for (int i = 0; i < dimensions[0]; i++)

		for (int j = 0; j < dimensions[1]; j++)

			returned_matrix.data[j*dimensions[0] + i] = data[i*dimensions[1] + j];

	clear();

	*this = returned_matrix;

	returned_matrix.clear();
}

template<typename T> inline Matrix<T> Matrix<T>::getMean ( Orientation orientation )
{
	Matrix<T> returned_matrix ( orientation == ROWS ? dimensions[ 0 ] : 1, orientation == ROWS ? 1 : dimensions[ 1 ] );

	for ( unsigned int i = 0; i < dimensions[ orientation == ROWS ? 0 : 1 ]; i++ )
	{
		for ( unsigned int j = 0; j < dimensions[ orientation == ROWS ? 1 : 0 ]; j++ )

			returned_matrix.data[ i ] += this->data[ ( ( orientation == ROWS ? i : j ) * dimensions[ 1 ] ) + ( orientation == ROWS ? j : i ) ];

		returned_matrix.data[ i ] /= dimensions[ orientation == ROWS ? 1 : 0 ];
	}

	return returned_matrix;
}

template<typename T> inline void Matrix<T>::clear()
{
	dimensions[0] = 0;

	dimensions[1] = 0;

	if (data != NULL)

		free(data);

	data = NULL;
}

template<typename T> inline void Matrix<T>::print(void)
{
	for (unsigned int i = 0; i < dimensions[0]; i++)
	{
		for (unsigned int j = 0; j < dimensions[1]; j++)

			cout << "| " << data[i*dimensions[1] + j] << " ";

		cout << "|\n";
	}
	cout << "\n";
}