diff --git a/linalg/Makefile b/linalg/Makefile new file mode 100644 index 0000000..1b8aaef --- /dev/null +++ b/linalg/Makefile @@ -0,0 +1,20 @@ +CXXFLAGS = -Wall -g -std=c++1z -fpic +LDLIBS = +OBJ = vector.o vector_io.o matrix.o matrix_io.o block.o + +all: linalg_test + +linalg_test: linalg_test.cc liblinalg.so + g++ -Wall -llinalg -lboost_unit_test_framework -o linalg_test linalg_test.cc + +liblinalg.so: ${OBJ} linalg.h + g++ -shared -o liblinalg.so ${OBJ} + +install: liblinalg.so linalg.h + cp liblinalg.so /usr/local/lib + cp linalg.h /usr/local/include + +clean: + rm -f liblinalg.so + rm -f linalg_test + rm -f *.o diff --git a/linalg/README b/linalg/README new file mode 100644 index 0000000..2b32c6e --- /dev/null +++ b/linalg/README @@ -0,0 +1,3 @@ +# Vectors & Matrices + +A toy implementation of vector and matrix classes for linear algebra operations. diff --git a/linalg/block.cc b/linalg/block.cc new file mode 100644 index 0000000..2e02b72 --- /dev/null +++ b/linalg/block.cc @@ -0,0 +1,78 @@ +#include +#include +#include "linalg.h" + +Block::Block() : size_(0), data_(nullptr) {} + +Block::Block(int size) : Block() { + size_ = size; + data_ = new double[size_]; +} + +Block::Block(const Block& b) : Block(b.size_) { + memcpy(data_, b.data_, sizeof(double) * size_); +} + +Block::~Block() { + if (data_) + delete[] data_; +} + +void Block::setAll(double value) { + for (int i = 0; i < size_; i++) + data_[i] = value; +} + +bool Block::isEqual(Block& m) { + if (size_ != m.size_) + return false; + else if (memcmp(data_, m.data_, size_ * sizeof(double)) != 0) + return false; + return true; +} + +void Block::assign(Block& b) { + if (size_ != b.size_) { + if (data_) + delete[] data_; + size_ = b.size_; + data_ = new double[size_]; + } + memcpy(data_, b.data_, size_ * TYPE_SIZE); +} + +void Block::operateScalar(double value, int op) { + for (int i = 0; i < size_; i++) + switch (op) { + case SUM: + data_[i] += value; + break; + case DIFF: + data_[i] -= value; + break; + case PROD: + data_[i] *= value; + break; + case DIV: + data_[i] /= value; + break; + } +} + +void Block::operateBlock(Block& b, int op) { + for (int i = 0; i < size_; i++) + switch (op) { + case SUM: + data_[i] += b.data_[i]; + break; + case DIFF: + data_[i] -= b.data_[i]; + break; + case PROD: + data_[i] *= b.data_[i]; + break; + case DIV: + data_[i] /= b.data_[i]; + break; + } +} diff --git a/linalg/linalg.h b/linalg/linalg.h new file mode 100644 index 0000000..5e54c74 --- /dev/null +++ b/linalg/linalg.h @@ -0,0 +1,115 @@ +#ifndef _LINALG_H +#define _LINALG_H 1 + +#include + +#define GET(a, i, j) ((a).data_[(i) * ((a).cols_) + (j)]) + +const int TYPE_SIZE = sizeof(double); + +enum { + SUM = 1, + DIFF, + PROD, + DIV +}; + +class Block { + public: + Block(); + Block(int size); + Block(const Block& b); + ~Block(); + int size() { return size_; } + double* data() { return data_; } + void operateScalar(double value, int op); + void operateBlock(Block& b, int op); + void assign(Block& b); + void setAll(double value); + bool isEqual(Block& m); + protected: + int size_; + double* data_; +}; + +class Vector; + +class Matrix : public Block { + public: + Matrix(); + Matrix(const Matrix& m); + Matrix(int rows, int cols); + Matrix(std::initializer_list> list); + Matrix(const char *filename); + Matrix& operator=(Matrix m); + Matrix& operator=(std::initializer_list> list); + Matrix operator+(double value); + Matrix operator-(double value); + Matrix operator*(double value); + Matrix operator/(double value); + Matrix& operator+=(double value); + Matrix& operator-=(double value); + Matrix& operator*=(double value); + Matrix& operator/=(double value); + bool operator==(Matrix m); + bool operator!=(Matrix m); + Matrix operator+(Matrix &m); + Matrix& operator+=(Matrix &m); + Matrix operator-(Matrix &m); + Matrix& operator-=(Matrix &m); + Matrix operator*(Matrix &m); + Matrix t(); + const double& operator()(int row, int col) const; + double& operator()(int row, int col); + void read(std::ifstream& ifs); + void write(std::ofstream& ofs); + void writetxt(std::ofstream& ofs); + void fread(const char *filename); + void fwrite(const char *filename); + int rows() { return rows_; } + int cols() { return cols_; } + friend std::ostream& operator<<(std::ostream& os, Matrix& m); + private: + int rows_; + int cols_; + friend class Vector; +}; + +class Vector : public Block { + public: + Vector(); + Vector(const Vector& m); + Vector(int size_); + Vector(std::initializer_list list); + Vector(const char *filename); + Vector& operator=(Vector m); + Vector& operator=(std::initializer_list list); + Vector& operator=(double value); + Vector operator+(double value); + Vector operator-(double value); + Vector operator*(double value); + Vector operator/(double value); + Vector& operator+=(double value); + Vector& operator-=(double value); + Vector& operator*=(double value); + Vector& operator/=(double value); + bool operator==(Vector m); + bool operator!=(Vector m); + Vector operator+(Vector &m); + Vector& operator+=(Vector &m); + Vector operator-(Vector &m); + Vector& operator-=(Vector &m); + double operator*(Vector &m); + Vector t(); + const double& operator()(int index) const; + double& operator()(int index); + void read(std::ifstream& ifs); + void write(std::ofstream& ofs); + void fread(const char *filename); + void fwrite(const char *filename); + friend std::ostream& operator<<(std::ostream& os, Vector& m); + private: + bool coltype_; +}; + +#endif // _LINALG_H diff --git a/linalg/linalg_test.cc b/linalg/linalg_test.cc new file mode 100644 index 0000000..8fe00ff --- /dev/null +++ b/linalg/linalg_test.cc @@ -0,0 +1,210 @@ +#include "linalg.h" +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE matrix +#include + +BOOST_AUTO_TEST_CASE(MatrixConstructor) +{ + Matrix a; + BOOST_CHECK_EQUAL(a.rows(), 0); + BOOST_CHECK_EQUAL(a.cols(), 0); + BOOST_CHECK(a.data() == nullptr); + + Matrix b(3, 4); + BOOST_CHECK_EQUAL(b.rows(), 3); + BOOST_CHECK_EQUAL(b.cols(), 4); + BOOST_CHECK(b.data() != nullptr); +} + +BOOST_AUTO_TEST_CASE(ValueAssignment) +{ + Matrix a(3, 3); + a(0, 0) = 1; + BOOST_CHECK_EQUAL(*a.data(), 1); +} + +BOOST_AUTO_TEST_CASE(MatrixComparison) +{ + Matrix a({{1, 2}, {3, 4}}); + Matrix b = {{1, 2}, {3, 4}}; + BOOST_CHECK(a == b); + b(0, 0) = 0; + BOOST_CHECK(a != b); +} + +BOOST_AUTO_TEST_CASE(MatrixInitialization) +{ + Matrix a({{1, 2}, {3, 4}}); + BOOST_CHECK_EQUAL(a(1, 1), 4); + Matrix b = {{1, 2}, {3, 4}}; + BOOST_CHECK_EQUAL(b(1, 0), 3); + Matrix c; + c = {{1, 2}, {3, 4}}; + BOOST_CHECK_EQUAL(c(0, 1), 2); + Matrix d; + d = c; + BOOST_CHECK_EQUAL(d(0, 0), 1); + Matrix e = d; + BOOST_CHECK_EQUAL(e(1, 1), 4); + e(0, 0) = 10; + BOOST_CHECK_EQUAL(e(0, 0), 10); +} + +BOOST_AUTO_TEST_CASE(MatrixFilesOps) +{ + Matrix a({{1, 2}, {3, 4}}); + a.fwrite("test.dat"); + Matrix b; + b.fread("test.dat"); + // BOOST_CHECK(a == b); + Matrix c("test.dat"); + // BOOST_CHECK(a == c); +} + +BOOST_AUTO_TEST_CASE(MatrixTranspose) +{ + Matrix a({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}); + Matrix b = a.t(); + Matrix c = b.t(); + BOOST_CHECK(a == c); +} + +BOOST_AUTO_TEST_CASE(MatrixScalarArithmetics) +{ + Matrix a({{1, 2}, {3, 4}}); + Matrix b = a + 1.0; + BOOST_CHECK_EQUAL(b(0, 0), 2); + b = a - 1.0; + BOOST_CHECK_EQUAL(b(0, 1), 1); + b = a * 2; + BOOST_CHECK_EQUAL(b(1, 0), 6); + b = a / 2; + BOOST_CHECK_EQUAL(b(1, 1), 2); + a += 2; + BOOST_CHECK_EQUAL(a(0, 0), 3); + a -= 3; + BOOST_CHECK_EQUAL(a(0, 1), 1); + a *= 2; + BOOST_CHECK_EQUAL(a(1, 0), 4); + a /= 2; + BOOST_CHECK_EQUAL(a(1, 1), 3); +} + +BOOST_AUTO_TEST_CASE(MatrixMatrixArithmetics) +{ + Matrix a{{1, 2}, {-1, 4}}; + Matrix b{{2, 2}, {1, -1}}; + Matrix c = a + b; + BOOST_CHECK(c == Matrix({{3, 4}, {0, 3}})); + c = a - b; + BOOST_CHECK(c == Matrix({{-1, 0}, {-2, 5}})); + a += b; + BOOST_CHECK(a == Matrix({{3, 4}, {0, 3}})); + a -= b; + BOOST_CHECK(a == Matrix({{1, 2}, {-1, 4}})); +} + +BOOST_AUTO_TEST_CASE(MatrixProduct) +{ + Matrix a{{1, 2}, {-1, 4}}; + Matrix b{{2, 2, -1}, {1, -1, 2}}; + Matrix c = a * b; + BOOST_CHECK(c == Matrix({{4, 0, 3}, {2, -6, 9}})); +} + +BOOST_AUTO_TEST_CASE(VectorConstructor) +{ + Vector a; + BOOST_CHECK_EQUAL(a.size(), 0); + BOOST_CHECK(a.data() == nullptr); + + Vector b(3); + BOOST_CHECK_EQUAL(b.size(), 3); + BOOST_CHECK(b.data() != nullptr); +} + +BOOST_AUTO_TEST_CASE(VectorValueAssignment) +{ + Vector a(3); + a(0) = 1; + BOOST_CHECK_EQUAL(*a.data(), 1); +} + +BOOST_AUTO_TEST_CASE(VectorComparison) +{ + Vector a({1, 2, 3, 4}); + Vector b = {1, 2, 3, 4}; + BOOST_CHECK(a == b); + b(0) = 0; + BOOST_CHECK(a != b); +} + +BOOST_AUTO_TEST_CASE(VectorInitialization) +{ + Vector a({1, 2, 3, 4}); + BOOST_CHECK_EQUAL(a(3), 4); + Vector b = {1, 2, 3, 4}; + BOOST_CHECK_EQUAL(b(2), 3); + Vector c; + c = {1, 2, 3, 4}; + BOOST_CHECK_EQUAL(c(1), 2); + Vector d; + d = b; + BOOST_CHECK_EQUAL(d(0), 1); + Vector e = d; + BOOST_CHECK_EQUAL(e(3), 4); + e(0) = 10; + BOOST_CHECK_EQUAL(e(0), 10); +} + +BOOST_AUTO_TEST_CASE(VectorTranspose) +{ + Vector a({1, 2, 3}); + Vector b = a.t(); + Vector c = b.t(); + BOOST_CHECK(a == c); +} + +BOOST_AUTO_TEST_CASE(VectorScalarArithmetics) +{ + Vector a({1, 2, 3, 4}); + Vector b = a + 1.0; + BOOST_CHECK_EQUAL(b(0), 2); + b = a - 1; + BOOST_CHECK_EQUAL(b(1), 1); + b = a * 2; + BOOST_CHECK_EQUAL(b(2), 6); + b = a / 2; + BOOST_CHECK_EQUAL(b(3), 2); + a += 2; + BOOST_CHECK_EQUAL(a(0), 3); + a -= 3; + BOOST_CHECK_EQUAL(a(1), 1); + a *= 2; + BOOST_CHECK_EQUAL(a(2), 4); + a /= 2; + BOOST_CHECK_EQUAL(a(3), 3); +} + +BOOST_AUTO_TEST_CASE(VectorVectorArithmetics) +{ + Vector a{1, 2, -1, 4}; + Vector b{2, 2, 1, -1}; + Vector c = a + b; + BOOST_CHECK(c == Vector({3, 4, 0, 3})); + c = a - b; + BOOST_CHECK(c == Vector({-1, 0, -2, 5})); + a += b; + BOOST_CHECK(a == Vector({3, 4, 0, 3})); + a -= b; + BOOST_CHECK(a == Vector({1, 2, -1, 4})); +} + +BOOST_AUTO_TEST_CASE(VectorProduct) +{ + Vector a{1, 2, -1, 4}; + Vector b{2, 2, -1, 2}; + double c = a * b; + BOOST_CHECK(c == 15); +} + diff --git a/linalg/matrix.cc b/linalg/matrix.cc new file mode 100644 index 0000000..213bd75 --- /dev/null +++ b/linalg/matrix.cc @@ -0,0 +1,157 @@ +#include "linalg.h" + +Matrix::Matrix() : Block(), rows_(0), cols_(0) { +} + +Matrix::Matrix(int rows, int cols) + : Block(rows * cols), rows_(rows), cols_(cols) { +} + +Matrix::Matrix(const Matrix& m) + : Block(m), rows_(m.rows_), cols_(m.cols_) { +} + +Matrix::Matrix(std::initializer_list> list) + : Matrix() { + *this = list; +} + +Matrix::Matrix(const char *filename) : Block() { + fread(filename); +} + +Matrix& Matrix::operator=(Matrix m) { + if (this == &m) + return *this; + rows_ = m.rows_; + cols_ = m.cols_; + assign(m); + return *this; +} + + +Matrix& Matrix::operator=(std::initializer_list> list) { + rows_ = list.size(); + cols_ = list.begin()->size(); + int size_new = rows_ * cols_; + if (size_ != size_new) { + if (data_ != nullptr) + delete[] data_; + data_ = new double[size_new]; + } + size_ = size_new; + int i = 0; + for (std::initializer_list row : list) + for (double value : row) + data_[i++] = value; + return *this; +} + +bool Matrix::operator==(Matrix m) { + if (rows_ != m.rows_) + return false; + if (cols_ != m.cols_) + return false; + return isEqual(m); +} + +bool Matrix::operator!=(Matrix m) { + return !(*this == m); +} + +const double& Matrix::operator()(int row, int col) const { + return GET(*this, row, col); +} + +double& Matrix::operator()(int row, int col) { + return GET(*this, row, col); +} + +Matrix Matrix::t() { + Matrix m(*this); + m.rows_ = cols_; + m.cols_ = rows_; + for (int i = 0; i < rows_; i++) + for (int j = 0; j < cols_; j++) + GET(m, i, j) = GET(*this, j, i); + return m; +} + +Matrix Matrix::operator+(double value) { + Matrix m(*this); + m.operateScalar(value, SUM); + return m; +} + +Matrix Matrix::operator-(double value) { + Matrix m(*this); + m.operateScalar(value, DIFF); + return m; +} + +Matrix Matrix::operator*(double value) { + Matrix m(*this); + m.operateScalar(value, PROD); + return m; +} + +Matrix Matrix::operator/(double value) { + Matrix m(*this); + m.operateScalar(value, DIV); + return m; +} + +Matrix& Matrix::operator+=(double value) { + operateScalar(value, SUM); + return *this; +} + +Matrix& Matrix::operator-=(double value) { + operateScalar(value, DIFF); + return *this; +} + +Matrix& Matrix::operator*=(double value) { + operateScalar(value, PROD); + return *this; +} + +Matrix& Matrix::operator/=(double value) { + operateScalar(value, DIV); + return *this; +} + +Matrix Matrix::operator+(Matrix &m) { + Matrix n(*this); + n.operateBlock(m, SUM); + return n; +} + +Matrix Matrix::operator-(Matrix &m) { + Matrix n(*this); + n.operateBlock(m, DIFF); + return n; +} + +Matrix& Matrix::operator+=(Matrix &m) { + operateBlock(m, SUM); + return *this; +} + +Matrix& Matrix::operator-=(Matrix &m) { + operateBlock(m, DIFF); + return *this; +} + +Matrix Matrix::operator*(Matrix& m) { + double sum; + Matrix n(rows_, m.cols_); + for (int i = 0; i < n.rows_; i++) + for (int j = 0; j < n.cols_; j++) { + sum = 0; + for (int k = 0; k < cols_; k++) + sum += GET(*this, i, k) * GET(m, k, j); + GET(n, i, j) = sum; + } + return n; +} diff --git a/linalg/matrix_io.cc b/linalg/matrix_io.cc new file mode 100644 index 0000000..f8732c9 --- /dev/null +++ b/linalg/matrix_io.cc @@ -0,0 +1,53 @@ +#include +#include "linalg.h" + +void Matrix::read(std::ifstream& ifs) { + ifs.read(reinterpret_cast(&rows_), sizeof(int)); + ifs.read(reinterpret_cast(&cols_), sizeof(int)); + int size__new = rows_ * cols_; + if (size_ != size__new) { + if (data_ != nullptr) + free(data_); + data_ = static_cast(malloc(sizeof(double) * size__new)); + } + size_ = size__new; + ifs.read(reinterpret_cast(data_), size_ * sizeof(double)); +} + +void Matrix::fread(const char *filename) { + std::ifstream ifs(filename, std::ios::binary); + read(ifs); + ifs.close(); +} + +void Matrix::fwrite(const char *filename) { + std::ofstream ofs(filename, std::ios::out); + writetxt(ofs); + ofs.close(); +} + +void Matrix::write(std::ofstream& ofs) { + ofs.write(reinterpret_cast(&rows_), sizeof(int)); + ofs.write(reinterpret_cast(&cols_), sizeof(int)); + ofs.write(reinterpret_cast(data_), size_ * sizeof(double)); +} + +void Matrix::writetxt(std::ofstream& ofs) { + for (int i = 0; i < rows_; i++) + for (int j = 0; j < cols_; j++) { + ofs << GET(*this, i, j); + ofs << (j < cols_ - 1 ? ' ' : '\n'); + } +} + +std::ostream& operator<<(std::ostream& os, Matrix& m) { + for (int i = 0; i < m.rows(); i++) { + for (int j = 0; j < m.cols(); j++) { + std::cout << m(i, j); + if (j < m.cols() - 1) + std::cout << ' '; + } + std::cout << std::endl; + } + return os; +} diff --git a/linalg/test.dat b/linalg/test.dat new file mode 100644 index 0000000..2c3f803 --- /dev/null +++ b/linalg/test.dat @@ -0,0 +1,2 @@ +1 2 +3 4 diff --git a/linalg/vector.cc b/linalg/vector.cc new file mode 100644 index 0000000..3af9079 --- /dev/null +++ b/linalg/vector.cc @@ -0,0 +1,155 @@ +#include "linalg.h" + +Vector Vector::operator+(double value) { + Vector m(*this); + m.operateScalar(value, SUM); + return m; +} + +Vector Vector::operator-(double value) { + Vector m(*this); + m.operateScalar(value, DIFF); + return m; +} + +Vector Vector::operator*(double value) { + Vector m(*this); + m.operateScalar(value, PROD); + return m; +} + +Vector Vector::operator/(double value) { + Vector m(*this); + m.operateScalar(value, DIV); + return m; +} + +Vector& Vector::operator+=(double value) { + operateScalar(value, SUM); + return *this; +} + +Vector& Vector::operator-=(double value) { + operateScalar(value, DIFF); + return *this; +} + +Vector& Vector::operator*=(double value) { + operateScalar(value, PROD); + return *this; +} + +Vector& Vector::operator/=(double value) { + operateScalar(value, DIV); + return *this; +} + +Vector Vector::operator+(Vector &m) { + Vector n(*this); + n.operateBlock(m, SUM); + return n; +} + +Vector Vector::operator-(Vector &m) { + Vector n(*this); + n.operateBlock(m, DIFF); + return n; +} + +Vector& Vector::operator+=(Vector &m) { + operateBlock(m, SUM); + return *this; +} + +Vector& Vector::operator-=(Vector &m) { + operateBlock(m, DIFF); + return *this; +} + +Vector::Vector() : Block() { + coltype_ = true; +} + +Vector::Vector(int size) : Block(size) { + coltype_ = true; +} + +Vector::Vector(const Vector& m) : Block(m) { + coltype_ = m.coltype_; +} + +Vector::Vector(std::initializer_list list) + : Vector(list.size()) { + *this = list; +} + +Vector& Vector::operator=(Vector m) { + if (this == &m) + return *this; + coltype_ = m.coltype_; + assign(m); + return *this; +} + +Vector& Vector::operator=(std::initializer_list list) { + int size_new = list.size(); + if (size_ != size_new) { + if (data_ != nullptr) + delete[] data_; + data_ = new double[size_new]; + } + size_ = size_new; + int i = 0; + for (double value : list) + data_[i++] = value; + return *this; +} + +bool Vector::operator==(Vector x) { + if (coltype_ != x.coltype_) + return false; + return isEqual(x); +} + +bool Vector::operator!=(Vector m) { + return !(*this == m); +} + +const double& Vector::operator()(int index) const { + return data_[index]; +} + +double& Vector::operator()(int index) { + return data_[index]; +} + +void Vector::write(std::ofstream& ofs) { + Matrix m; + if (coltype_) { + m.cols_ = 1; + m.rows_ = this->size_; + } + else { + m.rows_ = 1; + m.cols_ = this->size_; + } + m.size_ = this->size_; + m.data_ = this->data_; + m.write(ofs); + m.data_ = nullptr; + +} + +Vector Vector::t() { + Vector x(*this); + x.coltype_ = !(this->coltype_); + return x; +} + + +double Vector::operator*(Vector& x) { + double sum = 0; + for (int i = 0; i < size_; i++) + sum += data_[i] * x.data_[i]; + return sum; +} diff --git a/linalg/vector_io.cc b/linalg/vector_io.cc new file mode 100644 index 0000000..b4f43f2 --- /dev/null +++ b/linalg/vector_io.cc @@ -0,0 +1,29 @@ +#include +#include "linalg.h" + +using namespace std; + +void Vector::read(ifstream& ifs) { + Matrix m; + if (coltype_) { + m.cols_ = 1; + m.rows_ = this->size_; + } + else { + m.rows_ = 1; + m.cols_ = this->size_; + } + m.size_ = this->size_; + m.data_ = this->data_; + m.read(ifs); + m.data_ = nullptr; + +} + +ostream& operator<<(ostream& os, Vector& m) { + for (int i = 0; i < m.size(); i++) { + cout << m(i); + cout << (i < m.size() - 1 ? ' ' : '\n'); + } + return os; +}