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"; }