|
|
@@ -0,0 +1,804 @@ |
|
|
|
/** |
|
|
|
* \file matrix.hpp |
|
|
|
* \brief A matrix abstraction implementation |
|
|
|
* |
|
|
|
* \author |
|
|
|
* Christos Choutouridis AEM:8997 |
|
|
|
* <cchoutou@ece.auth.gr> |
|
|
|
*/ |
|
|
|
#ifndef MATRIX_HPP_ |
|
|
|
#define MATRIX_HPP_ |
|
|
|
|
|
|
|
#include <type_traits> |
|
|
|
#include <utility> |
|
|
|
#include <algorithm> |
|
|
|
#include <vector> |
|
|
|
#include <tuple> |
|
|
|
|
|
|
|
namespace mtx { |
|
|
|
|
|
|
|
using std::size_t; |
|
|
|
|
|
|
|
/* |
|
|
|
* Small helper to strip types |
|
|
|
*/ |
|
|
|
template<typename T> |
|
|
|
struct remove_cvref { |
|
|
|
typedef std::remove_cv_t<std::remove_reference_t<T>> type; |
|
|
|
}; |
|
|
|
template<typename T> |
|
|
|
using remove_cvref_t = typename remove_cvref<T>::type; |
|
|
|
|
|
|
|
/*! |
|
|
|
* Enumerator to denote the storage type of the array to use. |
|
|
|
*/ |
|
|
|
enum class MatrixType { |
|
|
|
DENSE, /*!< Matrix is dense */ |
|
|
|
SPARSE, /*!< Matrix is sparse */ |
|
|
|
}; |
|
|
|
|
|
|
|
/*! |
|
|
|
* Enumerator to denote the storage type of the array to use. |
|
|
|
*/ |
|
|
|
enum class MatrixOrder { |
|
|
|
COLMAJOR, /*!< Matrix is column major */ |
|
|
|
ROWMAJOR, /*!< Matrix is row major */ |
|
|
|
}; |
|
|
|
|
|
|
|
/* |
|
|
|
* Forward type declarations |
|
|
|
*/ |
|
|
|
|
|
|
|
template<typename MatrixType> struct MatCol; |
|
|
|
template<typename MatrixType> struct MatRow; |
|
|
|
template<typename MatrixType> struct MatVal; |
|
|
|
|
|
|
|
/*! |
|
|
|
* A 2-D matrix functionality over a 1-D array |
|
|
|
* |
|
|
|
* This is a very thin abstraction layer over a native array. |
|
|
|
* This is tested using compiler explorer and our template produce |
|
|
|
* almost identical assembly. |
|
|
|
* |
|
|
|
* The penalty hit we have is due to the fact that we use a one dimension array |
|
|
|
* and we have to calculate the actual position from an (i,j) pair. |
|
|
|
* The use of 1D array was our intention from the beginning, so the penalty |
|
|
|
* was pretty much unavoidable. |
|
|
|
* |
|
|
|
* \tparam DataType The underling data type of the array |
|
|
|
* \tparam IndexType The underling type for the index variables and sizes |
|
|
|
* \tparam Type The storage type of the array |
|
|
|
* \arg FULL For full matrix |
|
|
|
* \arg SYMMETRIC For symmetric matrix (we use only the lower part) |
|
|
|
*/ |
|
|
|
template<typename DataType, |
|
|
|
typename IndexType = size_t, |
|
|
|
MatrixType Type = MatrixType::DENSE, |
|
|
|
MatrixOrder Order = MatrixOrder::ROWMAJOR, |
|
|
|
bool Symmetric = false> |
|
|
|
struct Matrix { |
|
|
|
|
|
|
|
using dataType = DataType; //!< meta:export of underling data type |
|
|
|
using indexType = IndexType; //!< meta:export of underling index type |
|
|
|
static constexpr MatrixOrder matrixOrder = Order; //!< meta:export of array order |
|
|
|
static constexpr MatrixType matrixType = Type; //!< meta:export of array type |
|
|
|
static constexpr bool symmetric = Symmetric; //!< meta:export symmetric flag |
|
|
|
|
|
|
|
/*! |
|
|
|
* \name Obj lifetime |
|
|
|
*/ |
|
|
|
//! @{ |
|
|
|
|
|
|
|
//! Construct an empty matrix with dimensions rows x columns |
|
|
|
Matrix(IndexType rows = IndexType{}, IndexType columns = IndexType{}) noexcept |
|
|
|
: vector_storage_(capacity(rows, columns)), |
|
|
|
raw_storage_(nullptr), |
|
|
|
use_vector_(true), |
|
|
|
rows_(rows), |
|
|
|
cols_(columns) { |
|
|
|
data_ = vector_storage_.data(); |
|
|
|
} |
|
|
|
|
|
|
|
//! Construct a matrix by copying existing data with dimensions rows x columns |
|
|
|
Matrix(DataType* data, IndexType major_start, IndexType major_length, IndexType minor_length) noexcept |
|
|
|
: vector_storage_(), |
|
|
|
raw_storage_ (data + major_start * minor_length), |
|
|
|
use_vector_ (false) { |
|
|
|
if constexpr (Order == MatrixOrder::ROWMAJOR) { |
|
|
|
rows_ = major_length; |
|
|
|
cols_ = minor_length; |
|
|
|
} |
|
|
|
else { |
|
|
|
rows_ = minor_length; |
|
|
|
cols_ = major_length; |
|
|
|
} |
|
|
|
data_ = raw_storage_; |
|
|
|
} |
|
|
|
|
|
|
|
//! Construct a matrix using an initializer list |
|
|
|
Matrix(IndexType rows, IndexType columns, std::initializer_list<DataType> list) |
|
|
|
: vector_storage_(list), |
|
|
|
raw_storage_(nullptr), |
|
|
|
use_vector_(true), |
|
|
|
rows_(rows), |
|
|
|
cols_(columns) { |
|
|
|
if (list.size() != capacity(rows, columns)) { |
|
|
|
throw std::invalid_argument("Matrix initializer list size does not match matrix dimensions."); |
|
|
|
} |
|
|
|
data_ = vector_storage_.data(); |
|
|
|
} |
|
|
|
|
|
|
|
//! move ctor |
|
|
|
Matrix(Matrix&& m) noexcept { moves(std::move(m)); } |
|
|
|
//! move |
|
|
|
Matrix& operator=(Matrix&& m) noexcept { moves(std::move(m)); return *this; } |
|
|
|
Matrix(const Matrix& m) = delete; //!< No copy ctor |
|
|
|
Matrix& operator=(const Matrix& m) = delete; //!< No copy |
|
|
|
//Matrix(const Matrix& m); |
|
|
|
//Matrix& operator=(const Matrix& m) { copy(m); } |
|
|
|
|
|
|
|
//! @} |
|
|
|
|
|
|
|
//! \name Data exposure |
|
|
|
//! @{ |
|
|
|
|
|
|
|
|
|
|
|
//! Get/Set the size of each dimension |
|
|
|
IndexType rows() const noexcept { return rows_; } |
|
|
|
IndexType columns() const noexcept { return cols_; } |
|
|
|
|
|
|
|
//! Get the interface size of the Matrix (what appears to be the size) |
|
|
|
IndexType size() const { |
|
|
|
return rows_ * cols_; |
|
|
|
} |
|
|
|
//! Set the interface size of the Matrix (what appears to be the size) |
|
|
|
IndexType resize(IndexType rows, IndexType columns) { |
|
|
|
if (use_vector_) { |
|
|
|
rows_ = rows; |
|
|
|
cols_ = columns; |
|
|
|
vector_storage_.resize(capacity(rows_, cols_)); |
|
|
|
data_ = vector_storage_.data(); |
|
|
|
} |
|
|
|
return capacity(rows_, cols_); |
|
|
|
} |
|
|
|
|
|
|
|
//! Actual memory capacity of the symmetric matrix |
|
|
|
static constexpr IndexType capacity(IndexType M, IndexType N) { |
|
|
|
if constexpr (Symmetric) |
|
|
|
return (M+1)*N/2; |
|
|
|
else |
|
|
|
return M*N; |
|
|
|
} |
|
|
|
|
|
|
|
/* |
|
|
|
* virtual 2D accessors |
|
|
|
*/ |
|
|
|
DataType get (IndexType i, IndexType j) { |
|
|
|
if constexpr (Symmetric) { |
|
|
|
auto T = [](size_t i)->size_t { return i*(i+1)/2; }; // Triangular number of i |
|
|
|
if constexpr (Order == MatrixOrder::COLMAJOR) { |
|
|
|
// In column major we use the lower triangle of the matrix |
|
|
|
if (i>=j) return data_[j*rows_ - T(j) + i]; // Lower, use our notation |
|
|
|
else return data_[i*rows_ - T(i) + j]; // Upper, use opposite index |
|
|
|
} |
|
|
|
else { |
|
|
|
// In row major we use the upper triangle of the matrix |
|
|
|
if (i<=j) return data_[i*cols_ - T(i) + j]; // Upper, use our notation |
|
|
|
else return data_[j*cols_ - T(j) + i]; // Lower, use opposite index |
|
|
|
} |
|
|
|
} |
|
|
|
else { |
|
|
|
if constexpr (Order == MatrixOrder::COLMAJOR) |
|
|
|
return data_[i + j*rows_]; |
|
|
|
else |
|
|
|
return data_[i*cols_ + j]; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
/*! |
|
|
|
* \fn DataType set(DataType, IndexType, IndexType) |
|
|
|
* \param v |
|
|
|
* \param i |
|
|
|
* \param j |
|
|
|
* \return |
|
|
|
*/ |
|
|
|
DataType set (DataType v, IndexType i, IndexType j) { |
|
|
|
if constexpr (Symmetric) { |
|
|
|
auto T = [](size_t i)->size_t { return i*(i+1)/2; }; // Triangular number of i |
|
|
|
if constexpr (Order == MatrixOrder::COLMAJOR) { |
|
|
|
// In column major we use the lower triangle of the matrix |
|
|
|
if (i>=j) return data_[j*rows_ - T(j) + i] = v; // Lower, use our notation |
|
|
|
else return data_[i*rows_ - T(i) + j] = v; // Upper, use opposite index |
|
|
|
} |
|
|
|
else { |
|
|
|
// In row major we use the upper triangle of the matrix |
|
|
|
if (i<=j) return data_[i*cols_ - T(i) + j] = v; // Upper, use our notation |
|
|
|
else return data_[j*cols_ - T(j) + i] = v; // Lower, use opposite index |
|
|
|
} |
|
|
|
} |
|
|
|
else { |
|
|
|
if constexpr (Order == MatrixOrder::COLMAJOR) |
|
|
|
return data_[i + j*rows_] = v; |
|
|
|
else |
|
|
|
return data_[i*cols_ + j] = v; |
|
|
|
} |
|
|
|
} |
|
|
|
// DataType operator()(IndexType i, IndexType j) { return get(i, j); } |
|
|
|
/*! |
|
|
|
* Return a proxy MatVal object with read and write capabilities. |
|
|
|
* @param i The row number |
|
|
|
* @param j The column number |
|
|
|
* @return tHE MatVal object |
|
|
|
*/ |
|
|
|
MatVal<Matrix> operator()(IndexType i, IndexType j) noexcept { |
|
|
|
return MatVal<Matrix>(this, get(i, j), i, j); |
|
|
|
} |
|
|
|
|
|
|
|
// a basic serial iterator support |
|
|
|
DataType* data() noexcept { return data_; } |
|
|
|
DataType* begin() noexcept { return data_; } |
|
|
|
const DataType* begin() const noexcept { return data_; } |
|
|
|
DataType* end() noexcept { return data_ + capacity(rows_, cols_); } |
|
|
|
const DataType* end() const noexcept { return data_ + capacity(rows_, cols_); } |
|
|
|
|
|
|
|
// IndexType begin_idx() noexcept { return 0; } |
|
|
|
// IndexType end_idx() noexcept { return capacity(rows_, cols_); } |
|
|
|
|
|
|
|
const DataType* data() const noexcept { return data_; } |
|
|
|
const IndexType begin_idx() const noexcept { return 0; } |
|
|
|
const IndexType end_idx() const noexcept { return capacity(rows_, cols_); } |
|
|
|
//! @} |
|
|
|
|
|
|
|
/*! |
|
|
|
* \name Safe iteration API |
|
|
|
* |
|
|
|
* This api automates the iteration over the array based on |
|
|
|
* MatrixType |
|
|
|
*/ |
|
|
|
//! @{ |
|
|
|
template<typename F, typename... Args> |
|
|
|
void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) { |
|
|
|
for (IndexType it=begin ; it<end ; ++it) { |
|
|
|
std::forward<F>(lambda)(std::forward<Args>(args)..., it); |
|
|
|
} |
|
|
|
} |
|
|
|
//! @} |
|
|
|
|
|
|
|
// |
|
|
|
void swap(Matrix& src) noexcept { |
|
|
|
std::swap(vector_storage_, src.vector_storage_); |
|
|
|
std::swap(raw_storage_, src.raw_storage_); |
|
|
|
std::swap(data_, src.data_); |
|
|
|
std::swap(use_vector_, src.use_vector_); |
|
|
|
std::swap(rows_, src.rows_); |
|
|
|
std::swap(cols_, src.cols_); |
|
|
|
} |
|
|
|
|
|
|
|
private: |
|
|
|
//! move helper |
|
|
|
void moves(Matrix&& src) noexcept { |
|
|
|
vector_storage_ = std::move(src.vector_storage_); |
|
|
|
raw_storage_ = std::move(src.raw_storage_); |
|
|
|
data_ = std::move(src.data_); |
|
|
|
use_vector_ = std::move(src.use_vector_); |
|
|
|
rows_ = std::move(src.rows_); |
|
|
|
cols_ = std::move(src.cols_); |
|
|
|
} |
|
|
|
|
|
|
|
// Storage |
|
|
|
std::vector<DataType> |
|
|
|
vector_storage_; //!< Internal storage (if used). |
|
|
|
DataType* raw_storage_; //!< External storage (if used). |
|
|
|
DataType* data_; //!< Pointer to active storage. |
|
|
|
bool use_vector_; //!< True if using vector storage, false for raw pointer. |
|
|
|
IndexType rows_{}; //!< the virtual size of rows. |
|
|
|
IndexType cols_{}; //!< the virtual size of columns. |
|
|
|
}; |
|
|
|
|
|
|
|
|
|
|
|
/** |
|
|
|
* A simple sparse matrix specialization. |
|
|
|
* |
|
|
|
* We use CSC format and provide get/set functionalities for each (i,j) item |
|
|
|
* on the matrix. We also provide a () overload using a proxy MatVal object. |
|
|
|
* This way the user can: |
|
|
|
* \code |
|
|
|
* auto v = A(3,4); |
|
|
|
* A(3, 4) = 7; |
|
|
|
* \endcode |
|
|
|
* |
|
|
|
* We also provide getCol() and getRow() functions witch return a viewer/iterator to rows and |
|
|
|
* columns of the matrix. In the case of a symmetric matrix instead of a row we return the |
|
|
|
* equivalent column. This way we gain speed due to CSC format nature. |
|
|
|
* |
|
|
|
* @tparam DataType The type for values |
|
|
|
* @tparam IndexType The type for indexes |
|
|
|
* @tparam Type The Matrix type (FULL or SYMMETRIC) |
|
|
|
*/ |
|
|
|
template<typename DataType, typename IndexType, |
|
|
|
MatrixOrder Order, |
|
|
|
bool Symmetric> |
|
|
|
struct Matrix<DataType, IndexType, MatrixType::SPARSE, Order, Symmetric> { |
|
|
|
|
|
|
|
using dataType = DataType; //!< meta:export of underling data type |
|
|
|
using indexType = IndexType; //!< meta:export of underling index type |
|
|
|
static constexpr MatrixOrder matrixOrder = Order; //!< meta:export of array order |
|
|
|
static constexpr MatrixType matrixType = MatrixType::SPARSE; //!< meta:export of array type |
|
|
|
static constexpr bool symmetric = Symmetric; //!< meta:export symmetric flag |
|
|
|
|
|
|
|
friend struct MatCol<Matrix>; |
|
|
|
friend struct MatRow<Matrix>; |
|
|
|
friend struct MatVal<Matrix>; |
|
|
|
|
|
|
|
/*! |
|
|
|
* \name Obj lifetime |
|
|
|
*/ |
|
|
|
//! @{ |
|
|
|
|
|
|
|
//! Default ctor with optional memory allocations |
|
|
|
Matrix(IndexType n=IndexType{}) noexcept: |
|
|
|
values{}, |
|
|
|
rows{}, |
|
|
|
col_ptr((n)? n+1:2, IndexType{}), |
|
|
|
N(n), |
|
|
|
NNZ(0) { } |
|
|
|
|
|
|
|
//! A ctor using csc array data |
|
|
|
Matrix(IndexType n, IndexType nnz, const IndexType* row, const IndexType* col) noexcept: |
|
|
|
values(nnz, 1), |
|
|
|
rows(row, row+nnz), |
|
|
|
col_ptr(col, col+n+1), |
|
|
|
N(n), |
|
|
|
NNZ(nnz) { } |
|
|
|
|
|
|
|
//! ctor using csc array data with value array |
|
|
|
Matrix(IndexType n, IndexType nnz, const DataType* v, const IndexType* row, const IndexType* col) noexcept: |
|
|
|
values(v, v+nnz), |
|
|
|
rows(row, row+nnz), |
|
|
|
col_ptr(col, col+n+1), |
|
|
|
N(n), |
|
|
|
NNZ(nnz) { } |
|
|
|
|
|
|
|
//! ctor vectors of row/col and default value for values array |
|
|
|
Matrix(IndexType n, IndexType nnz, const DataType v, |
|
|
|
const std::vector<IndexType>& row, const std::vector<IndexType>& col) noexcept: |
|
|
|
values(nnz, v), |
|
|
|
rows (row), |
|
|
|
col_ptr(col), |
|
|
|
N(n), |
|
|
|
NNZ(nnz) { } |
|
|
|
|
|
|
|
//! move ctor |
|
|
|
Matrix(Matrix&& m) noexcept { moves(std::move(m)); } |
|
|
|
//! move |
|
|
|
Matrix& operator=(Matrix&& m) noexcept { moves(std::move(m)); return *this; } |
|
|
|
Matrix(const Matrix& m) = delete; //!< make sure there are no copies |
|
|
|
Matrix& operator=(const Matrix& m) = delete; //!< make sure there are no copies |
|
|
|
//! @} |
|
|
|
|
|
|
|
//! \name Data exposure |
|
|
|
//! @{ |
|
|
|
|
|
|
|
//! \return the dimension of the matrix |
|
|
|
IndexType size() noexcept { return N; } |
|
|
|
//! After construction size configuration tool |
|
|
|
IndexType resize(IndexType n) { |
|
|
|
col_ptr.resize(n+1); |
|
|
|
return N = n; |
|
|
|
} |
|
|
|
//! \return the NNZ of the matrix |
|
|
|
IndexType capacity() noexcept { return NNZ; } |
|
|
|
//! After construction NNZ size configuration tool |
|
|
|
IndexType capacity(IndexType nnz) noexcept { |
|
|
|
values.reserve(nnz); |
|
|
|
rows.reserve(nnz); |
|
|
|
return NNZ; |
|
|
|
} |
|
|
|
// getters for row arrays of the struct (unused) |
|
|
|
std::vector<DataType>& getValues() noexcept { return values; } |
|
|
|
std::vector<IndexType>& getRows() noexcept { return rows; } |
|
|
|
std::vector<IndexType>& getCols() noexcept { return col_ptr; } |
|
|
|
|
|
|
|
/*! |
|
|
|
* Return a proxy MatVal object with read and write capabilities. |
|
|
|
* @param i The row number |
|
|
|
* @param j The column number |
|
|
|
* @return tHE MatVal object |
|
|
|
*/ |
|
|
|
MatVal<Matrix> operator()(IndexType i, IndexType j) noexcept { |
|
|
|
return MatVal<Matrix>(this, get(i, j), i, j); |
|
|
|
} |
|
|
|
|
|
|
|
/*! |
|
|
|
* A read item functionality using binary search to find the correct row |
|
|
|
* |
|
|
|
* @param i The row number |
|
|
|
* @param j The column number |
|
|
|
* @return The value of the item or DataType{} if is not present. |
|
|
|
*/ |
|
|
|
DataType get(IndexType i, IndexType j) noexcept { |
|
|
|
IndexType idx; bool found; |
|
|
|
std::tie(idx, found) =find_idx(rows, col_ptr[j], col_ptr[j+1], i); |
|
|
|
return (found) ? values[idx] : 0; |
|
|
|
} |
|
|
|
|
|
|
|
/*! |
|
|
|
* A write item functionality. |
|
|
|
* |
|
|
|
* First we search if the matrix has already a value in (i, j) position. |
|
|
|
* If so we just change it to a new value. If not we add the item on the matrix. |
|
|
|
* |
|
|
|
* @note |
|
|
|
* When change a value, we don't increase the NNZ value of the struct. We expect the user has already |
|
|
|
* change the NNZ value to the right one using @see capacity() function. When adding a value we |
|
|
|
* increase the NNZ. |
|
|
|
* |
|
|
|
* @param i The row number |
|
|
|
* @param j The column number |
|
|
|
* @return The new value of the item . |
|
|
|
*/ |
|
|
|
DataType set(DataType v, IndexType i, IndexType j) { |
|
|
|
IndexType idx; bool found; |
|
|
|
std::tie(idx, found) = find_idx(rows, col_ptr[j], col_ptr[j+1], i); |
|
|
|
if (found) |
|
|
|
return values[idx] = v; // we don't change NNZ even if we write "0" |
|
|
|
else { |
|
|
|
values.insert(values.begin()+idx, v); |
|
|
|
rows.insert(rows.begin()+idx, i); |
|
|
|
std::transform(col_ptr.begin()+j+1, col_ptr.end(), col_ptr.begin()+j+1, [](IndexType it) { |
|
|
|
return ++it; |
|
|
|
}); |
|
|
|
++NNZ; // we increase the NNZ even if we write "0" |
|
|
|
return v; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
/*! |
|
|
|
* Get a view of a CSC column |
|
|
|
* @param j The column to get |
|
|
|
* @return The MatCol object @see MatCol |
|
|
|
*/ |
|
|
|
MatCol<Matrix> getCol(IndexType j) noexcept { |
|
|
|
return MatCol<Matrix>(this, col_ptr[j], col_ptr[j+1]); |
|
|
|
} |
|
|
|
|
|
|
|
/*! |
|
|
|
* Get a view of a CSC row |
|
|
|
* |
|
|
|
* In case of a SYMMETRIC matrix we can return a column instead. |
|
|
|
* |
|
|
|
* @param j The row to get |
|
|
|
* @return On symmetric matrix MatCol otherwise a MatRow |
|
|
|
*/ |
|
|
|
|
|
|
|
MatCol<Matrix> getRow(IndexType i) noexcept { |
|
|
|
if constexpr (Symmetric) |
|
|
|
return getCol(i); |
|
|
|
else |
|
|
|
return MatRow<Matrix>(this, i); |
|
|
|
} |
|
|
|
|
|
|
|
// values only iterator support |
|
|
|
DataType* begin() noexcept { return values.begin(); } |
|
|
|
DataType* end() noexcept { return values.end(); } |
|
|
|
//! @} |
|
|
|
|
|
|
|
//! A small iteration helper |
|
|
|
template<typename F, typename... Args> |
|
|
|
void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) { |
|
|
|
for (IndexType it=begin ; it<end ; ++it) { |
|
|
|
std::forward<F>(lambda)(std::forward<Args>(args)..., it); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
private: |
|
|
|
/*! |
|
|
|
* A small binary search implementation using index for begin-end instead of iterators. |
|
|
|
* |
|
|
|
* \param v Reference to vector to search |
|
|
|
* \param begin The vector's index to begin |
|
|
|
* \param end The vector's index to end |
|
|
|
* \param match What to search |
|
|
|
* \return An <index, status> pair. |
|
|
|
* index is the index of the item or end if not found |
|
|
|
* status is true if found, false otherwise |
|
|
|
*/ |
|
|
|
std::pair<IndexType, bool> find_idx(const std::vector<IndexType>& v, IndexType begin, IndexType end, IndexType match) { |
|
|
|
if (v.capacity() != 0 && begin < end) { |
|
|
|
IndexType b = begin, e = end-1; |
|
|
|
while (b <= e) { |
|
|
|
IndexType m = (b+e)/2; |
|
|
|
if (v[m] == match) return std::make_pair(m, true); |
|
|
|
else if (b >= e) return std::make_pair(end, false); |
|
|
|
else { |
|
|
|
if (v[m] < match) b = m +1; |
|
|
|
else e = m -1; |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
return std::make_pair(end, false); |
|
|
|
} |
|
|
|
|
|
|
|
// move helper |
|
|
|
void moves(Matrix&& src) noexcept { |
|
|
|
values = std::move(src.values); |
|
|
|
rows = std::move(src.rows); |
|
|
|
col_ptr = std::move(src.col_ptr); |
|
|
|
N = std::move(src.N); // redundant for primitives |
|
|
|
NNZ = std::move(src.NNZ); // |
|
|
|
} |
|
|
|
//! \name Data |
|
|
|
//! @{ |
|
|
|
std::vector<DataType> values {}; //!< vector to store the values of the matrix |
|
|
|
std::vector<IndexType> rows{}; //!< vector to store the row information |
|
|
|
std::vector<IndexType> col_ptr{1,0}; //!< vector to store the column pointers |
|
|
|
IndexType N{0}; //!< The dimension of the matrix (square) |
|
|
|
IndexType NNZ{0}; //!< The NNZ (capacity of the matrix) |
|
|
|
//! @} |
|
|
|
}; |
|
|
|
|
|
|
|
|
|
|
|
/*! |
|
|
|
* A view/iterator hybrid object for Matrix columns. |
|
|
|
* |
|
|
|
* This object provides access to a column of a Matrix. The public functionalities |
|
|
|
* allow data access using indexes instead of iterators. We prefer indexes over iterators |
|
|
|
* because we can apply the same index to different inner vector of Matrix without conversion. |
|
|
|
* |
|
|
|
* @tparam DataType |
|
|
|
* @tparam IndexType |
|
|
|
*/ |
|
|
|
template<typename MatrixType> |
|
|
|
struct MatCol { |
|
|
|
using owner_t = MatrixType; |
|
|
|
|
|
|
|
using DataType = typename MatrixType::dataType; |
|
|
|
using IndexType = typename MatrixType::indexType; |
|
|
|
|
|
|
|
/*! |
|
|
|
* ctor using column pointers for begin-end. own is pointer to Matrix. |
|
|
|
*/ |
|
|
|
MatCol(owner_t* own, const IndexType begin, const IndexType end) noexcept : |
|
|
|
owner_(own), index_(begin), begin_(begin), end_(end) { |
|
|
|
vindex_ = vIndexCalc(index_); |
|
|
|
} |
|
|
|
MatCol() = default; |
|
|
|
MatCol(const MatCol&) = delete; //!< make sure there are no copies |
|
|
|
MatCol& operator=(const MatCol&)= delete; //!< make sure there are no copies |
|
|
|
MatCol(MatCol&&) = default; |
|
|
|
MatCol& operator=(MatCol&&) = default; |
|
|
|
|
|
|
|
//! a simple dereference operator, like an iterator |
|
|
|
DataType operator* () { |
|
|
|
return get(); |
|
|
|
} |
|
|
|
//! Increment operator acts on index(), like an iterator |
|
|
|
MatCol& operator++ () { advance(); return *this; } |
|
|
|
MatCol& operator++ (int) { MatCol& p = *this; advance(); return p; } |
|
|
|
|
|
|
|
//! () operator acts as member access (like a view) |
|
|
|
DataType operator()(IndexType x) { |
|
|
|
return (x == index())? get() : DataType{}; |
|
|
|
} |
|
|
|
//! = operator acts as member assignment (like a view) |
|
|
|
DataType operator= (DataType v) { return owner_->values[index_] = v; } |
|
|
|
// iterator like handlers |
|
|
|
// these return a virtual index value based on the items position on the full matrix |
|
|
|
// but the move of the index is just a ++ away. |
|
|
|
IndexType index() noexcept { return vindex_; } |
|
|
|
const IndexType index() const noexcept { return vindex_; } |
|
|
|
IndexType begin() noexcept { return vIndexCalc(begin_); } |
|
|
|
const IndexType begin() const noexcept { return vIndexCalc(begin_); } |
|
|
|
IndexType end() noexcept { return owner_->N; } |
|
|
|
const IndexType end() const noexcept { return owner_->N; } |
|
|
|
|
|
|
|
/*! |
|
|
|
* Multiplication operator |
|
|
|
* |
|
|
|
* We follow only the non-zero values and multiply only the common indexes. |
|
|
|
* |
|
|
|
* @tparam C Universal reference for the type right half site column |
|
|
|
* |
|
|
|
* @param c The right hand site matrix |
|
|
|
* @return The value of the inner product of two vectors |
|
|
|
* @note The time complexity is \$ O(nnz1+nnz2) \$. |
|
|
|
* Where the nnz is the max NNZ elements of the column of the matrix |
|
|
|
*/ |
|
|
|
template <typename C> |
|
|
|
DataType operator* (C&& c) { |
|
|
|
static_assert(std::is_same<remove_cvref_t<C>, MatCol<MatrixType>>(), ""); |
|
|
|
DataType v{}; |
|
|
|
while (index() != end() && c.index() != c.end()) { |
|
|
|
if (index() < c.index()) advance(); // advance me |
|
|
|
else if (index() > c.index()) ++c; // advance other |
|
|
|
else { //index() == c.index() |
|
|
|
v += get() * *c; // multiply and advance both |
|
|
|
++c; |
|
|
|
advance(); |
|
|
|
} |
|
|
|
} |
|
|
|
return v; |
|
|
|
} |
|
|
|
|
|
|
|
private: |
|
|
|
//! small tool to increase the index pointers to Matrix |
|
|
|
void advance() noexcept { |
|
|
|
++index_; |
|
|
|
vindex_ = vIndexCalc(index_); |
|
|
|
} |
|
|
|
//! tool to translate between col_ptr indexes and Matrix "virtual" full matrix indexes |
|
|
|
IndexType vIndexCalc(IndexType idx) { |
|
|
|
return (idx < end_) ? owner_->rows[idx] : end(); |
|
|
|
} |
|
|
|
//! small get tool |
|
|
|
DataType get() { return owner_->values[index_]; } |
|
|
|
|
|
|
|
owner_t* owner_ {nullptr}; //!< Pointer to owner Matrix. MatCol is just a view |
|
|
|
IndexType vindex_ {IndexType{}}; //!< Virtual index of full matrix |
|
|
|
IndexType index_ {IndexType{}}; //!< index to Matrix::rows |
|
|
|
IndexType begin_ {IndexType{}}; //!< beginning index of the column in Matrix::rows |
|
|
|
IndexType end_ {IndexType{}}; //!< ending index of the column in Matrix::rows |
|
|
|
}; |
|
|
|
|
|
|
|
/*! |
|
|
|
* A view/iterator hybrid object for Matrix rows. |
|
|
|
* |
|
|
|
* This object provides access to a column of a Matrix. The public functionalities |
|
|
|
* allow data access using indexes instead of iterators. We prefer indexes over iterators |
|
|
|
* because we can apply the same index to different inner vector of Matrix without conversion. |
|
|
|
* |
|
|
|
* @tparam DataType |
|
|
|
* @tparam IndexType |
|
|
|
*/ |
|
|
|
template<typename MatrixType> |
|
|
|
struct MatRow { |
|
|
|
using owner_t = MatrixType; |
|
|
|
|
|
|
|
using DataType = typename MatrixType::dataType; |
|
|
|
using IndexType = typename MatrixType::indexType; |
|
|
|
|
|
|
|
/*! |
|
|
|
* ctor using virtual full matrix row index. own is pointer to Matrix. |
|
|
|
*/ |
|
|
|
MatRow(owner_t* own, const IndexType row) noexcept : |
|
|
|
owner_(own), vindex_(IndexType{}), row_(row), index_(IndexType{}), |
|
|
|
begin_(IndexType{}), end_(owner_->NNZ) { |
|
|
|
// place begin |
|
|
|
while(begin_ != end_ && owner_->rows[begin_] != row_) |
|
|
|
++begin_; |
|
|
|
// place index_ and vindex_ |
|
|
|
if (owner_->rows[index_] != row_) |
|
|
|
advance(); |
|
|
|
} |
|
|
|
MatRow() = default; |
|
|
|
MatRow(const MatRow&) = delete; //!< make sure there are no copies |
|
|
|
MatRow& operator=(const MatRow&)= delete; //!< make sure there are no copies |
|
|
|
MatRow(MatRow&&) = default; |
|
|
|
MatRow& operator=(MatRow&&) = default; |
|
|
|
|
|
|
|
//! a simple dereference operator, like an iterator |
|
|
|
DataType operator* () { |
|
|
|
return get(); |
|
|
|
} |
|
|
|
//! Increment operator acts on index(), like an iterator |
|
|
|
//! here the increment is a O(N) process. |
|
|
|
MatRow& operator++ () { advance(); return *this; } |
|
|
|
MatRow& operator++ (int) { MatRow& p = *this; advance(); return p; } |
|
|
|
|
|
|
|
//! () operator acts as member access (like a view) |
|
|
|
DataType operator()(IndexType x) { |
|
|
|
return (x == index())? get() : DataType{}; |
|
|
|
} |
|
|
|
//! = operator acts as member assignment (like a view) |
|
|
|
DataType operator= (DataType v) { return owner_->values[index_] = v; } |
|
|
|
// iterator like handlers |
|
|
|
// these return a virtual index value based on the items position on the full matrix |
|
|
|
// but the move of the index is just a ++ away. |
|
|
|
IndexType index() noexcept { return vindex_; } |
|
|
|
const IndexType index() const noexcept { return vindex_; } |
|
|
|
IndexType begin() noexcept { return vIndexCalc(begin_); } |
|
|
|
const IndexType begin() const noexcept { return vIndexCalc(begin_); } |
|
|
|
IndexType end() noexcept { return owner_->N; } |
|
|
|
const IndexType end() const noexcept { return owner_->N; } |
|
|
|
|
|
|
|
/*! |
|
|
|
* Multiplication operator |
|
|
|
* |
|
|
|
* We follow only the non-zero values and multiply only the common indexes. |
|
|
|
* |
|
|
|
* @tparam C Universal reference for the type right half site column |
|
|
|
* |
|
|
|
* @param c The right hand site matrix |
|
|
|
* @return The value of the inner product of two vectors |
|
|
|
* @note The time complexity is \$ O(N+nnz2) \$ and way heavier the ColxCol multiplication. |
|
|
|
* Where the nnz is the max NNZ elements of the column of the matrix |
|
|
|
*/ |
|
|
|
template <typename C> |
|
|
|
DataType operator* (C&& c) { |
|
|
|
static_assert(std::is_same<remove_cvref_t<C>, MatCol<MatrixType>>(), ""); |
|
|
|
DataType v{}; |
|
|
|
while (index() != end() && c.index() != c.end()) { |
|
|
|
if (index() < c.index()) advance(); // advance me |
|
|
|
else if (index() > c.index()) ++c; // advance other |
|
|
|
else { //index() == c.index() |
|
|
|
v += get() * *c; // multiply and advance both |
|
|
|
++c; |
|
|
|
advance(); |
|
|
|
} |
|
|
|
} |
|
|
|
return v; |
|
|
|
} |
|
|
|
private: |
|
|
|
//! small tool to increase the index pointers to Matrix matrix |
|
|
|
//! We have to search the entire rows vector in Matrix to find the next |
|
|
|
//! virtual row position. |
|
|
|
//! time complexity O(N) |
|
|
|
void advance() noexcept { |
|
|
|
do |
|
|
|
++index_; |
|
|
|
while(index_ != end_ && owner_->rows[index_] != row_); |
|
|
|
vindex_ = vIndexCalc(index_); |
|
|
|
} |
|
|
|
//! tool to translate between col_ptr indexes and Matrix "virtual" full matrix indexes |
|
|
|
IndexType vIndexCalc(IndexType idx) { |
|
|
|
for(IndexType i =0 ; i<(owner_->N+1) ; ++i) |
|
|
|
if (idx < owner_->col_ptr[i]) |
|
|
|
return i-1; |
|
|
|
return end(); |
|
|
|
} |
|
|
|
//! small get tool |
|
|
|
DataType get() { return owner_->values[index_]; } |
|
|
|
|
|
|
|
owner_t* owner_ {nullptr}; //!< Pointer to owner Matrix. MatCol is just a view |
|
|
|
IndexType vindex_ {IndexType{}}; //!< Virtual index of full matrix |
|
|
|
IndexType row_ {IndexType{}}; //!< The virtual full matrix row of the object |
|
|
|
IndexType index_ {IndexType{}}; //!< index to Matrix::rows |
|
|
|
IndexType begin_ {IndexType{}}; //!< beginning index of the column in Matrix::rows |
|
|
|
IndexType end_ {IndexType{}}; //!< ending index of the column in Matrix::rows |
|
|
|
}; |
|
|
|
|
|
|
|
/*! |
|
|
|
* A proxy Matrix value object/view. |
|
|
|
* |
|
|
|
* This object acts as proxy to provide read/write access to an Matrix item. |
|
|
|
* |
|
|
|
* @tparam DataType The type of the values of the Matrix matrix |
|
|
|
* @tparam IndexType The type of the indexes of the Matrix matrix |
|
|
|
*/ |
|
|
|
template<typename MatrixType> |
|
|
|
struct MatVal { |
|
|
|
using owner_t = MatrixType; |
|
|
|
|
|
|
|
using DataType = typename MatrixType::dataType; |
|
|
|
using IndexType = typename MatrixType::indexType; |
|
|
|
|
|
|
|
//!< ctor using all value-row-column data, plus a pointer to owner Matrix object |
|
|
|
MatVal(owner_t* own, DataType v, IndexType i, IndexType j) : |
|
|
|
owner_(own), v_(v), i_(i), j_(j) { } |
|
|
|
MatVal() = default; |
|
|
|
MatVal(const MatVal&) = delete; //!< make sure there are no copies |
|
|
|
MatVal& operator=(const MatVal&) = delete; //!< make sure there are no copies |
|
|
|
MatVal(MatVal&&) = default; |
|
|
|
MatVal& operator=(MatVal&&) = default; |
|
|
|
|
|
|
|
//! Operator to return the DataType value implicitly |
|
|
|
operator DataType() { return v_; } |
|
|
|
//! Operator to write back to owner the assigned value |
|
|
|
//! for ex: A(2,3) = 5; |
|
|
|
MatVal& operator=(DataType v) { |
|
|
|
v_ = v; |
|
|
|
owner_->set(v_, i_, j_); |
|
|
|
return *this; |
|
|
|
} |
|
|
|
private: |
|
|
|
owner_t* owner_{nullptr}; //!< Pointer to owner Matrix. MatVal is just a view. |
|
|
|
DataType v_{DataType{}}; //!< The value of the row-column pair (for speed) |
|
|
|
IndexType i_{IndexType{}}; //!< The row |
|
|
|
IndexType j_{IndexType{}}; //!< the column |
|
|
|
}; |
|
|
|
|
|
|
|
|
|
|
|
} // namespace mtx |
|
|
|
|
|
|
|
|
|
|
|
#endif /* MATRIX_HPP_ */ |