ell_matrix.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cstddef>
#include <valarray>
#include <vector>
 
class SliceView {
    double* const _data;
    const std::size_t _stride;
    const std::size_t _len;
 
  public:
    SliceView(double* data, std::size_t stride, std::size_t len) noexcept
        : _data{data}, _stride{stride}, _len{len} {}
 
    double& operator[](std::size_t i) { return this->_data[i * this->_stride]; }
    const double& operator[](std::size_t i) const { return this->_data[i * this->_stride]; }
 
    SliceView& operator=(const std::valarray<double>& v) {
        for (std::size_t i = 0; i < this->_len; ++i) this->_data[i * this->_stride] = v[i];
        return *this;
    }
 
    SliceView& operator=(double scalar) {
        for (std::size_t i = 0; i < this->_len; ++i) this->_data[i * this->_stride] = scalar;
        return *this;
    }
 
    operator std::valarray<double>() const {
        std::valarray<double> res(this->_len);
        for (std::size_t i = 0; i < this->_len; ++i) res[i] = this->_data[i * this->_stride];
        return res;
    }
 
    double sum() const {
        double s = 0.0;
        for (std::size_t i = 0; i < this->_len; ++i) s += this->_data[i * this->_stride];
        return s;
    }
};
 
class Matrix {
    std::size_t ndim;
    std::vector<double> data;
 
  public:
    explicit Matrix(std::size_t ndim, double init = 0.0) : ndim{ndim}, data(ndim * ndim, init) {}
 
    void clear(double value = 0.0) { std::ranges::fill(this->data, value); }
 
    double& operator()(std::size_t row, std::size_t col) {
        return this->data[row * this->ndim + col];
    }
 
    const double& operator()(std::size_t row, std::size_t col) const {
        return this->data[row * this->ndim + col];
    }
 
    Matrix& operator*=(double alpha) {
        for (auto& v : this->data) v *= alpha;
        return *this;
    }
 
    Matrix operator*(double alpha) const {
        Matrix res(*this);
        return res *= alpha;
    }
 
    void identity() {
        this->clear();
        auto d = this->diagonal();
        for (std::size_t i = 0; i < this->ndim; ++i) d[i] = 1.0;
    }
 
    SliceView diagonal() { return SliceView(this->data.data(), this->ndim + 1, this->ndim); }
 
    SliceView row(std::size_t r) {
        return SliceView(this->data.data() + r * this->ndim, 1, this->ndim);
    }
 
    SliceView column(std::size_t c) {
        return SliceView(this->data.data() + c, this->ndim, this->ndim);
    }
 
    double trace() const {
        double s = 0.0;
        for (std::size_t i = 0; i < this->ndim; ++i) {
            s += this->operator()(i, i);
        }
        return s;
    }
};