fixed memory leaks in stim::matrix

David Mayerich
1 parent c495d065
Showing 1 changed file with 150 additions and 34 deletions Show diff stats
stim/math/matrix.h
@@ -17,42 +17,79 @@ class matrix {
 	size_t R;							//number of rows
 	size_t C;							//number of colums
-	void init(size_t rows, size_t cols){
+	/*void init(size_t rows, size_t cols){
 		R = rows;
 		C = cols;
-		M = (T*) malloc (R * C * sizeof(T));	//allocate space for the matrix
+		if (R == 0 || C == 0) M = NULL;
+		else
+			M = (T*)malloc(R * C * sizeof(T));	//allocate space for the matrix
+	}*/
+
+	T get(const size_t row, const size_t col) const {
+		if (row >= R || col >= C) {
+			std::cout << "ERROR: row or column out of range." << std::endl;
+			exit(1);
+		}
+		return M[col * R + row];
 	}
 	T& at(size_t row, size_t col){
+		if (row >= R || col >= C) {
+			std::cout << "ERROR: row or column out of range." << std::endl;
+			exit(1);
+		}
 		return M[col * R + row];
 	}
 public:
-	
+	matrix() {
+		R = 0;
+		C = 0;
+		M = NULL;
+	}
+
 	matrix(size_t rows, size_t cols) {
-		init(rows, cols);									//initialize memory
+		R = rows;
+		C = cols;
+		if (R * C == 0) 
+			M = NULL;
+		else
+			M = new T[R * C];
 	}
 	matrix(size_t rows, size_t cols, T* data) {
-		init(rows, cols);
+		R = rows;
+		C = cols;
+		if (R * C == 0)
+			M = NULL;
+		else
+			M = new T[R * C];
 		memcpy(M, data, R * C * sizeof(T));
 	}
 	matrix(const matrix<T>& cpy){
-		init(cpy.R, cpy.C);
-		memcpy(M, cpy.M, R * C * sizeof(T));
+		
+		if (cpy.R * cpy.C == 0)
+			M = NULL;
+		else
+			M = new T[cpy.R * cpy.C];
+		memcpy(M, cpy.M, cpy.R * cpy.C * sizeof(T));
+
+		R = cpy.R;
+		C = cpy.C;
 	}
 	~matrix() {
-		R = C = 0;
 		if(M) free(M);
+		M = NULL;
+		R = C = 0;
 	}
-	size_t rows(){
+	size_t rows() const {
 		return R;
 	}
-	size_t cols(){
+	size_t cols() const {
 		return C;
 	}
@@ -60,10 +97,8 @@ public:
 		return at(row, col);
 	}
-	matrix<T> operator=(T rhs) {
-		if (&rhs == this)
-			return *this;
-		init(R, C);
+	matrix<T>& operator=(const T rhs) {
+		//init(R, C);
 		size_t N = R * C;
 		for(size_t n=0; n<N; n++)
 			M[n] = rhs;
@@ -71,14 +106,21 @@ public:
 		return *this;
 	}
-	matrix<T>& operator=(matrix<T> rhs){
-		init(rhs.R, rhs.C);
-		memcpy(M, rhs.M, R * C * sizeof(T));
+	matrix<T>& operator=(const matrix<T>& rhs){
+		if (this != &rhs) {											//if the matrix isn't self-assigned
+			T* new_matrix = new T[rhs.R * rhs.C];					//allocate new resources
+			memcpy(new_matrix, rhs.M, rhs.R * rhs.C * sizeof(T));	//copy the matrix
+
+			delete[] M;												//delete the previous array
+			M = new_matrix;
+			R = rhs.R;
+			C = rhs.C;
+		}
 		return *this;
 	}
 	//element-wise operations
-	matrix<T> operator+(T rhs) {
+	matrix<T> operator+(const T rhs) const {
 		matrix<T> result(R, C);					//create a result matrix
 		size_t N = R * C;
@@ -88,11 +130,39 @@ public:
 		return result;
 	}
-	matrix<T> operator-(T rhs) {
+	matrix<T> operator+(const matrix<T> rhs) const {
+		if (R != rhs.R || C != rhs.C) {
+			std::cout << "ERROR: addition is only defined for matrices that are the same size." << std::endl;
+			exit(1);
+		}
+		matrix<T> result(R, C);					//create a result matrix
+		size_t N = R * C;
+
+		for (int i = 0; i < N; i++)
+			result.M[i] = M[i] + rhs.M[i];			//calculate the operation and assign to result
+
+		return result;
+	}
+
+	matrix<T> operator-(const T rhs) const {
 		return operator+(-rhs);					//add the negative of rhs
 	}
-	matrix<T> operator*(T rhs) {
+	matrix<T> operator-(const matrix<T> rhs) const {
+		return operator+(-rhs);
+	}
+
+	matrix<T> operator-() const {
+		matrix<T> result(R, C);					//create a result matrix
+		size_t N = R * C;
+
+		for (int i = 0; i < N; i++)
+			result.M[i] = -M[i];			//calculate the operation and assign to result
+
+		return result;
+	}
+
+	matrix<T> operator*(const T rhs) const {
 		matrix<T> result(R, C);					//create a result matrix
 		size_t N = R * C;
@@ -102,7 +172,7 @@ public:
 		return result;
 	}
-	matrix<T> operator/(T rhs) {
+	matrix<T> operator/(const T rhs) const {
 		matrix<T> result(R, C);					//create a result matrix
 		size_t N = R * C;
@@ -113,7 +183,7 @@ public:
 	}
 	//matrix multiplication
-	matrix<T> operator*(matrix<T> rhs){
+	matrix<T> operator*(const matrix<T> rhs) const {
 		if(C != rhs.R){
 			std::cout<<"ERROR: matrix multiplication is undefined for matrices of size ";
 			std::cout<<"[ "<<R<<" x "<<C<<" ] and [ "<<rhs.R<<" x "<<rhs.C<<"]"<<std::endl;
@@ -127,7 +197,7 @@ public:
 			for(r = 0; r < R; r++){
 				inner = (T)0;
 				for(i = 0; i < C; i++){
-					inner += at(r, i) * rhs.at(i, c);
+					inner += get(r, i) * rhs.get(i, c);
 				}
 				result.M[c * R + r] = inner;
 			}
@@ -141,7 +211,7 @@ public:
 	}
 	//return a transposed matrix
-	matrix<T> transpose(){
+	matrix<T> transpose() const {
 		matrix<T> result(C, R);
 		size_t c, r;
 		for(c = 0; c < C; c++){
@@ -152,8 +222,31 @@ public:
 		return result;
 	}
+	///Calculate and return the determinant of the matrix
+	T det() const {
+		if (R != C) {
+			std::cout << "ERROR: a determinant can only be calculated for a square matrix." << std::endl;
+			exit(1);
+		}
+		if (R == 1) return M[0];			//if the matrix only contains one value, return it
+
+		int r, c, ri, cia, cib;
+		T a = 0;
+		T b = 0;
+		for (c = 0; c < C; c++) {
+			for (r = 0; r < R; r++) {
+				ri = r;
+				cia = (r + c) % C;
+				cib = (C - 1 - r) % C;
+				a += get(ri, cia);
+				b += get(ri, cib);
+			}
+		}
+		return a - b;
+	}
+
 	/// Sort rows of the matrix by the specified indices
-	matrix<T> sort_rows(size_t* idx) {
+	matrix<T> sort_rows(size_t* idx) const {
 		matrix<T> result(C, R);					//create the output matrix
 		size_t r, c;
 		for (c = 0; c < C; c++) {								//for each column
@@ -165,7 +258,7 @@ public:
 	}
 	/// Sort columns of the matrix by the specified indices
-	matrix<T> sort_cols(size_t* idx) {
+	matrix<T> sort_cols(size_t* idx) const {
 		matrix<T> result(C, R);
 		size_t c;
 		for (c = 0; c < C; c++) {											//for each column
@@ -174,7 +267,7 @@ public:
 		return result;
 	}
-	std::string toStr() {
+	std::string toStr() const {
 		std::stringstream ss;
 		for(int r = 0; r < R; r++) {
@@ -187,25 +280,35 @@ public:
 		return ss.str();
 	}
-	std::string csv() {
-		std::stringstream csvss;
+	void csv(std::ostream& out) const {
+		//std::stringstream csvss;
 		for (size_t i = 0; i < R; i++) {
-			csvss << M[i];
-			for (size_t j = 0; j < C; j++) csvss << ", " << M[j * R + i];
-			csvss << std::endl;
+			out << std::fixed << M[i];
+			for (size_t j = 1; j < C; j++)
+				out << ", " << std::fixed << M[j * R + i];
+			out << std::endl;
 		}
+		//return csvss.str();
+	}
+
+	std::string csv() const {
+		std::stringstream csvss;
+		int digits = std::numeric_limits<double>::max_digits10;
+		csvss.precision(digits);
+		csv(csvss);
 		return csvss.str();
 	}
+
+
 	//save the data as a CSV file
-	void csv(std::string filename) {
+	void csv(std::string filename) const {
 		ofstream basisfile(filename.c_str());
 		basisfile << csv();
 		basisfile.close();
 	}
 	static matrix<T> I(size_t N) {
-
 		matrix<T> result(N, N);							//create the identity matrix
 		memset(result.M, 0, N * N * sizeof(T));			//set the entire matrix to zero
 		for (size_t n = 0; n < N; n++) {
@@ -214,6 +317,19 @@ public:
 		return result;
 	}
+	//loads a matrix from a stream in CSV format
+	void csv(std::istream& in) {
+		size_t c, r;
+		T v;
+		for (r = 0; r < R; r++) {
+			for (c = 0; c < C; c++) {
+				in >> v;
+				if (in.peek() == ',') in.seekg(1, std::ios::cur);
+				at(r, c) = v;;
+			}
+		}
+	}
+
 };
 }	//end namespace rts