matrix utility

.gitignore

@@ -0,0 +1,50 @@
+/.idea/
+
+
+# common
+*.csv
+*.obj
+*.exe
+*.lst
+*.tmp
+*.pdb
+
+# tex
+*.aux
+*.lof
+*.log
+*.lot
+*.fls
+*.out
+*.toc
+*.fmt
+*.fot
+*.cb
+*.cb2
+.*.lb
+*.pdf
+
+/cmake-*/
+
+/sprawozdanie/*-figure*
+
+## Build tool auxiliary files:
+*.fdb_latexmk
+*.synctex
+*.synctex(busy)
+*.synctex.gz
+*.synctex.gz(busy)
+*.pdfsync
+*.auxlock
+
+# hyperref
+*.brf
+
+# listings
+*.lol
+
+# makeidx
+*.idx
+*.ilg
+*.ind
+*.ist

CMakeLists.txt

@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 3.10)
+project(P02)
+
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_C_STANDARD 11)
+
+add_executable(main main.cpp)

main.cpp

@@ -0,0 +1,29 @@
+#include "matrix.h"
+#include <iostream>
+#include <utility>
+
+const size_t N = 10;
+
+std::pair<Matrix<double>, Matrix<double>> prepare(size_t index, size_t n) {
+    double a1, a2, a3;
+    a2 = a3 = -1;
+    a1 = 5 + (index / 100) % 10;
+
+    auto M = Matrix<double>::diag(n, a1)
+           + Matrix<double>::diag(n, a2, 1) + Matrix<double>::diag(n, a2, -1)
+           + Matrix<double>::diag(n, a3, 2) + Matrix<double>::diag(n, a3, -2);
+
+    Matrix<double> b(n, 1);
+    for (size_t i = 0; i < n; ++i) {
+        b(i, 0) = sin((i + 1)*((index / 1000) % 10 + 1));
+    }
+
+    return std::make_pair(M, b);
+}
+
+int main(int argc, const char argv[])
+{
+    auto tuple = prepare(165581, N);
+    auto M = std::get<0>(tuple);
+    auto b = std::get<1>(tuple);
+}

matrix.h

@@ -0,0 +1,233 @@
+#ifndef P02_MATRIX_H
+#define P02_MATRIX_H
+
+#include <memory>
+#include <vector>
+#include <cassert>
+#include <ostream>
+#include <iostream>
+#include <utility>
+
+#define NORM_INF INT_MAX
+
+using std::size_t;
+
+template <typename T>
+class Matrix {
+    using self_t = Matrix<T>;
+
+    std::unique_ptr<T[]> values;
+    size_t _rows, _cols;
+
+public:
+    explicit Matrix(size_t n);
+    Matrix(std::initializer_list<std::initializer_list<T>> values);
+    Matrix(size_t n, size_t m);
+
+    Matrix(const self_t& m);
+
+    T& operator()(size_t n, size_t m);
+    T operator()(size_t n, size_t m) const;
+
+    self_t operator*(const self_t& rhs);
+    self_t operator+(const self_t& rhs);
+
+    self_t operator*(T rhs);
+    self_t operator+(T rhs);
+
+    self_t& operator=(self_t&& rhs) noexcept;
+
+    size_t rows() const;
+    size_t cols() const;
+
+    template <typename U>
+    friend std::ostream &operator<<(std::ostream &os, const Matrix<U> &matrix);
+
+    static self_t diag(size_t n, T value, int offset = 0);
+    static self_t hvec(std::initializer_list<T> values);
+    static self_t vvec(std::initializer_list<T> values);
+//    static self_t diag(std::initializer_list values, size_t offset = 0);
+};
+
+template<typename T>
+Matrix<T>::Matrix(size_t n, size_t m) : _rows(n), _cols(m), values(new T[m*n]) {
+    for(size_t i = 0; i < m*n; i++) {
+        this->values.get()[i] = T();
+    }
+}
+
+template<typename T>
+Matrix<T>::Matrix(size_t n) : Matrix(n, n) { }
+
+template<typename T>
+Matrix<T>::Matrix(const Matrix<T> &m) : Matrix(m._rows, m._cols) {
+    std::copy(m.values.get(), m.values.get() + m._cols*m._rows, values.get());
+}
+
+template<typename T>
+T &Matrix<T>::operator()(size_t n, size_t m) {
+    return values.get()[n*this->_cols + m];
+}
+
+template<typename T>
+T Matrix<T>::operator()(size_t n, size_t m) const {
+    return values.get()[n*this->_cols + m];
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::operator+(const Matrix<T> &rhs) {
+    assert(this->_cols == rhs._cols && this->_rows == rhs._rows);
+    self_t result(this->_rows, this->_cols);
+
+    for (size_t i = 0; i < _rows; i++)
+        for (size_t j = 0; j < _cols; j++)
+            result(i, j) = (*this)(i, j) + rhs(i, j);
+
+    return result;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::operator*(const Matrix<T> &rhs) {
+    assert(_cols == rhs._rows);
+    self_t result(_rows, rhs._cols);
+
+    for (size_t i = 0; i < _rows; i++) {
+        for (size_t j = 0; j < rhs._cols; j++) {
+            T accumulator = 0;
+            for(size_t k = 0; k < _cols; k++)
+                accumulator += (*this)(i, k) * rhs(k, j);
+
+            result(i, j) = std::move(accumulator);
+        }
+    }
+
+    return result;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::operator*(const T rhs) {
+    self_t result(*this);
+
+    for (size_t i = 0; i < _rows; i++)
+        for (size_t j = 0; j < _cols; j++)
+            result(i, j) *= rhs;
+
+    return result;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::operator+(const T rhs) {
+    self_t result(*this);
+
+    for (size_t i = 0; i < _rows; i++)
+        for (size_t j = 0; j < _cols; j++)
+            result(i, j) += rhs;
+
+    return result;
+}
+
+template<typename T>
+std::ostream &operator<<(std::ostream &os, const Matrix<T> &matrix) {
+    os << "[ ";
+
+    for (size_t i = 0; i < matrix._rows; i++) {
+        for (size_t j = 0; j < matrix._cols; j++) {
+            os << matrix(i, j) << " ";
+        }
+        os << (i == matrix._rows - 1 ? "]\n" : "\n  ");
+    }
+
+    return os;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::diag(size_t n, T value, int offset) {
+    self_t result(n);
+    if (offset >= 0)
+        for(size_t i = 0; i < n - offset; i++)
+            result(i,i+offset) = value;
+    else
+        for(size_t i = 0; i < n - offset; i++)
+            result(i-offset, i) = value;
+
+    return result;
+}
+
+template<typename T>
+Matrix<T>::Matrix(std::initializer_list<std::initializer_list<T>> values) : Matrix(values.size(), values.begin()->size()) {
+    size_t i = 0;
+    for (auto h : values) {
+        if (i > _rows) break;
+
+        size_t j = 0;
+        for (auto x : h) {
+            if (j > _cols) break;
+
+            (*this)(i, j) = x;
+            j++;
+        }
+        i++;
+    }
+}
+
+template<typename T>
+Matrix<T>& Matrix<T>::operator=(Matrix<T> &&rhs) noexcept {
+    this->_cols = rhs._cols;
+    this->_rows = rhs._rows;
+    std::swap(values, rhs.values);
+
+    return *this;
+}
+
+template<typename T>
+inline size_t Matrix<T>::rows() const {
+    return this->_rows;
+}
+
+template<typename T>
+inline size_t Matrix<T>::cols() const {
+    return this->_cols;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::hvec(std::initializer_list<T> values) {
+    Matrix<T> result(1, values.size());
+
+    size_t i = 0;
+    for(auto val : values) {
+        result(0, i++) = val;
+    }
+
+    return result;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::vvec(std::initializer_list<T> values) {
+    Matrix<T> result(values.size(), 1);
+
+    size_t i = 0;
+    for(auto val : values) {
+        result(i++, 0) = val;
+    }
+
+    return result;
+}
+
+template <typename T>
+double norm(const Matrix<T>& matrix, unsigned norm = 2)
+{
+//    needed?
+//    assert(matrix.cols() == 1 || matrix.rows() == 1);
+
+    double accumulator = 0.;
+
+    for (size_t i = 0; i < matrix.cols(); ++i) {
+        for (size_t j = 0; j < matrix.rows(); ++j) {
+            accumulator += std::pow(matrix(j, i), norm);
+        }
+    }
+
+    return std::pow(accumulator, 1. / norm);
+}
+
+#endif //P02_MATRIX_H