Matrix.h

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2018-2025 ISciences, LLC
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation.
 * See the COPYING file for more information.
 *
 **********************************************************************/
 
#pragma once
 
#include <cstring>
#include <iomanip>
#include <memory>
#include <string>
#include <vector>
 
namespace geos::operation::grid {
 
template<typename T>
class Matrix
{
 
  public:
    using value_type = T;
 
    Matrix(size_t rows, size_t cols, std::shared_ptr<T[]> data)
        : m_data{ data }
        , m_rows{ rows }
        , m_cols{ cols }
    {}
 
    Matrix(size_t rows, size_t cols)
      : m_rows{ rows }
      , m_cols{ cols }
    {
        if (m_rows > 0 && m_cols > 0) {
            // new T[]() initializes to zero
            m_data = std::shared_ptr<T[]>(new T[m_rows * m_cols]());
        }
    }
 
    Matrix(size_t rows, size_t cols, T value)
      : m_rows{ rows }
      , m_cols{ cols }
    {
        if (m_rows > 0 && m_cols > 0) {
            // new T[] does not initialize
            m_data = std::shared_ptr<T[]>(new T[m_rows * m_cols]);
        }
 
        std::fill(m_data.get(), m_data.get() + m_rows * m_cols, value);
    }
 
    explicit Matrix(const std::vector<std::vector<T>>& data)
      : m_rows{ data.size() }
      , m_cols{ data[0].size() }
    {
        m_data = std::shared_ptr<T[]>(new T[m_rows * m_cols]());
 
        auto lastpos = m_data.get();
        for (auto& row : data) {
            lastpos = std::copy(row.begin(), row.end(), lastpos);
        }
    }
 
    Matrix(Matrix<T>&& m) noexcept
      : m_rows{ m.getNumRows() }
      , m_cols{ m.getNumCols() }
    {
        m_data = std::move(m.m_data);
    }
 
    T& operator()(size_t row, size_t col)
    {
        check(row, col);
        return m_data.get()[row * m_cols + col];
    }
 
    const T& operator()(size_t row, size_t col) const
    {
        check(row, col);
        return m_data.get()[row * m_cols + col];
    }
 
    bool operator!=(const Matrix<T>& other) const {
        return !(*this == other);
    }
 
    bool operator==(const Matrix<T>& other) const
    {
        if (m_rows != other.m_rows) {
            return false;
        }
        if (m_cols != other.m_cols) {
            return false;
        }
 
        return 0 == memcmp(m_data.get(), other.m_data.get(), m_rows * m_cols * sizeof(T));
    }
 
    void increment(size_t row, size_t col, const T& val)
    {
        check(row, col);
        m_data.get()[row * m_cols + col] += val;
    }
 
    size_t getNumRows() const { return m_rows; }
    size_t getNumCols() const { return m_cols; }
 
    T* row(size_t row)
    {
        return &(m_data[row * m_cols]);
    }
 
    T* data()
    {
        return m_data.get();
    }
 
    const T* data() const
    {
        return m_data.get();
    }
 
    T* release()
    {
        return m_data.release();
    }
 
    const T* begin() const {
        return m_data.get();
    }
 
    const T* end() const {
        return m_data.get() + m_rows * m_cols;
    }
 
#ifdef MATRIX_CHECK_BOUNDS
    void check(size_t row, size_t col) const
    {
        if (row + 1 > m_rows) {
            throw std::out_of_range("Row " + std::to_string(row) + " is out of range.");
        }
        if (col + 1 > m_cols) {
            throw std::out_of_range("Col " + std::to_string(col) + " is out of range.");
        }
    }
#else
    void check(size_t, size_t) const
    {
    }
#endif
 
  private:
    std::shared_ptr<T[]> m_data;
 
    size_t m_rows;
    size_t m_cols;
};
 
template<typename T>
std::ostream&
operator<<(std::ostream& os, const Matrix<T>& m)
{
    for (size_t i = 0; i < m.getNumRows(); i++) {
        for (size_t j = 0; j < m.getNumCols(); j++) {
            if (m(i, j) != 0) {
                os << std::right << std::fixed << std::setw(10) << std::setprecision(6) << m(i, j) << " ";
            } else {
                os << "           ";
            }
        }
        os << std::endl;
    }
 
    return os;
}
 
}