MNP02/matrix.h

#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) const;
    self_t operator+(const self_t& rhs) const;

    self_t operator*(T rhs) const;
    self_t operator+(T rhs) const;

    self_t operator-() const;
    self_t operator-(const self_t& rhs) const;

    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) const {
    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) const {
    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 {
    self_t result(*this);
    return result * -1;
}

template<typename T>
Matrix<T> Matrix<T>::operator-(const Matrix<T>& rhs) const {
    self_t result(*this);

    for (size_t i = 0; i < _rows; i++)
        for (size_t j = 0; j < _cols; j++)
            result(i, j) -= rhs(i, j);

    return result;
}

template<typename T>
Matrix<T> Matrix<T>::operator*(const T rhs) const {
    self_t result(*this);

    for (size_t i = 0; i < _rows; i++)
        for (size_t j = 0; j < _cols; j++)
            result(i, j) *= rhs(i, j);

    return result;
}

template<typename T>
Matrix<T> Matrix<T>::operator+(const T rhs) const {
    self_t result(*this);

    for (size_t i = 0; i < _rows; i++)
        for (size_t j = 0; j < _cols; j++)
            result(i, j) += rhs(i, j);

    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