723 lines
26 KiB
C++
723 lines
26 KiB
C++
/**
|
|
* \file impl.hpp
|
|
* \brief Implementation common header file
|
|
*
|
|
* \author
|
|
* Christos Choutouridis AEM:8997
|
|
* <cchoutou@ece.auth.gr>
|
|
*/
|
|
#ifndef IMPL_HPP_
|
|
#define IMPL_HPP_
|
|
|
|
#include <type_traits>
|
|
#include <utility>
|
|
#include <algorithm>
|
|
#include <vector>
|
|
#include <tuple>
|
|
#include <cstddef>
|
|
#include <iostream>
|
|
#include <fstream>
|
|
|
|
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 {
|
|
FULL, /*!< Matrix is asymmetric */
|
|
SYMMETRIC, /*!< Matrix is symmetric */
|
|
};
|
|
|
|
/*
|
|
* Forward type declerations
|
|
*/
|
|
template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC> struct Matrix;
|
|
template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC> struct SpMat;
|
|
template<typename DataType, typename IndexType> struct SpMatCol;
|
|
template<typename DataType, typename IndexType> struct SpMatRow;
|
|
template<typename DataType, typename IndexType> struct SpMatVal;
|
|
|
|
/*!
|
|
* 2D-array wrapper for v1 and v2 part of the exercise used as RAII
|
|
* and copy-prevention.
|
|
*
|
|
* 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 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, MatrixType Type>
|
|
struct Matrix {
|
|
|
|
using dataType = DataType; //!< meta:export of underling data type
|
|
using indexType = IndexType; //!< meta:export of underling index type
|
|
static constexpr MatrixType matrixType = Type; //!< export of array type
|
|
|
|
/*!
|
|
* \name Obj lifetime
|
|
*/
|
|
//! @{
|
|
|
|
//! Constructor using data type and size
|
|
Matrix(DataType* t, IndexType s) noexcept : m_(t), size_(s) { }
|
|
|
|
//! RAII deleter
|
|
~Matrix() { delete m_; }
|
|
//! move ctor
|
|
Matrix(Matrix&& a) noexcept { a.swap(*this); }
|
|
//! move
|
|
Matrix& operator=(Matrix&& a) noexcept { a.swap(*this); return *this; }
|
|
Matrix(const Matrix& a) = delete; //!< No copy ctor
|
|
Matrix& operator=(const Matrix& a) = delete; //!< No copy
|
|
//! @}
|
|
|
|
//! \name Data exposure
|
|
//! @{
|
|
|
|
//! memory capacity of the matrix with diagonal data
|
|
template<MatrixType AT=Type> std::enable_if_t<AT==MatrixType::SYMMETRIC, IndexType>
|
|
static constexpr capacity(IndexType N) noexcept { return (N+1)*N/2; }
|
|
//! memory capacity for full matrix
|
|
template<MatrixType AT=Type> std::enable_if_t<AT==MatrixType::FULL, IndexType>
|
|
static constexpr capacity(IndexType N) noexcept { return N*N; }
|
|
|
|
//! Get the size of each dimension
|
|
IndexType size() noexcept { return size_; }
|
|
|
|
/*
|
|
* virtual 2D accessors
|
|
*/
|
|
template<MatrixType AT=Type>
|
|
std::enable_if_t <AT==MatrixType::SYMMETRIC, DataType>
|
|
get (IndexType i, IndexType j) {
|
|
if (i<j) return m_[j*(j+1)/2 + i]; // Upper, use opposite index
|
|
else return m_[i*(i+1)/2 + j]; // Lower, use our notation
|
|
}
|
|
template<MatrixType AT=Type>
|
|
std::enable_if_t <AT==MatrixType::FULL, DataType>
|
|
get(IndexType i, IndexType j) {
|
|
return m_[i*size_ + j];
|
|
}
|
|
template<MatrixType AT=Type>
|
|
std::enable_if_t <AT==MatrixType::SYMMETRIC, DataType>
|
|
set(DataType v, IndexType i, IndexType j) {
|
|
if (i<j) return m_[j*(j+1)/2 + i] =v; // Upper, use opposite index
|
|
else return m_[i*(i+1)/2 + j] =v; // Lower, use our notation
|
|
}
|
|
template<MatrixType AT=Type>
|
|
std::enable_if_t <AT==MatrixType::FULL, DataType>
|
|
set(DataType v, IndexType i, IndexType j) {
|
|
return m_[i*size_ + j] =v;
|
|
}
|
|
DataType operator()(IndexType i, IndexType j) { return get(i, j); }
|
|
|
|
// a basic serial iterator support
|
|
DataType* begin() noexcept { return m_; }
|
|
DataType* end() noexcept { return m_ + Matrix::capacity(size_); }
|
|
//! @}
|
|
|
|
/*!
|
|
* \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);
|
|
}
|
|
}
|
|
//! @}
|
|
|
|
// move helper
|
|
void swap(Matrix& src) noexcept {
|
|
std::swap(m_, src.m_);
|
|
std::swap(size_, src.size_);
|
|
}
|
|
private:
|
|
DataType* m_; //!< Pointer to actual data.
|
|
IndexType size_; //!< the virtual size of each dimension.
|
|
};
|
|
|
|
/*!
|
|
* RAII allocation helper, the smart_ptr way.
|
|
*/
|
|
template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC>
|
|
Matrix<DataType, IndexType, Type> make_Matrix(IndexType s) {
|
|
return Matrix<DataType, IndexType, Type>(new DataType[Matrix<DataType, IndexType, Type>::capacity(s)], s);
|
|
}
|
|
|
|
|
|
/**
|
|
* A simple sparse matrix implementation.
|
|
*
|
|
* 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 SpMatVal 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, MatrixType Type>
|
|
struct SpMat {
|
|
using dataType = DataType; //!< meta:export of underling data type
|
|
using indexType = IndexType; //!< meta:export of underling index type
|
|
static constexpr MatrixType matrixType = Type; //!< export of array type
|
|
|
|
friend struct SpMatCol<DataType, IndexType>;
|
|
friend struct SpMatRow<DataType, IndexType>;
|
|
friend struct SpMatVal<DataType, IndexType>;
|
|
|
|
/*!
|
|
* \name Obj lifetime
|
|
*/
|
|
//! @{
|
|
|
|
//! Default ctor with optional memory allocations
|
|
SpMat(IndexType n=IndexType{}, IndexType nnz=IndexType{}) :
|
|
values(nnz, DataType{}),
|
|
rows(nnz, IndexType{}),
|
|
col_ptr((n)? n+1:2, IndexType{}),
|
|
N(n),
|
|
NNZ(nnz) { }
|
|
|
|
//! A ctor using csc array data
|
|
SpMat(IndexType n, IndexType nnz, const IndexType* row, const IndexType* col) :
|
|
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
|
|
SpMat(IndexType n, IndexType nnz, const DataType* v, const IndexType* row, const IndexType* col) :
|
|
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
|
|
SpMat(IndexType n, IndexType nnz, const DataType v, const std::vector<IndexType>& row, const std::vector<IndexType>& col) :
|
|
values(nnz, v),
|
|
rows (row),
|
|
col_ptr(col),
|
|
N(n),
|
|
NNZ(nnz) { }
|
|
|
|
//! move ctor
|
|
SpMat(SpMat&& a) noexcept { moves(std::move(a)); }
|
|
//! move
|
|
SpMat& operator=(SpMat&& a) noexcept { moves(std::move(a)); return *this; }
|
|
SpMat(const SpMat& a) = delete; //!< make sure there are no copies
|
|
SpMat& operator=(const SpMat& a) = 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 size(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) {
|
|
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 SpMatVal object with read and write capabilities.
|
|
* @param i The row number
|
|
* @param j The column number
|
|
* @return tHE SpMatVal object
|
|
*/
|
|
SpMatVal<DataType, IndexType> operator()(IndexType i, IndexType j) {
|
|
return SpMatVal<DataType, IndexType>(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) {
|
|
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 read item functionality using linear 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_lin(IndexType i, IndexType j) {
|
|
IndexType idx; bool found;
|
|
std::tie(idx, found) =find_lin_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
|
|
* 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.
|
|
*
|
|
* @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_lin_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 SpMatCol object @see SpMatCol
|
|
*/
|
|
SpMatCol<DataType, IndexType> getCol(IndexType j) {
|
|
return SpMatCol<DataType, IndexType>(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 The SpMatCol object @see SpMatCol
|
|
*/
|
|
template<MatrixType AT= Type>
|
|
std::enable_if_t<AT==MatrixType::SYMMETRIC, SpMatCol<DataType, IndexType>>
|
|
getRow(IndexType i) {
|
|
return getCol(i);
|
|
}
|
|
|
|
/*!
|
|
* Get a view of a CSC row
|
|
*
|
|
* @param j The row to get
|
|
* @return The SpMatRow object @see SpMatRow
|
|
*/
|
|
template<MatrixType AT= Type>
|
|
std::enable_if_t<AT==MatrixType::FULL, SpMatCol<DataType, IndexType>>
|
|
getRow(IndexType i) {
|
|
return SpMatRow<DataType, IndexType>(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);
|
|
}
|
|
}
|
|
|
|
// friend operations for printing
|
|
template<typename D, typename I, MatrixType T> friend void print(SpMat<D, I, T>& mat);
|
|
template<typename D, typename I, MatrixType T> friend void print_dense(SpMat<D, I, T>& mat);
|
|
|
|
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 The index of the item or end on failure.
|
|
*/
|
|
std::pair<IndexType, bool> find_idx(const std::vector<IndexType>& v, IndexType begin, IndexType end, IndexType match) {
|
|
IndexType b = begin, e = end-1;
|
|
while (true) {
|
|
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);
|
|
}
|
|
/*!
|
|
* find helper for set using index for begin-end instead of iterators.
|
|
*
|
|
* We search for the item or a place to add the item.
|
|
* So we return the index if we find it. Otherwise we return the place to add it.
|
|
*
|
|
* \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
|
|
*/
|
|
std::pair<IndexType, bool> find_lin_idx(const std::vector<IndexType>& v, IndexType begin, IndexType end, IndexType match) {
|
|
for ( ; begin < end ; ++begin) {
|
|
if (match == v[begin]) return std::make_pair(begin, true);
|
|
else if (match < v[begin]) return std::make_pair(begin, false);
|
|
}
|
|
return std::make_pair(end, false);
|
|
}
|
|
// move helper
|
|
void moves(SpMat&& 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 stor 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 SpMat columns.
|
|
*
|
|
* This object provides access to a column of a SpMat. 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 SpMat without conversion.
|
|
*
|
|
* @tparam DataType
|
|
* @tparam IndexType
|
|
*/
|
|
template<typename DataType, typename IndexType>
|
|
struct SpMatCol {
|
|
using owner_t = SpMat<DataType, IndexType>;
|
|
|
|
/*!
|
|
* ctor using column pointers for begin-end. own is pointer to SpMat.
|
|
*/
|
|
SpMatCol(owner_t* own, const IndexType begin, const IndexType end) noexcept :
|
|
owner_(own), index_(begin), begin_(begin), end_(end) {
|
|
vindex_ = vIndexCalc(index_);
|
|
}
|
|
SpMatCol() = default;
|
|
SpMatCol(const SpMatCol&) = delete; //!< make sure there are no copies
|
|
SpMatCol& operator=(const SpMatCol&)= delete; //!< make sure there are no copies
|
|
SpMatCol(SpMatCol&&) = default;
|
|
SpMatCol& operator=(SpMatCol&&) = default;
|
|
|
|
//! a simple dereference operator, like an iterator
|
|
DataType operator* () {
|
|
return get();
|
|
}
|
|
//! Increment operator acts on index(), like an iterator
|
|
SpMatCol& operator++ () { advance(); return *this; }
|
|
SpMatCol& operator++ (int) { SpMatCol& 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
|
|
* @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
|
|
*/
|
|
template <typename C>
|
|
DataType operator* (C&& c) {
|
|
static_assert(std::is_same<remove_cvref_t<C>, SpMatCol<DataType, IndexType>>(), "");
|
|
DataType v{};
|
|
while (index() != end() && c.index() != c.end()) {
|
|
if (index() < c.index()) advance();
|
|
else if (index() > c.index()) ++c;
|
|
else { //index() == c.index()
|
|
v += get() * *c;
|
|
++c;
|
|
advance();
|
|
}
|
|
}
|
|
return v;
|
|
}
|
|
|
|
private:
|
|
//! small tool to increase the index pointers to SpMat matrix
|
|
void advance() noexcept {
|
|
++index_;
|
|
vindex_ = vIndexCalc(index_);
|
|
}
|
|
//! tool to translate between col_ptr indexes and SpMat "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 SpMat. SpMatCol is just a view
|
|
IndexType vindex_ {IndexType{}}; //!< Virtual index of full matrix
|
|
IndexType index_ {IndexType{}}; //!< index to SpMat::rows
|
|
IndexType begin_ {IndexType{}}; //!< beginning index of the column in SpMat::rows
|
|
IndexType end_ {IndexType{}}; //!< ending index of the column in SpMat::rows
|
|
};
|
|
|
|
/*!
|
|
* A view/iterator hybrid object for SpMat rows.
|
|
*
|
|
* This object provides access to a column of a SpMat. 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 SpMat without conversion.
|
|
*
|
|
* @tparam DataType
|
|
* @tparam IndexType
|
|
*/
|
|
template<typename DataType, typename IndexType>
|
|
struct SpMatRow {
|
|
using owner_t = SpMat<DataType, IndexType>;
|
|
|
|
/*!
|
|
* ctor using virtual full matrix row index. own is pointer to SpMat.
|
|
*/
|
|
SpMatRow(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();
|
|
}
|
|
SpMatRow() = default;
|
|
SpMatRow(const SpMatRow&) = delete; //!< make sure there are no copies
|
|
SpMatRow& operator=(const SpMatRow&)= delete; //!< make sure there are no copies
|
|
SpMatRow(SpMatRow&&) = default;
|
|
SpMatRow& operator=(SpMatRow&&) = 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.
|
|
SpMatRow& operator++ () { advance(); return *this; }
|
|
SpMatRow& operator++ (int) { SpMatRow& 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
|
|
* @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
|
|
*/
|
|
template <typename C>
|
|
DataType operator* (C&& c) {
|
|
static_assert(std::is_same<remove_cvref_t<C>, SpMatCol<DataType, IndexType>>(), "");
|
|
DataType v{};
|
|
while (index() != end() && c.index() != c.end()) {
|
|
if (index() < c.index()) advance();
|
|
else if (index() > c.index()) ++c;
|
|
else { //index() == c.index()
|
|
v += get()* *c;
|
|
++c;
|
|
advance();
|
|
}
|
|
}
|
|
return v;
|
|
}
|
|
private:
|
|
//! small tool to increase the index pointers to SpMat matrix
|
|
//! We have to search the entire rows vector in SpMat 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 SpMat "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 SpMat. SpMatCol 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 SpMat::rows
|
|
IndexType begin_ {IndexType{}}; //!< beginning index of the column in SpMat::rows
|
|
IndexType end_ {IndexType{}}; //!< ending index of the column in SpMat::rows
|
|
};
|
|
|
|
/*!
|
|
* A proxy SpMat value object/view.
|
|
*
|
|
* This object acts as proxy to provide read/write access to an SpMat item.
|
|
*
|
|
* @tparam DataType The type of the values of the SpMat matrix
|
|
* @tparam IndexType The type of the indexes of the SpMat matrix
|
|
*/
|
|
template<typename DataType, typename IndexType>
|
|
struct SpMatVal {
|
|
using owner_t = SpMat<DataType, IndexType>;
|
|
|
|
//!< ctor using all value-row-column data, plus a pointer to owner SpMat object
|
|
SpMatVal(owner_t* own, DataType v, IndexType i, IndexType j) :
|
|
owner_(own), v_(v), i_(i), j_(j) { }
|
|
SpMatVal() = default;
|
|
SpMatVal(const SpMatVal&) = delete; //!< make sure there are no copies
|
|
SpMatVal& operator=(const SpMatVal&) = delete; //!< make sure there are no copies
|
|
SpMatVal(SpMatVal&&) = default;
|
|
SpMatVal& operator=(SpMatVal&&) = 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;
|
|
SpMatVal& operator=(DataType v) {
|
|
v_ = v;
|
|
owner_->set(v_, i_, j_);
|
|
return *this;
|
|
}
|
|
private:
|
|
owner_t* owner_{nullptr}; //!< Pointer to owner SpMat. SpMatVal 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
|
|
};
|
|
|
|
|
|
//! enumerator for input matrix type.
|
|
enum class InputMatrix{ UNSPECIFIED, GENERATE, MTX };
|
|
//! enumerator for output handling
|
|
enum class OutputMode{ STD, FILE };
|
|
|
|
/*!
|
|
* Session option for each invocation of the executable
|
|
*/
|
|
struct session_t {
|
|
InputMatrix inputMatrix {InputMatrix::UNSPECIFIED}; //!< Source of the matrix
|
|
std::ifstream mtxFile {}; //!< matrix file in MatrixMarket format
|
|
std::size_t gen_size {}; //!< size of the matrix if we select to generate a random one
|
|
double gen_prob {}; //!< probability of the binomial distribution for the matrix
|
|
//!< if we generate one
|
|
OutputMode outputMode {OutputMode::STD}; //!< Type of the output file
|
|
std::ofstream outFile {}; //!< File to use for output
|
|
std::size_t max_threads {}; //!< Maximum threads to use
|
|
std::size_t repeat {1}; //!< How many times we execute the calculations part
|
|
bool timing {false}; //!< Enable timing prints of the program
|
|
bool verbose {false}; //!< Flag to enable verbose output to stdout
|
|
#if CODE_VERSION == 3
|
|
bool makeSymmetric {false}; //!< symmetric matrix creation flag (true by default)
|
|
#else
|
|
bool makeSymmetric {true}; //!< symmetric matrix creation flag (true by default)
|
|
#endif
|
|
bool validate_mtx {false}; //!< Flag to request mtx input data triangular validation.
|
|
bool print_count {false}; //!< Flag to request total count printing
|
|
bool mtx_print {false}; //!< matrix print flag
|
|
std::size_t mtx_print_size {}; //!< matrix print size
|
|
bool dynamic {false}; //!< Selects dynamic scheduling for OpenMP and pthreads.
|
|
};
|
|
|
|
extern session_t session;
|
|
|
|
|
|
#endif /* IMPL_HPP_ */
|