diff --git a/homework_2/Makefile b/homework_2/Makefile index c573b2e..d890a86 100644 --- a/homework_2/Makefile +++ b/homework_2/Makefile @@ -25,8 +25,9 @@ TARGET := distbitonic # Source directories list(space seperated). Makefile-relative path, UNDER current directory. SRC_DIR_LIST := src # Include directories list(space seperated). Makefile-relative path. -INC_DIR_LIST := inc \ - src \ +INC_DIR_LIST := include \ + /usr/lib/x86_64-linux-gnu/openmpi/include/ \ + src # Exclude files list(space seperated). Filenames only. # EXC_FILE_LIST := bad.cpp old.cpp @@ -77,8 +78,8 @@ DOCKER := CSIZE := size CFLAGS := $(DEB_CFLAGS) CXXFLAGS := $(DEB_CXXFLAGS) -CXX := g++ -CC := gcc +CXX := g++ #mpic++ +CC := gcc #mpicc # # =========== Main body and Patterns =========== @@ -180,11 +181,15 @@ hpc-clean: # make IMAGE=hpcimage v0 # make IMAGE=hpcimage v1_cilk # +dist_v05: CC := mpicc +dist_v05: CXX := mpic++ dist_v05: CFLAGS := $(DEB_CFLAGS) -DCODE_VERSION=50 dist_v05: CXXFLAGS := $(DEB_CXXFLAGS) -DCODE_VERSION=50 dist_v05: TARGET := dist_v05 dist_v05: $(BUILD_DIR)/$(TARGET) +dist_v1: CC := mpicc +dist_v1: CXX := mpic++ dist_v1: CFLAGS := $(DEB_CFLAGS) -DCODE_VERSION=100 dist_v1: CXXFLAGS := $(DEB_CXXFLAGS) -DCODE_VERSION=100 dist_v1: TARGET := dist_v1 diff --git a/homework_2/exersize.md b/homework_2/exersize.md index a61f4a1..b8bf998 100644 --- a/homework_2/exersize.md +++ b/homework_2/exersize.md @@ -1,4 +1,4 @@ -Parallel & Distributed Computer Systems +Parallel & Distributed Computer Systems HW2 December 6, 2024 diff --git a/homework_2/include/config.h b/homework_2/include/config.h new file mode 100644 index 0000000..433c74d --- /dev/null +++ b/homework_2/include/config.h @@ -0,0 +1,39 @@ +/*! + * \file config,h + * \brief Build configuration file. + * + * \author + * Christos Choutouridis AEM:8997 + * + */ + +#ifndef CONFIG_H_ +#define CONFIG_H_ + +#include +#include + +/* + * Defines for different version of the exercise + */ +#define V50 (50) +#define V100 (100) + + +// Fail-safe version selection +#if !defined CODE_VERSION +#define CODE_VERSION V1 +#endif + + +/*! + * Session option for each invocation of the executable + */ +struct session_t { + bool timing {false}; + bool verbose {false}; //!< Flag to enable verbose output to stdout +}; + +extern session_t session; + +#endif /* CONFIG_H_ */ diff --git a/homework_2/include/matrix.hpp b/homework_2/include/matrix.hpp new file mode 100644 index 0000000..33b3266 --- /dev/null +++ b/homework_2/include/matrix.hpp @@ -0,0 +1,804 @@ +/** + * \file matrix.hpp + * \brief A matrix abstraction implementation + * + * \author + * Christos Choutouridis AEM:8997 + * + */ +#ifndef MATRIX_HPP_ +#define MATRIX_HPP_ + +#include +#include +#include +#include +#include + +namespace mtx { + +using std::size_t; + +/* + * Small helper to strip types + */ +template +struct remove_cvref { + typedef std::remove_cv_t> type; +}; +template +using remove_cvref_t = typename remove_cvref::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 struct MatCol; +template struct MatRow; +template 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 +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 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 operator()(IndexType i, IndexType j) noexcept { + return MatVal(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 + void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) { + for (IndexType it=begin ; it(lambda)(std::forward(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 + 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 +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 = MatrixType::SPARSE; //!< meta:export of array type + static constexpr bool symmetric = Symmetric; //!< meta:export symmetric flag + + friend struct MatCol; + friend struct MatRow; + friend struct MatVal; + + /*! + * \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& row, const std::vector& 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& getValues() noexcept { return values; } + std::vector& getRows() noexcept { return rows; } + std::vector& 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 operator()(IndexType i, IndexType j) noexcept { + return MatVal(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 getCol(IndexType j) noexcept { + return MatCol(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 getRow(IndexType i) noexcept { + if constexpr (Symmetric) + return getCol(i); + else + return MatRow(this, i); + } + + // values only iterator support + DataType* begin() noexcept { return values.begin(); } + DataType* end() noexcept { return values.end(); } + //! @} + + //! A small iteration helper + template + void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) { + for (IndexType it=begin ; it(lambda)(std::forward(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 pair. + * index is the index of the item or end if not found + * status is true if found, false otherwise + */ + std::pair find_idx(const std::vector& 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 values {}; //!< vector to store the values of the matrix + std::vector rows{}; //!< vector to store the row information + std::vector 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 +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 + DataType operator* (C&& c) { + static_assert(std::is_same, MatCol>(), ""); + 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 +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 + DataType operator* (C&& c) { + static_assert(std::is_same, MatCol>(), ""); + 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 +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_ */ diff --git a/homework_2/include/utils.hpp b/homework_2/include/utils.hpp new file mode 100644 index 0000000..c6385ac --- /dev/null +++ b/homework_2/include/utils.hpp @@ -0,0 +1,87 @@ +/** + * \file utils.hpp + * \brief Utilities header + * + * \author + * Christos Choutouridis AEM:8997 + * + */ +#ifndef UTILS_HPP_ +#define UTILS_HPP_ + +#include +#include +#include + +#include "matrix.hpp" +#include "config.h" + +/*! + * A Logger for entire program. + */ +struct Log { + struct Endl {} endl; //!< a tag object to to use it as a new line request. + + //! We provide logging via << operator + template + Log& operator<< (T&& t) { + if (session.verbose) { + if (line_) { + std::cout << "[Log]: " << t; + line_ = false; + } + else + std::cout << t; + } + return *this; + } + // overload for special end line handling + Log& operator<< (Endl e) { (void)e; + if (session.verbose) { + std::cout << '\n'; + line_ = true; + } + return *this; + } +private: + bool line_ {true}; +}; + +extern Log logger; + +/*! + * A small timing utility based on chrono. + */ +struct Timing{ + using Tpoint = std::chrono::steady_clock::time_point; + using microseconds = std::chrono::microseconds; + using milliseconds = std::chrono::milliseconds; + using seconds = std::chrono::seconds; + + //! tool to mark the starting point + Tpoint start () noexcept { return start_ = std::chrono::steady_clock::now(); } + //! tool to mark the ending point + Tpoint stop () noexcept { return stop_ = std::chrono::steady_clock::now(); } + + auto dt () noexcept { + return std::chrono::duration_cast(stop_ - start_).count(); + } + //! tool to print the time interval + void print_dt (const char* what) noexcept { + if (session.timing) { + auto t = stop_ - start_; + if (std::chrono::duration_cast(t).count() < 10000) + std::cout << "[Timing]: " << what << ": " << std::to_string(std::chrono::duration_cast(t).count()) << " [usec]\n"; + else if (std::chrono::duration_cast(t).count() < 10000) + std::cout << "[Timing]: " << what << ": " << std::to_string(std::chrono::duration_cast(t).count()) << " [msec]\n"; + else + std::cout << "[Timing]: " << what << ": " << std::to_string(std::chrono::duration_cast(t).count()) << " [sec]\n"; + } + } +private: + Tpoint start_; + Tpoint stop_; +}; + + +#endif /* UTILS_HPP_ */ diff --git a/homework_2/src/main.cpp b/homework_2/src/main.cpp index a51103c..0bfa754 100644 --- a/homework_2/src/main.cpp +++ b/homework_2/src/main.cpp @@ -1,3 +1,78 @@ -int main(void) { +/*! + * \file main.cpp + * \brief Main application file for PDS HW2 (MPI) + * + * \author + * Christos Choutouridis AEM:8997 + * + */ +#include +#include + +#include + +#include "matrix.hpp" +#include "utils.hpp" +#include "config.h" + +// Global session data +session_t session; + +/*! + * A small command line argument parser + * \return The status of the operation + */ +bool get_options(int argc, char* argv[]){ + bool status =true; + + // iterate over the passed arguments + for (int i=1 ; i [-v]\n"; + std::cout << '\n'; + std::cout << "Options:\n\n"; + std::cout << " -v | --verbose\n"; + std::cout << " Request a more verbose output to stdout.\n\n"; + std::cout << " -h | --help\n"; + std::cout << " Prints this and exit.\n\n"; + std::cout << "Examples:\n\n"; + std::cout << " ...Example case...:\n"; + std::cout << " > distbitonic -x \n\n"; + + exit(0); + } + else { // parse error + std::cout << "Invocation error. Try -h for details.\n"; + status = false; + } + } + + return status; +} + +int main(int argc, char* argv[]) try { + // try to read command line + if (!get_options(argc, argv)) + exit(1); + return 0; +} +catch (std::exception& e) { + //we probably pollute the user's screen. Comment `cerr << ...` if you don't like it. + std::cerr << "Error: " << e.what() << '\n'; + exit(1); } \ No newline at end of file