Vectors & Matrices

This commit is contained in:
Jaime Lopez 2021-11-25 13:55:46 -05:00
parent d8de858ecc
commit a6bc879c83
10 changed files with 822 additions and 0 deletions

20
linalg/Makefile Normal file
View File

@ -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

3
linalg/README Normal file
View File

@ -0,0 +1,3 @@
# Vectors & Matrices
A toy implementation of vector and matrix classes for linear algebra operations.

78
linalg/block.cc Normal file
View File

@ -0,0 +1,78 @@
#include <cstdlib>
#include <cstring>
#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;
}
}

115
linalg/linalg.h Normal file
View File

@ -0,0 +1,115 @@
#ifndef _LINALG_H
#define _LINALG_H 1
#include <fstream>
#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<std::initializer_list<double>> list);
Matrix(const char *filename);
Matrix& operator=(Matrix m);
Matrix& operator=(std::initializer_list<std::initializer_list<double>> 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<double> list);
Vector(const char *filename);
Vector& operator=(Vector m);
Vector& operator=(std::initializer_list<double> 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

210
linalg/linalg_test.cc Normal file
View File

@ -0,0 +1,210 @@
#include "linalg.h"
#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE matrix
#include <boost/test/unit_test.hpp>
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);
}

157
linalg/matrix.cc Normal file
View File

@ -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<std::initializer_list<double>> 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<std::initializer_list<double>> 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<double> 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;
}

53
linalg/matrix_io.cc Normal file
View File

@ -0,0 +1,53 @@
#include <iostream>
#include "linalg.h"
void Matrix::read(std::ifstream& ifs) {
ifs.read(reinterpret_cast<char *>(&rows_), sizeof(int));
ifs.read(reinterpret_cast<char *>(&cols_), sizeof(int));
int size__new = rows_ * cols_;
if (size_ != size__new) {
if (data_ != nullptr)
free(data_);
data_ = static_cast<double*>(malloc(sizeof(double) * size__new));
}
size_ = size__new;
ifs.read(reinterpret_cast<char *>(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<char *>(&rows_), sizeof(int));
ofs.write(reinterpret_cast<char *>(&cols_), sizeof(int));
ofs.write(reinterpret_cast<char *>(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;
}

2
linalg/test.dat Normal file
View File

@ -0,0 +1,2 @@
1 2
3 4

155
linalg/vector.cc Normal file
View File

@ -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<double> 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<double> 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;
}

29
linalg/vector_io.cc Normal file
View File

@ -0,0 +1,29 @@
#include <iostream>
#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;
}