Merge branch 'main' into first_half

make x column

.gitignore

@@ -0,0 +1,33 @@
+# Prerequisites
+*.d
+CMakeLists.txt
+.idea
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app

main.cpp

@@ -4,36 +4,16 @@
 #include "simplex.h"
 
 int main() {
-    // TODO: Initially should be positive
+    ColumnVector C(6);
-
+    Matrix A(6, 4);
-    Vector C = {-5, -4, 0, 0, 0, 0};
+    int Carr[C.getColumns()] = {-5, -4, 0, 0, 0, 0};
-
+    for (int i = 0; i < C.getRows(); ++i) {
-    Matrix A = {
+        C[i] = Carr[i];
-             {6, 4, 1, 0, 0, 0},
-             {1, 2, 0, 1, 0, 0},
-             {-1, 1, 0, 0, 1, 0},
-             {0, 1, 0, 0, 0, 1}
-                };
-
-    Vector b = {24, 6, 1, 2};
-
-    //Matrix test = {{1, -1, -2}, {1, 1, -2}, {1, -1, 2}};
-
-    //showMatrix(test);
-    
-
-    Result result = Simplex(C, A, b);
-    if(result.state == bounded) {
-
-        if(result.solution == nullptr) {
-            // TODO: error
-        }
-        std::cout << *result.solution << std::endl;
-        std::cout << result.objective_fucntion_value << std::endl;
-
-        delete result.solution;
     }
-    //std::cout << "asdasd" << std::endl;
+    A = {{6, 4, 1, 0, 0, 0}, {1, 2, 0, 1, 0, 0}, {-1, 1, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 1}};
-
+    ColumnVector b(A.getRows());
+    b = {24, 6, 1, 2};
+    std::cout << A << std::endl;
+    std::cout << C << std::endl;
     return 0;
 }

simplex.cpp

@@ -2,7 +2,6 @@
 #include <iostream>
 #include "tools/matrix.h"
 #include "tools/math.h"
-#include "tools/elimination.h"
 
 enum solver_state {
 	unbounded,
@@ -11,77 +10,55 @@ enum solver_state {
 
 struct Result {
 	solver_state state;
-	Vector *solution;
+	ColumnVector *solution;
 	double objective_function_value;
 };
 
-
+struct FracturedMatrix {
-
+	Matrix A;
-Result Simplex(Vector C, Matrix A, Vector b, double eps = 0.01, bool maximize=true) {
+	ColumnVector C;
-	if (maximize == false) {
+	ColumnVector b;
-		for (int i = 0; i < C.size(); i++) {
+	int pivot_column_index;
-		    C[i] = -C[i];
+	int pivot_row_index;
-		    }
+	FracturedMatrix();
+	FracturedMatrix(const FracturedMatrix& other) : A(other.A), C(other.C), b(other.b), pivot_column_index(other.pivot_column_index), pivot_row_index(other.pivot_row_index) {}
+	FracturedMatrix(Matrix A, ColumnVector C, ColumnVector b, int pivot_column_index, int pivot_row_index) : A(A), C(C), b(b), pivot_column_index(pivot_column_index), pivot_row_index(pivot_row_index) {}
+	FracturedMatrix& FracturedMatrix::operator=(const FracturedMatrix& other) {
+		A = other.A;
+        C = other.C;
+        b = other.b;
+        pivot_column_index = other.pivot_column_index;
+        pivot_row_index = other.pivot_row_index;
+		return *this;
 	}
- 	Result result{};
+};
-  	Matrix generalMatrix = createGeneralMatrix(A, C, b);
-	std::cout << "Before:" << std::endl;
-	showMatrix(generalMatrix);
-	std::vector<int> basicVars(generalMatrix.getRows());
-	basicVars[0] = -1;
 
-	for (size_t i = 1; i < basicVars.size(); i++) {
+Result Simplex(ColumnVector C, Matrix A, ColumnVector b, double eps = 0.01, bool maximize=true) {
-		basicVars[i] = static_cast<int>(basicVars.size()) + i;
-	}
 
-	while (true) {
+	std::cout << "C: " << C << std::endl;
+	std::cout << "A: " << A << std::endl;
+    std::cout << "b: " << b << std::endl;
+
+    int n = C.getRows();
+    int m = A.getColumns();
+
+	while (b[0] < eps) {
 		//3
 		int pivot_column_index = 0;
-		pivot_column_index = min_index(generalMatrix[0]);
+		pivot_column_index = Math::max_index(C);
-
-		if (generalMatrix[0][pivot_column_index] >= 0) {
-			DestroyMatrix destroyedGeneralMatrix = destroyGeneralMatrix(generalMatrix);
-			Matrix _A = destroyedGeneralMatrix.A;
-			Vector _C = destroyedGeneralMatrix.C;
-			Vector _b = destroyedGeneralMatrix.b;
-
-			result.state = bounded;
-			result.solution = new Vector(C.size());
-			for (int i = 0; i < C.size(); i++) {
-				result.solution->operator[](i) = 0;
-			}
-			for (size_t i = 1; i < basicVars.size(); i++) {
-				if (basicVars[i] <= C.size()) {
-					result.solution->operator[](basicVars[i]) = _b[i];
-				}
-			}
-			result.objective_function_value = b[0];
-			return result;
-		}
 
 		//4
-		Vector ratio_vector(generalMatrix.getRows());
+		ColumnVector ratio_vector(m);
-		for (int i = 1; i < generalMatrix.getRows(); i++) {
+		for (int i = 0; i < m; i++) {
-			if (generalMatrix[i][pivot_column_index] != 0) {
+			ratio_vector[i] = b[i] / A[i][pivot_column_index];
-            	ratio_vector[i] = generalMatrix[i][generalMatrix.getColumns() - 1] / generalMatrix[i][pivot_column_index];
-        	} else {
-            	ratio_vector[i] = 0;
-        	}
 		}
-		ratio_vector[0] = 0;
+		int pivot_row_index = Math::min_index(ratio_vector);
-
-		int pivot_row_index = min_index_positive(ratio_vector);
-
-		basicVars[pivot_row_index] = pivot_column_index;
 
 		//5
-		elimination(generalMatrix, pivot_row_index, pivot_column_index);
+		struct FracturedMatrix fractured_matrix();
-		std::cout << "After:" << std::endl;
+		fractured_matrix = eleminate(A, C, b, pivot_column_index, pivot_row_index);
-		showMatrix(generalMatrix);
 	}
 
-	return result;
-
 	/*
 	Result result;
 	std::vector<int> basicVars(A.getColumns() - A.getRows());
@@ -102,7 +79,7 @@ Result Simplex(Vector C, Matrix A, Vector b, double eps = 0.01, bool maximize=tr
 
 	if (A[0][kc] >= 0) {
 		result.state = unbounded;
-		result.solution = new Vector(C.getRows());
+		result.solution = new ColumnVector(C.getRows());
 		for (int i = 0; i < C.getRows(); i++) {
 			result.solution->operator[](i) = 0;
 		}

simplex.h

@@ -11,10 +11,10 @@ enum solver_state {
 
 struct Result {
   solver_state state;
-  Vector *solution;
+  ColumnVector *solution;
   double objective_fucntion_value;
 };
 
-Result Simplex(Vector C, Matrix A, Vector b, double eps = 0.01, bool maximize = true);
+Result Simplex(ColumnVector C, Matrix A, ColumnVector b, double eps = 0.01, bool maximize = true);
 
 #endif // SIMPLEX_H

tools/elimination.cpp

@@ -1,24 +1,28 @@
 #include "elimination.h"
+#include "matrix.h"
 #include "math.h"
 
-FracturedMatrix::FracturedMatrix(const FracturedMatrix& other) : A(other.A), C(other.C), b(other.b), pivot_column_index(other.pivot_column_index), pivot_row_index(other.pivot_row_index) {}
+struct FracturedMatrix {
-FracturedMatrix::FracturedMatrix(Matrix A, Vector C, Vector b, int pivot_column_index, int pivot_row_index) : A(A), C(C), b(b), pivot_column_index(pivot_column_index), pivot_row_index(pivot_row_index) {}
+    Matrix A;
-FracturedMatrix& FracturedMatrix::operator=(const FracturedMatrix& other) {
+    ColumnVector C;
-    A = other.A;
+    ColumnVector b;
-    C = other.C;
+    int pivot_column_index;
-    b = other.b;
+    int pivot_row_index;
-    pivot_column_index = other.pivot_column_index;
+};
-    pivot_row_index = other.pivot_row_index;
+
-    return *this;
+struct DestroyMatrix {
-}
+    Matrix A;
+    ColumnVector C;
+    ColumnVector b;
+};
 
 DestroyMatrix destroyGeneralMatrix(Matrix& generalMatrix) {
     int rows = generalMatrix.getRows();
     int cols = generalMatrix.getColumns();
 
     Matrix A(rows - 1, cols - 1);
-    Vector C(cols - 1);
+    ColumnVector C(cols - 1);
-    Vector b(rows - 1);
+    ColumnVector b(rows - 1);
 
     for (int j = 0; j < cols - 1; ++j) {
         C[j] = generalMatrix[0][j];
@@ -37,15 +41,14 @@ DestroyMatrix destroyGeneralMatrix(Matrix& generalMatrix) {
     return {A, C, b};
 }
 
-Matrix createGeneralMatrix(Matrix& A, Vector& C, Vector& b) {
+Matrix createGeneralMatrix(Matrix& A, ColumnVector& C, ColumnVector& b) {
     Matrix generalMatrix(A.getRows() + 1, A.getColumns() + 1);
     for (int i = 0; i < A.getColumns(); i++) {
         generalMatrix[0][i] = C[i];
     }
 
-    // For objective function (j=0) the value is set to zero automatically
+    for (int j = 0; j < A.getRows(); j++) {
-    for (int j = 1; j < A.getRows() + 1; j++) {
+        generalMatrix[j][A.getColumns()] = b[j];
-        generalMatrix[j][A.getColumns()] = b[j-1];
     }
 
     for (int i = 0; i < A.getRows(); i++) {
@@ -57,7 +60,8 @@ Matrix createGeneralMatrix(Matrix& A, Vector& C, Vector& b) {
     return generalMatrix;
 }
 
-void elimination(Matrix& generalMatrix, int pivot_row_index, int pivot_column_index) {
+FracturedMatrix elimination(Matrix A, ColumnVector C, ColumnVector b, int pivot_column_index, int pivot_row_index) {
+    Matrix generalMatrix = createGeneralMatrix(A, C, b);
 
     int rows = generalMatrix.getRows();
     int cols = generalMatrix.getColumns();
@@ -78,5 +82,20 @@ void elimination(Matrix& generalMatrix, int pivot_row_index, int pivot_column_in
             generalMatrix[i][j] -= pivotColumnCoefficient * generalMatrix[pivot_row_index][j];
         }
     }
+
+    DestroyMatrix destroyedMatrix = destroyGeneralMatrix(generalMatrix);
+
+    pivot_column_index = Math::max_index(destroyedMatrix.C);
+    
+    ColumnVector ratio_vector(A.getRows());
+
+    for (int i = 0; i < A.getRows(); i++) {
+      ratio_vector[i] = b[i] / A[i][pivot_column_index];
+    }
+
+    int pivot_row_index = Math::min_index(ratio_vector);
+
+    return {A, C, b, pivot_column_index, pivot_row_index};
 }
 
+

tools/elimination.h

@@ -1,30 +1,13 @@
 #ifndef ELIMINATION_H
 #define ELIMINATION_H
 
-#include "matrix.h"
 
 
-struct FracturedMatrix {
+class Elimination {
-	Matrix A;
+    public:
-	Vector C;
+    Matrix elimination(Matrix A, ColumnVector C, ColumnVector b, double eps);
-	Vector b;
-	int pivot_column_index;
-	int pivot_row_index;
-	FracturedMatrix() = default;
-	FracturedMatrix(const FracturedMatrix& other);
-	FracturedMatrix(Matrix A, Vector C, Vector b, int pivot_column_index, int pivot_row_index);
-	FracturedMatrix& operator=(const FracturedMatrix& other);
 };
 
-struct DestroyMatrix {
-    Matrix A;
-    Vector C;
-    Vector b;
-};
-
-DestroyMatrix destroyGeneralMatrix(Matrix& generalMatrix);
-void elimination(Matrix&, int pivot_row_index, int pivot_column_index);
-Matrix createGeneralMatrix(Matrix& A, Vector& C, Vector& b);
 
 
 #endif //ELIMINATION_H

tools/math.cpp

@@ -1,64 +1,45 @@
 #include "math.h"
+#include "matrix.h"
 
-double min(Vector vector) {
+double min(ColumnVector columnVector) {
-    double temp = vector[0];
+    double temp = columnVector[0];
-    for (int j = 0; j < vector.size(); j++) {
+    for (int j = 0; j < columnVector.getColumns(); j++) {
-        if (vector[j] < temp) {
+        if (columnVector[j] < temp) {
-            temp = vector[j];
+            temp = columnVector[j];
         }
     }
     return temp;
 }
 
-double max(Vector vector) {
+double max(ColumnVector columnVector) {
-    double temp = vector[0];
+    double temp = columnVector[0];
-    for (int j = 0; j < vector.size(); j++) {
+    for (int j = 0; j < columnVector.getColumns(); j++) {
-        if (vector[j] > temp) {
+        if (columnVector[j] > temp) {
-            temp = vector[j];
+            temp = columnVector[j];
         }
     }
     return temp;
 }
 
-int min_index_positive(Vector vector) {
+int min_index(ColumnVector columnVector) {
-    int i = 0;
+    int temp_index = 0;
-    while (i < vector.size() && vector[i] <= 0) {
+    double temp = columnVector[0];
-        ++i;
+    for (int j = 0; j < columnVector.getColumns(); j++) {
-    }
+        if (columnVector[j] < temp) {
-    if (i >= vector.size()) {
-        throw std::runtime_error("No positive min found");
-    }
-    int temp_index = i;
-    double temp = vector[i];
-
-    for (int j = i + 1; j < vector.size(); j++) {
-        if (vector[j] < temp && vector[j] > 0) {
             temp_index = j;
-            temp = vector[j];
+            temp = columnVector[j];
         }
     }
     return temp_index;
 }
 
-int min_index(Vector vector) {
+int max_index(ColumnVector columnVector) {
     int temp_index = 0;
-    double temp = vector[0];
+    double temp = columnVector[0];
-    for (int j = 0; j < vector.size(); j++) {
+    for (int j = 0; j < columnVector.getColumns(); j++) {
-        if (vector[j] < temp) {
+        if (columnVector[j] > temp) {
             temp_index = j;
-            temp = vector[j];
+            temp = columnVector[j];
-        }
-    }
-    return temp_index;
-}
-
-int max_index(Vector vector) {
-    int temp_index = 0;
-    double temp = vector[0];
-    for (int j = 0; j < vector.size(); j++) {
-        if (vector[j] > temp) {
-            temp_index = j;
-            temp = vector[j];
         }
     }
     return temp_index;

tools/math.h

@@ -2,13 +2,11 @@
 #define MATH_H
 
 #include <vector>
-#include "matrix.h"
 
-double min(Vector vector);
+double min(ColumnVector columnVector);
-double max(Vector vector);
+double max(ColumnVector columnVector);
-int min_index(Vector vector);
+int min_index(ColumnVector columnVector);
-int min_index_positive(Vector vector);
+int max_index(ColumnVector columnVector);
-int max_index(Vector vector);
 double min(std::vector<double> array);
 double max(std::vector<double> array);
 int min_index(std::vector<double> array);

tools/matrix.cpp

@@ -1,116 +1,68 @@
 #include "matrix.h"
 
-void showMatrix(Matrix matrix) {
+ColumnVector::ColumnVector(int n) {
-  int maxNumberLength = 1;
-  for (size_t y = 0; y < matrix.getRows(); y++) {
-    for (size_t x = 0; x < matrix.getColumns(); x++) {
-      if (std::to_string(matrix[y][x]).length() > maxNumberLength) {
-        maxNumberLength = std::to_string(matrix[y][x]).length();
-      }
-    }
-  }
-  std::cout << maxNumberLength << std::endl;
-  for (size_t y = 0; y < matrix.getRows(); y++) {
-    std::string row = "";
-    for (size_t x = 0; x < matrix.getColumns(); x++) {
-      std::string strNumber = std::to_string(matrix[y][x]);
-      bool spaceRight = true;
-      while (strNumber.length() < maxNumberLength) {
-        if (spaceRight) {
-          strNumber += " ";
-        }
-        else {
-          strNumber = " " + strNumber;
-        }
-        spaceRight = !spaceRight;
-        strNumber = "" + strNumber;
-
-      }
-      if (matrix[y][x] == 0) {
-        row += "[\033[0m" + strNumber + "\033[0m] ";
-      }
-      else if(matrix[y][x] > 0) {
-        row += "[\033[32m" + strNumber + "\033[0m] ";
-      }else {
-        row += "[\033[31m" + strNumber + "\033[0m] ";
-      }
-
-    }
-    std::cout << row << std::endl;
-  }
-}
-
-Vector::Vector(int n) {
   rows = n;
   columns = 1;
-  vector.resize(n);
+  columnVector.resize(n);
 }
 
-Vector::Vector(const Vector& other) {
+ColumnVector::ColumnVector(const ColumnVector& other) {
   rows = other.rows;
   columns = other.columns;
-  vector = other.vector;
+  columnVector = other.columnVector;
 }
 
-Vector::Vector(std::initializer_list<double> init) {
+int ColumnVector::getRows() const {
-  rows = init.size();
-  columns = 1;
-  vector = std::vector<double>(rows);
-
-  auto it = init.begin();
-  
-  for (size_t i = 0; i < init.size(); ++i) {
-    vector[i] = *it++;
-  }
-}
-
-int Vector::size() const {
   return rows;
 }
 
-double& Vector::operator[](int row) {
+int ColumnVector::getColumns() const {
-  return vector[row];
+  return columns;
 }
 
-Vector& Vector::operator=(const Vector& other) {
+double& ColumnVector::operator[](int row) {
+  return columnVector[row];
+}
+
+ColumnVector& ColumnVector::operator=(const ColumnVector& other) {
   rows = other.rows;
   columns = other.columns;
-  vector = other.vector;
+  columnVector = other.columnVector;
   return *this;
 }
 
-Vector Vector::operator+(Vector& other) {
+ColumnVector ColumnVector::operator+(ColumnVector& other) {
   if (rows != other.rows) {
     throw std::runtime_error("Error: the dimensional problem occurred");
   }
-  Vector result(rows);
+  ColumnVector result(rows);
   for (int i = 0; i < rows; ++i) {
-    result[i] = vector[i] + other[i];
+    result[i] = columnVector[i] + other[i];
   }
 
   return result;
 }
 
-Vector Vector::operator-(Vector& other) {
+ColumnVector ColumnVector::operator-(ColumnVector& other) {
   if (rows != other.rows) {
     throw std::runtime_error("Error: the dimensional problem occurred");
   }
-  Vector result(rows);
+  ColumnVector result(rows);
   for (int i = 0; i < rows; ++i) {
-    result[i] = vector[i] - other[i];
+    result[i] = columnVector[i] - other[i];
   }
 
   return result;
 }
 
-std::istream& operator>>(std::istream& cin, Vector& vectorObj) {
+std::istream& operator>>(std::istream& cin, ColumnVector& vectorObj) {
   for (int i = 0; i < vectorObj.rows; ++i) {
     cin >> vectorObj[i];
   }
   return cin;
 }
 
-std::ostream& operator<<(std::ostream& cout, Vector& vectorObj) {
+std::ostream& operator<<(std::ostream& cout, ColumnVector& vectorObj) {
   for (int i = 0; i < vectorObj.rows; ++i) {
     if (i == vectorObj.rows - 1) {
       cout << vectorObj[i] << std::endl;
@@ -125,9 +77,9 @@ std::ostream& operator<<(std::ostream& cout, Vector& vectorObj) {
 Matrix::Matrix(int n, int m) {
   rows = n;
   columns = m;
-  matrix.resize(n, Vector(m));
+  matrix.resize(n, ColumnVector(m));
   for (auto& row : matrix) {
-    row = Vector(m);
+    row = ColumnVector(m);
   }
 }
 
@@ -137,20 +89,6 @@ Matrix::Matrix(const Matrix& other) {
   matrix = other.matrix;
 }
 
-Matrix::Matrix(std::initializer_list<std::vector<double>> init) {
-  rows = init.size();
-  auto it = init.begin();
-  columns = it->size();
-  matrix = std::vector<Vector>(rows, Vector(columns));
-
-  for (int i = 0; i < rows; ++i) {
-    for (int j = 0; j < columns; ++j) {
-      matrix[i][j] = it->operator[](j);
-    }
-    ++it;
-  }
-}
-
 int Matrix::getRows() const {
   return rows;
 }
@@ -159,7 +97,7 @@ int Matrix::getColumns() const {
   return columns;
 }
 
-Vector& Matrix::operator[](int row) {
+ColumnVector& Matrix::operator[](int row) {
   return matrix[row];
 }
 
@@ -222,12 +160,12 @@ Matrix Matrix::operator*(Matrix& other) const {
   return result;
 }
 
-Vector Matrix::operator*(Vector other) const {
+ColumnVector Matrix::operator*(ColumnVector other) const {
-  if (columns != other.size()) {
+  if (columns != other.getRows()) {
     throw std::runtime_error("Error: the dimensional problem occurred");
   }
 
-  Vector result(rows);
+  ColumnVector result(rows);
   for (int i = 0; i < rows; ++i) {
     result[i] = 0;
     for (int k = 0; k < columns; ++k) {

tools/matrix.h

@@ -6,44 +6,43 @@
 #include <vector>
 
 /**
-* Vector is a class to represent
+* ColumnVector is a class to represent
 * a column vector with n rows.
 */
-
+class ColumnVector {
-
-class Vector {
 protected:
 	// Number of rows in vector
 	int rows;
 	// Number of columns in vector
 	int columns;
 	// Matrix representation as vector of vectors of integers
-	std::vector<double> vector;
+	std::vector<double> columnVector;
 public:
-	Vector(int n);
+	ColumnVector(int n);
 
-	Vector(const Vector& other);
+	ColumnVector(const ColumnVector& other);
-
-	Vector(std::initializer_list<double>);
 
 	/* Getter for the number of rows */
-	int size() const;
+	int getRows() const;
+
+	/* Getter for the number of columns */
+	int getColumns() const;
 
 	double& operator[](int row);
 
-	Vector& operator=(const Vector& other);
+	ColumnVector& operator=(const ColumnVector& other);
 
-	Vector operator+(Vector& other);
+	ColumnVector operator+(ColumnVector& other);
 
-	Vector operator-(Vector& other);
+	ColumnVector operator-(ColumnVector& other);
 
 	/* Input operator reads element of vector */
-	friend std::istream& operator>>(std::istream& cin, Vector& vectorObj);
+	friend std::istream& operator>>(std::istream& cin, ColumnVector& vectorObj);
 
 	/* Output operator prints elements of the vector
 	* in a row separated with a space (no space at the end of the line)
 	*/
-	friend std::ostream& operator<<(std::ostream& cout, Vector& vectorObj);
+	friend std::ostream& operator<<(std::ostream& cout, ColumnVector& vectorObj);
 };
 
 /**
@@ -58,21 +57,19 @@ protected:
 	// Number of columns in matrix
 	int columns;
 	// Matrix representation as vector of vectors of integers
-	std::vector<Vector> matrix;
+	std::vector<ColumnVector> matrix;
 public:
 	Matrix(int n, int m);
 
 	Matrix(const Matrix& other);
 
-	Matrix(std::initializer_list<std::vector<double>>);
-
 	/* Getter for the number of rows */
 	int getRows() const;
 
 	/* Getter for the number of columns */
 	int getColumns() const;
 
-	Vector& operator[](int row);
+	ColumnVector& operator[](int row);
 
 	Matrix& operator=(const Matrix& other);
 
@@ -83,7 +80,7 @@ public:
 	Matrix operator*(Matrix& other) const;
 
 	/* Matrix-Vector multiplication */
-	Vector operator*(Vector other) const;
+	ColumnVector operator*(ColumnVector other) const;
 
 	/* Produces transposed version of the matrix */
 	Matrix transpose() const;
@@ -97,6 +94,4 @@ public:
 	friend std::ostream& operator<<(std::ostream& cout, Matrix& matrixObj);
 };
 
-void showMatrix(Matrix);
-
 #endif  // TOOLS_MATRIX_H