Verified Commit 9b53b2cd authored by Tim Gymnich's avatar Tim Gymnich
Browse files

added matrix.h

parent b00263fe
......@@ -117,6 +117,7 @@ set(PAIN_HEADERS
src/global.h
src/abstract_state.h
src/hash_utils.h
src/matrix.h
)
include_directories(${LLVM_INCLUDE_DIRS})
......@@ -157,6 +158,16 @@ target_link_libraries(normalized_conjunction_test
${LLVM_AVAILABLE_LIBS}
)
add_llvm_executable(matrix_test
test/matrix_test.cpp
${PAIN_HEADERS}
${PAIN_SOURCES}
)
target_link_libraries(matrix_test
${LLVM_AVAILABLE_LIBS}
)
enable_testing()
add_test(NAME intervalAnalysisTest
......@@ -171,6 +182,10 @@ add_test(NAME normalizedConjunctionTest
COMMAND normalized_conjunction_test
)
add_test(NAME matrixTest
COMMAND matrix_test
)
set(SAMPLES
add-1-float
add-1
......
#pragma once
#include "global.h"
#include <vector>
using namespace std;
namespace pcpo {
/// Matrix. Row and column are indexed beginning at 0
template<typename T> class Matrix {
public:
//MARK: - Initializers
/// Creates a matrix with dimensions height x width and initalizes its values to `value`
/// @param width Width of the matrix
/// @param height Height of the matrix
/// @param value Initial value for each element
Matrix<T>(int height, int width, T value) {
this->width = width;
this->height = height;
this->vectors.reserve(width);
for (int i = 0; i < height; i++) {
vector<T> vector(width,value);
vectors.push_back(vector);
}
};
/// Creates a matrix with dimensions height x width and initalizes its values to 0
/// @param width Width of the matrix
/// @param height Height of the matrix
Matrix<T>(int height, int width) : Matrix<T>(height, width, 0) {};
Matrix(Matrix const& matrix) = default;
/// Creates an identity matrix with dimension eye x eye
/// @param eye dimension of the matrix
Matrix(int eye){
this->width = eye;
this->height = eye;
this->vectors.reserve(width);
for (int i = 0; i < height; i++) {
vector<T> vector(width,0);
vector[i] = 1;
vectors[i] = vector;
}
};
/// Creates a matrix from a 2D vector
/// @param vectors 2D vector containing columns with rows
Matrix(vector<vector<T>> const &vectors) {
this->height = vectors.size();
this->width = vectors.front().size();
this->vectors = vectors;
};
/// Creates a vector from a std::vector
/// @param vector the vector
Matrix(vector<T> const &vector) {
std::vector<std::vector<T>> vectors = {vector};
this->vectors = vectors;
};
// MARK: - Properties
/// The height of the matrix (number of rows)
int getHeight() const { return height; };
/// The width of the matrix (number of columns)
int getWidth() const { return width; };
// MARK: - Matrix operations
/// Transposes the matrix
Matrix transpose() const {
Matrix result = Matrix(width, height);
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
result.vectors[j][i] = vectors[i][j];
}
}
return result;
};
/// Transforms the matrix to reduced row echelon form
Matrix echelon() const {
Matrix result = Matrix(*this);
int pivot = 0;
for (int row = 0; row < height; row++) {
if (pivot >= width) { return result; }
int i = row;
while (result.vectors[i][pivot] == 0) {
if (++i >= height) {
i = row;
if (++pivot >= width) { return result; }
}
}
result.swap_rows(i, row);
result.divide_row(row, result.vectors[row][pivot]);
for (i = 0; i < height; i++) {
if (i != row) {
result.add_multiple_row(i, row, -result.vectors[i][pivot]);
}
}
}
return result;
};
/// Linear span of the matrix ... fixme
Matrix span() const;
/// Returns a vector with the elements of the row at index i. The returned row cannot be modified.
/// @param i Index of the row to return.
vector<T>& row(int i) { return vectors[i]; };
/// Returns the column at index i. The returned column can be modified.
/// @param i Index of the column to return
vector<T> column(int i) const {
vector<T> row;
row.reserve(width);
for (vector<T> const& x : vectors) {
row.push_back(x[i]);
}
return row;
}
// MARK: - Operators
T& operator()(int i, int j) { return vectors[i][j]; };
Matrix<T>& operator *=(Matrix<T> const& rhs) {
assert(width == rhs.height);
Matrix result = Matrix(height,rhs.width);
for (int i = 0; i < height; i++) {
for (int k = 0; k < width; k++) {
for (int j = 0; j < rhs.width; j++) {
result.vectors[i][j] += vectors[i][k] * rhs.vectors[k][j];
}
}
}
this->vectors = result.vectors;
this->width = result.width;
this->height = result.height;
return *this;
};
Matrix<T>& operator *=(T scalar) {
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
vectors[i][j] *= scalar;
}
}
return *this;
};
Matrix<T>& operator +=(Matrix<T> const& rhs) {
assert(rhs.width == width && rhs.height == height);
for (int i=0;i<height;i++) {
for (int j = 0; j < width; j++) {
vectors[i][j] += rhs.vectors[i][j];
}
}
return *this;
};
Matrix<T>& operator +=(T scalar) {
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
vectors[i][j] += scalar;
}
}
return *this;
};
Matrix<T>& operator -=(Matrix<T> const& rhs) {
assert(rhs.width == width && rhs.height == height);
for (int i = 0; i < height; i++) {
for (int j= 0; j < width; j++) {
vectors[i][j] -= rhs.vectors[i][j];
}
}
return *this;
};
Matrix<T>& operator -=(int scalar) {
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
vectors[i][j] -= scalar;
}
}
return *this;
};
bool operator==(Matrix<T> const& rhs) const {
return rhs.vectors == rhs.vectors && width == rhs.width && height == rhs.height;
};
protected:
vector<vector<T>> vectors;
int width;
int height;
// MARK: - Echelon helpers
/// Swaps two rows
/// @param a index of the first row
/// @param b index of the second row
void swap_rows(int a, int b) { vectors[a].swap(vectors[b]); };
/// Divides a row by a constant
/// @param row index of the row to divide
/// @param quotient quotient to divide the row by
void divide_row(int row, int quotient) {
for (int column = 0; column < width; column++) {
vectors[row][column] /= quotient;
}
};
/// Adds a multiple of row b to a row a
/// @param a Row to add a multiple of b to
/// @param b Row to be added to row a
/// @param factor Factor to multiply row b with when adding it to row a
void add_multiple_row(int a, int b, int factor) {
for (int column = 0; column < width; column++) {
vectors[a][column] += vectors[b][column] * factor;
}
};
// MARK: - Utils
/// Greates common divisor.
/// @param lhs
/// @param rhs
static int ggT(int lhs, int rhs) {
int h;
if (lhs == 0) { return abs(rhs); }
if (rhs == 0) { return abs(lhs); }
do {
h = lhs % rhs;
lhs = rhs;
rhs = h;
} while (rhs != 0);
return abs(lhs);
};
/// Least common multiple.
/// @param lhs
/// @param rhs
static int kgV(int lhs, int rhs) { return (lhs * rhs) / ggT(lhs, rhs); };
};
template <typename T>
inline Matrix<T> operator*(Matrix<T> lhs, Matrix<T> const& rhs) { return lhs *= rhs; };
template <typename T>
inline Matrix<T> operator*(Matrix<T> lhs, T scalar) { return lhs *= scalar; };
template <typename T>
inline Matrix<T> operator+(Matrix<T> lhs, Matrix<T> const& rhs) { return lhs += rhs; };
template <typename T>
inline Matrix<T> operator+(Matrix<T> lhs, T scalar) { return lhs += scalar; };
template <typename T>
inline Matrix<T> operator-(Matrix<T> lhs, Matrix<T> const& rhs) { return lhs -= rhs; };
template <typename T>
inline Matrix<T> operator-(Matrix<T> lhs, T scalar) { return lhs -= scalar; };
}
#include <cstdio>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>
#include "../src/matrix.h"
using namespace std;
using namespace pcpo;
template <typename T>
class MatrixTest: Matrix<T> {
public:
static bool runTestMul1();
static bool runTestMul2();
static bool runTestTranspose1();
static bool runTestTranspose2();
static bool runTestEchelon1();
static bool runTestEchelon2();
static bool runTestEchelon3();
};
template <typename T>
bool MatrixTest<T>::runTestMul1() {
std::cout << "Testing multiplication 1: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,4,7},
{2,5,8},
{3,6,9}
};
std::vector<std::vector<T>> b = {
{4,29,0},
{-1,27,2},
{100,5,3}
};
std::vector<std::vector<T>> expected = {
{700,172,29},
{803,233,34},
{906,294,39}
};
auto actual = Matrix(a) * Matrix(b);
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
template <typename T>
bool MatrixTest<T>::runTestMul2() {
std::cout << "Testing multiplication 2: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,6,11},
{2,7,12},
{3,8,13},
{4,9,14},
{5,10,-9}
};
std::vector<std::vector<T>> b = {
{43,45,1,9},
{224,7,-2,24},
{12,1,13,-6}
};
std::vector<std::vector<T>> expected = {
{1519,87,132,87},
{1798,139,144,114},
{2077,191,156,141},
{2356,243,168,168},
{2347,295,-132,339}
};
auto actual = Matrix(a) * Matrix(b);
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
template<typename T>
bool MatrixTest<T>::runTestTranspose1() {
std::cout << "Testing transpose 1: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,2,3},
{4,5,6},
{7,8,9}
};
std::vector<std::vector<T>> expected = {
{1,4,7},
{2,5,8},
{3,6,9}
};
auto matrix = Matrix(a);
auto actual = matrix.transpose();
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
template<typename T>
bool MatrixTest<T>::runTestTranspose2() {
std::cout << "Testing transpose 2: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,3},
{2,4},
{3,5}
};
std::vector<std::vector<T>> expected = {
{1,2,3},
{4,5,6}
};
auto matrix = Matrix(a);
auto actual = matrix.transpose();
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
template<typename T>
bool MatrixTest<T>::runTestEchelon1() {
std::cout << "Testing echelon 1: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,4,7},
{2,5,8},
{3,6,9}
};
std::vector<std::vector<T>> expected = {
{1,0,-1},
{0,1,2},
{0,0,0}
};
auto matrix = Matrix(a);
auto actual = matrix.echelon();
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
template<typename T>
bool MatrixTest<T>::runTestEchelon2() {
std::cout << "Testing echelon 2: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,2,1},
{1,4,8},
{1,6,3}
};
std::vector<std::vector<T>> expected = {
{1,0,0},
{0,1,0},
{0,0,1}
};
auto matrix = Matrix(a);
auto actual = matrix.echelon();
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
template<typename T>
bool MatrixTest<T>::runTestEchelon3() {
std::cout << "Testing echelon 3: ";
bool result = false;
std::vector<std::vector<T>> a = {
{1,2,4},
{2,4,8},
{4,8,16}
};
std::vector<std::vector<T>> expected = {
{1,2,4},
{0,0,0},
{0,0,0}
};
auto matrix = Matrix(a);
auto actual = matrix.echelon();
auto x = Matrix(expected);
result = actual == x;
std::cout << (result? "success" : "failed") << "\n";
return result;
}
int main() {
return !(MatrixTest<int>::runTestMul1()
&& MatrixTest<int>::runTestMul2()
&& MatrixTest<int>::runTestTranspose1()
&& MatrixTest<int>::runTestTranspose2()
&& MatrixTest<int>::runTestEchelon1()
&& MatrixTest<int>::runTestEchelon2()
&& MatrixTest<int>::runTestEchelon3()
);
};
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment