Small changes to HW1

homework_1/inc/config.h

@@ -51,13 +51,13 @@ struct session_t {
    std::string   queryMtxFile   {};                        //!< optional query matrix file name in HDF5 format
    std::string   queryDataSet   {};                        //!< optional query dataset name in HDF5 matrix file
    bool          queryMtx       {false};                   //!< Flag to indicate that there is a separate query matrix
-   size_t        k              {1};                       //!< The number of nearest neighbors to find
+   size_t        k              {1UL};                     //!< The number of nearest neighbors to find
    std::string   outMtxFile     {"out.hdf5"};              //!< output matrix file name in HDF5 format
    std::string   outMtxIdxDataSet {"/Idx"};                //!< Index output dataset name in HDF5 matrix file
    std::string   outMtxDstDataSet {"/Dst"};                //!< Distance output dataset name in HDF5 matrix file
-   std::size_t   max_threads    {0};                       //!< Maximum threads to use
-   std::size_t   slices         {0};                       //!< Slices/threads to use
-   std::size_t   accuracy       {100};                     //!< The neighbor finding accuracy
+   std::size_t   max_threads    {0UL};                     //!< Maximum threads to use
+   std::size_t   slices         {0UL};                     //!< Slices/threads to use
+   std::size_t   accuracy       {100UL};                   //!< The neighbor finding accuracy
    bool          timing         {false};                   //!< Enable timing prints of the program
    bool          verbose        {false};                   //!< Flag to enable verbose output to stdout
 };

homework_1/inc/v0.hpp

@@ -39,7 +39,8 @@ void pdist2(const Matrix& X, const Matrix& Y, Matrix& D2) {
    int d = X.columns();
 
    // Compute the squared norms of each row in X and Y
-   std::vector<DataType> X_norms(M), Y_norms(N);
+   std::vector<DataType> X_norms(M);
+   std::vector<DataType> Y_norms(N);
    for (int i = 0; i < M ; ++i) {
       X_norms[i] = cblas_ddot(d, X.data() + i * d, 1, X.data() + i * d, 1);
    }
@@ -58,7 +59,6 @@ void pdist2(const Matrix& X, const Matrix& Y, Matrix& D2) {
          D2.set(std::sqrt(D2.get(i, j)),                i, j); // Take the square root of each
       }
    }
-   M++;
 }
 
 /*!
@@ -92,7 +92,7 @@ void quickselect(std::vector<std::pair<DataType, IndexType>>& vec, int k) {
  *             point of Q
  */
 template<typename MatrixD, typename MatrixI>
-void knnsearch(MatrixD& C, MatrixD& Q, size_t idx_offset, size_t k, size_t m, MatrixI& idx, MatrixD& dst) {
+void knnsearch(MatrixD& C, MatrixD& Q, size_t idx_offset, size_t k, [[maybe_unused]] size_t m, MatrixI& idx, MatrixD& dst) {
 
    using DstType = typename MatrixD::dataType;
    using IdxType = typename MatrixI::dataType;
@@ -104,10 +104,10 @@ void knnsearch(MatrixD& C, MatrixD& Q, size_t idx_offset, size_t k, size_t m, Ma
 
    pdist2(C, Q, D);
 
-   for (size_t j = 0; j < N; ++j) {
+   for (size_t j = 0UL; j < N; ++j) {
       // Create a vector of pairs (distance, index) for the j-th query
       std::vector<std::pair<DstType, IdxType>> dst_idx(M);
-      for (size_t i = 0; i < M; ++i) {
+      for (size_t i = 0UL; i < M; ++i) {
          dst_idx[i] = {D.data()[i * N + j], i};
       }
       // Find the k smallest distances using quickSelectKSmallest
@@ -117,7 +117,7 @@ void knnsearch(MatrixD& C, MatrixD& Q, size_t idx_offset, size_t k, size_t m, Ma
       std::sort(dst_idx.begin(), dst_idx.end());
 
       // Store the indices and distances
-      for (size_t i = 0; i < k; ++i) {
+      for (size_t i = 0UL; i < k; ++i) {
          dst.set(dst_idx[i].first, j, i);
          idx.set(dst_idx[i].second + idx_offset, j, i);
       }

homework_1/inc/v1.hpp

@@ -60,17 +60,17 @@ void mergeResultsWithM(mtx::Matrix<IndexType>& N1, mtx::Matrix<DataType>& D1,
                        size_t k, size_t m,
                        mtx::Matrix<IndexType>& N, mtx::Matrix<DataType>& D) {
    size_t numQueries = N1.rows();
-   size_t maxCandidates = std::min((IndexType)m, (IndexType)(N1.columns() + N2.columns()));
+   size_t maxCandidates = std::min(static_cast<IndexType>(m), static_cast<IndexType>(N1.columns() + N2.columns()));
 
-   for (size_t q = 0; q < numQueries; ++q) {
+   for (size_t q = 0UL; q < numQueries; ++q) {
       // Combine distances and neighbors
       std::vector<std::pair<DataType, IndexType>> candidates(N1.columns() + N2.columns());
 
       // Concatenate N1 and N2 rows
-      for (size_t i = 0; i < N1.columns(); ++i) {
+      for (size_t i = 0UL; i < N1.columns(); ++i) {
          candidates[i] = {D1.get(q, i), N1.get(q, i)};
       }
-      for (size_t i = 0; i < N2.columns(); ++i) {
+      for (size_t i = 0UL; i < N2.columns(); ++i) {
          candidates[i + N1.columns()] = {D2.get(q, i), N2.get(q, i)};
       }
 
@@ -81,7 +81,7 @@ void mergeResultsWithM(mtx::Matrix<IndexType>& N1, mtx::Matrix<DataType>& D1,
       std::sort(candidates.begin(), candidates.begin() + maxCandidates);
 
       // If m < k, pad the remaining slots with invalid values
-      for (size_t i = 0; i < k; ++i) {
+      for (size_t i = 0UL; i < k; ++i) {
          if (i < maxCandidates) {
             D.set(candidates[i].first, q, i);
             N.set(candidates[i].second, q, i);
@@ -110,7 +110,7 @@ void worker_body (std::vector<MatrixD>& corpus_slices,
    using IdxType = typename MatrixI::dataType;
 
 
-   for (size_t ci = 0; ci < num_slices; ++ci) {
+   for (size_t ci = 0UL; ci < num_slices; ++ci) {
       size_t idx_offset = ci * corpus_slice_size;
 
       // Intermediate matrixes for intermediate results
@@ -121,8 +121,8 @@ void worker_body (std::vector<MatrixD>& corpus_slices,
       v0::knnsearch(corpus_slices[ci], query_slices[slice], idx_offset, k, m, temp_idx, temp_dst);
 
       // Merge temporary results to final results
-      MatrixI idx_slice((IdxType*)idx.data(), slice * query_slice_size, query_slices[slice].rows(), k);
-      MatrixD dst_slice((DstType*)dst.data(), slice * query_slice_size, query_slices[slice].rows(), k);
+      MatrixI idx_slice(static_cast<IdxType*>(idx.data()), slice * query_slice_size, query_slices[slice].rows(), k);
+      MatrixD dst_slice(static_cast<DstType*>(dst.data()), slice * query_slice_size, query_slices[slice].rows(), k);
 
       mergeResultsWithM(idx_slice, dst_slice, temp_idx, temp_dst, k, m, idx_slice, dst_slice);
    }
@@ -145,29 +145,29 @@ void knnsearch(MatrixD& C, MatrixD& Q, size_t num_slices, size_t k, size_t m, Ma
     using IdxType = typename MatrixI::dataType;
 
     //Slice calculations
-    size_t corpus_slice_size = C.rows() / ((num_slices == 0)? 1:num_slices);
-    size_t query_slice_size = Q.rows() / ((num_slices == 0)? 1:num_slices);
+    size_t corpus_slice_size = C.rows() / ((num_slices == 0UL)? 1UL:num_slices);
+    size_t query_slice_size = Q.rows() / ((num_slices == 0UL)? 1UL:num_slices);
 
     // Make slices
-    std::vector<MatrixD> corpus_slices;
-    std::vector<MatrixD> query_slices;
+    std::vector<MatrixD> corpus_slices{};
+    std::vector<MatrixD> query_slices{};
 
-   for (size_t i = 0; i < num_slices; ++i) {
+   for (size_t i = 0UL; i < num_slices; ++i) {
       corpus_slices.emplace_back(
-         (DstType*)C.data(),
+         static_cast<DstType*>(C.data()),
          i * corpus_slice_size,
-         (i == num_slices - 1 ? C.rows() - i * corpus_slice_size : corpus_slice_size),
+         (i == num_slices - 1UL ? C.rows() - i * corpus_slice_size : corpus_slice_size),
          C.columns());
       query_slices.emplace_back(
-         (DstType*)Q.data(),
+              static_cast<DstType*>(Q.data()),
          i * query_slice_size,
-         (i == num_slices - 1 ? Q.rows() - i * query_slice_size : query_slice_size),
+         (i == num_slices - 1UL ? Q.rows() - i * query_slice_size : query_slice_size),
          Q.columns());
    }
 
    // Initialize results
-   for (size_t i = 0; i < dst.rows(); ++i) {
-      for (size_t j = 0; j < dst.columns(); ++j) {
+   for (size_t i = 0UL; i < dst.rows(); ++i) {
+      for (size_t j = 0UL; j < dst.columns(); ++j) {
          dst.set(std::numeric_limits<DstType>::infinity(), i, j);
          idx.set(static_cast<IdxType>(-1), i, j);
       }
@@ -176,7 +176,7 @@ void knnsearch(MatrixD& C, MatrixD& Q, size_t num_slices, size_t k, size_t m, Ma
     // Main loop
 #if defined OMP
    #pragma omp parallel for
-   for (size_t qi = 0; qi < num_slices; ++qi) {
+   for (size_t qi = 0UL; qi < num_slices; ++qi) {
          worker_body (corpus_slices, query_slices, idx, dst, qi, num_slices, corpus_slice_size, query_slice_size, k, m);
    }
 #elif defined CILK
@@ -185,6 +185,7 @@ void knnsearch(MatrixD& C, MatrixD& Q, size_t num_slices, size_t k, size_t m, Ma
    }
 #elif defined PTHREADS
    std::vector<std::thread> workers;
+   workers.reserve(num_slices);
 for (size_t qi = 0; qi < num_slices; ++qi) {
       workers.push_back(
          std::thread (worker_body<MatrixD, MatrixI>,

homework_2/.gitignore

@@ -1,23 +0,0 @@
-# project
-bin/
-out/
-mat/
-mtx/
-.unused/
-various/
-
-# hpc
-
-# IDEs
-.idea/
-.clangd
-
-# eclipse
-.project
-.cproject
-.settings/
-
-.vs/
-.vscode/
-
-

homework_2/Makefile

@@ -1,221 +0,0 @@
-#
-# PDS homework_1 Makefile
-#
-# Copyright (C) 2024 Christos Choutouridis <christos@choutouridis.net>
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as
-# published by the Free Software Foundation, either version 3
-# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-# ============== Project settings ==============
-# Project's name
-PROJECT         := PDS_homework_2
-# Excecutable's name
-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    := include \
-				   /usr/lib/x86_64-linux-gnu/openmpi/include/ \
-                   src
-
-# Exclude files list(space seperated). Filenames only.
-# EXC_FILE_LIST := bad.cpp old.cpp
-
-# Build directories
-BUILD_DIR       := bin
-OBJ_DIR         := $(BUILD_DIR)/obj
-DEP_DIR         := $(BUILD_DIR)/.dep
-
-
-# ========== Compiler settings ==========
-# Compiler flags for debug and release
-DEB_CFLAGS      := -DDEBUG -g3 -Wall -Wextra -std=c11
-REL_CFLAGS      := -Wall -Wextra -O3 -std=c11
-DEB_CXXFLAGS    := -DDEBUG -g3 -Wall -Wextra -std=c++17
-REL_CXXFLAGS    := -Wall -Wextra -O3 -std=c++17
-
-# Pre-defines
-# PRE_DEFS := MYCAB=1729 SUPER_MODE
-PRE_DEFS		:=
-
-# ============== Linker settings ==============
-# Linker flags (example: -pthread -lm)
-LDFLAGS         := -pthread
-
-# Map output file
-MAP_FILE        := output.map
-MAP_FLAG        := -Xlinker -Map=$(BUILD_DIR)/$(MAP_FILE)
-
-# ============== Docker settings ==============
-# We need:
-#  - Bind the entire project directory(the dir that icludes all the code) as volume.
-#  - In docker instance, change to working directory(where the makefile is).
-DOCKER_VOL_DIR  := $(shell pwd)
-DOCKER_WRK_DIR  :=
-DOCKER_RUN      := docker run --rm
-DOCKER_FLAGS    := -v $(DOCKER_VOL_DIR):/usr/src/$(PROJECT) -w /usr/src/$(PROJECT)/$(DOCKER_WRK_DIR)
-
-# docker invoke mechanism (edit with care)
-#   note:
-#   Here, `DOCKER` variable is empty. Rules can assign `DOCKER := DOCKER_CMD` when docker
-#   functionality is needed.
-DOCKER_CMD      = $(DOCKER_RUN) $(DOCKER_FLAGS) $(IMAGE)
-DOCKER          :=
-
-# ============== Tool selection ==============
-# compiler and compiler flags.
-CSIZE           := size
-CFLAGS          := $(DEB_CFLAGS)
-CXXFLAGS		:= $(DEB_CXXFLAGS)
-CXX             := g++ #mpic++
-CC              := gcc #mpicc
-
-#
-# =========== Main body and Patterns ===========
-#
-
-#ifeq ($(OS), Windows_NT)
-#	TARGET := $(TARGET).exe
-#endif
-INC     := $(foreach dir,$(INC_DIR_LIST),-I$(dir))
-DEF     := $(foreach def,$(PRE_DEFS),-D$(def))
-EXC     := $(foreach fil,$(EXC_FILE_LIST),                              \
-               $(foreach dir,$(SRC_DIR_LIST),$(wildcard $(dir)/$(fil))) \
-           )
-# source files. object and dependencies list
-# recursive search into current and source directories
-SRC     := $(wildcard *.cpp)
-SRC     += $(foreach dir,$(SRC_DIR_LIST),$(wildcard $(dir)/*.cpp))
-SRC     += $(foreach dir,$(SRC_DIR_LIST),$(wildcard $(dir)/**/*.cpp))
-SRC     := $(filter-out $(EXC),$(SRC))
-#SRC     := $(abspath $(SRC))
-
-OBJ     := $(foreach file,$(SRC:%.cpp=%.o),$(OBJ_DIR)/$(file))
-DEP     := $(foreach file,$(SRC:%.cpp=%.d),$(DEP_DIR)/$(file))
-
-
-# Make Dependencies pattern.
-# This little trick enables recompilation only when dependencies change
-# and it does so for changes both in source AND header files ;)
-# 
-# It is based on Tom Tromey's method.
-# 
-# Invoke cpp to create makefile rules with dependencies for each source file
-$(DEP_DIR)/%.d: %.c
-	@mkdir -p $(@D)
-	@$(DOCKER) $(CC) -E $(CFLAGS) $(INC) $(DEF) -MM -MT $(OBJ_DIR)/$(<:.c=.o) -MF $@ $<
-
-# c file objects depent on .c AND dependency files, which have an empty recipe 
-$(OBJ_DIR)/%.o: %.c $(DEP_DIR)/%.d
-	@mkdir -p $(@D)
-	@$(DOCKER) $(CC) -c $(CFLAGS) $(INC) $(DEF) -o $@ $<
-
-$(DEP_DIR)/%.d: %.cpp
-	@mkdir -p $(@D)
-	@$(DOCKER) $(CXX) -E $(CXXFLAGS) $(INC) $(DEF) -MM -MT $(OBJ_DIR)/$(<:.cpp=.o) -MF $@ $<
-
-# cpp file objects depent on .cpp AND dependency files, which have an empty recipe 
-$(OBJ_DIR)/%.o: %.cpp $(DEP_DIR)/%.d
-	@mkdir -p $(@D)
-	@$(DOCKER) $(CXX) -c $(CXXFLAGS) $(INC) $(DEF) -o $@ $<
-
-# empty recipe for dependency files. This prevents make errors
-$(DEP):
-
-# now include all dependencies
-# After all they are makefile dependency rules ;)
-include $(wildcard $(DEP))
-
-# main target rule
-$(BUILD_DIR)/$(TARGET): $(OBJ)
-	@mkdir -p $(@D)
-	@echo Linking to target: $(TARGET)
-	@echo $(DOCKER) $(CXX) '$$(OBJ)' $(LDFLAGS) $(MAP_FLAG) -o $(@D)/$(TARGET)
-	@$(DOCKER) $(CXX) $(OBJ) $(LDFLAGS) $(MAP_FLAG) -o $(@D)/$(TARGET)
-	@echo
-	@echo Print size information
-	@$(CSIZE) $(@D)/$(TARGET)
-	@echo Done
-
-.PHONY: clean
-clean:
-	@echo Cleaning build directories
-	@rm -rf $(OBJ_DIR)
-	@rm -rf $(DEP_DIR)
-	@rm -rf $(BUILD_DIR)
-
-#
-# ================ Local build rules =================
-# example:
-# make debug
-
-debug: CFLAGS := $(DEB_CFLAGS)
-debug: $(BUILD_DIR)/$(TARGET)
-
-release: CFLAGS := $(REL_CFLAGS)
-release: $(BUILD_DIR)/$(TARGET)
-
-all: release
-
-hpc-results/post:
-	$(CXX) $(CFLAGS) -o $@ hpc-results/main.cpp 
-
-hpc-clean:
-	rm hpc-results/post
-
-#
-# ================ Local (and/or) via docker build rules =================
-#
-# examples:
-# 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
-dist_v1: $(BUILD_DIR)/$(TARGET)
-
-#
-# ========= Inside CSAL Image build rules ===========
-#
-# 1) first jump into image (make sure you are in the directory where Makefile is):
-#    > docker run -it -v ${PWD}:/usr/src/exercise_1 -w /usr/src/exercise_1/ hpcimage
-# 2) Clean binaries first **important**
-#    > make clean
-# 3) for v4 cilk for example:
-#    > make csal_v4_cilk
-# 4) run executables from `bin/`. Examples:
-#    > ./bin/tcount_ompv3 -i mtx/NACA0015.mtx --timing -r 3 -o /dev/null
-#    > ./bin/tcount_pthv4 -i mtx/com_Youtube.mtx --timing --dynamic --print_count
-
-
-
-#
-# ======== Run from container =========
-#
-# examples:
-#
-# make IMAGE=hpcimage EXEC=knnsearch_v1 run
-# make IMAGE=hpcimage EXEC=knnsearch_v1 run
-#
-

homework_2/exersize.md

@@ -1,33 +0,0 @@
-Parallel & Distributed Computer Systems HW2
-
-December 6, 2024
-
-Write a distributed program that sorts $N$ integers in ascending order, using MPI. The inter-process communications should be defined by the Bitonic sort algorithm as presented in our class notes.
-
-The program must perform the following tasks:
-
-- The user specifies two positive integers $q$ and $p$.
-
-- Start $2^p$ processes with an array of $2^q$ random integers is each processes.
-
-- Sort all $N = \left ( 2^{(q + p)} \right)$ elements int ascending order.
-
-- Check the correctness of the final result.
-
-Your implementation should be based on the following steps:
-- Start p processes, each process gets n/p data and sorts (ascending or descending depending on the process id) them using sort from any library. Use two buffers, to sort from one to the other, to send from one and receive in the other.
-- Repeat $O((log_2(p))^2)$:
-    * Exchange data with the corresponding partner and keep the min or max elements, depending on the process id and phase of the computation
-    * Sort locally the bitonic sequence using a modification of merge sort, starting from the "elbow" of the bitonic.
-
-You may use the C standard library function stdlib qsort() to check the correctness of your results and as the initial local sorting routine for each process.
-
-You must deliver:
-- A report (about $3-4$ pages) that describes your parallel algorithm and implementation.
-
-- Your comments on the speed of your parallel program compared to the serial sort, after trying you program on aristotelis for $p = [1:7]$ and $q = [20:27]$.
-
-- The source code of your program uploaded online.
-
-Ethics: If you use code found on the web or by an LLM, you should mention your source and the changes you made. You may work in pairs; both partners must submit a single report with both names.
-Deadline: 7 January, $2025$.

homework_2/include/config.h

@@ -1,39 +0,0 @@
-/*!
- * \file    config,h
- * \brief   Build configuration file.
- *
- * \author
- *    Christos Choutouridis AEM:8997
- *    <cchoutou@ece.auth.gr>
- */
-
-#ifndef CONFIG_H_
-#define CONFIG_H_
-
-#include <iostream>
-#include <string>
-
-/*
- * 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_ */

homework_2/include/matrix.hpp

@@ -1,804 +0,0 @@
-/**
- * \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_ */

homework_2/include/utils.hpp

@@ -1,87 +0,0 @@
-/**
- * \file    utils.hpp
- * \brief   Utilities header
- *
- * \author
- *    Christos Choutouridis AEM:8997
- *    <cchoutou@ece.auth.gr>
- */
-#ifndef UTILS_HPP_
-#define UTILS_HPP_
-
-#include <iostream>
-#include <chrono>
-#include <unistd.h>
-
-#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<typename T>
-   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<std::chrono::microseconds>(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<microseconds>(t).count() < 10000)
-            std::cout << "[Timing]: " << what << ": " << std::to_string(std::chrono::duration_cast<microseconds>(t).count()) << " [usec]\n";
-         else if (std::chrono::duration_cast<milliseconds>(t).count() < 10000)
-            std::cout << "[Timing]: " << what << ": " << std::to_string(std::chrono::duration_cast<milliseconds>(t).count()) << " [msec]\n";
-         else
-            std::cout << "[Timing]: " << what << ": " << std::to_string(std::chrono::duration_cast<seconds>(t).count()) << " [sec]\n";
-      }
-   }
-private:
-   Tpoint start_;
-   Tpoint stop_;
-};
-
-
-#endif /* UTILS_HPP_ */

homework_2/julia/bitonic_v05.jl

@@ -1,114 +0,0 @@
-#
-# ---------------------------------------------
-# Bitonic v0.5 functionality
-#
-
-function partner(node, step)
-    partner = node + ((((node+step) % 2) == 0) ? 1 : -1)
-end
-
-function active(id, p)
-    ret =  (id >= 0) && (id < p)
-end
-
-function exchange(localid,  remoteid)
-    if verbose
-        println("Exchange local data from $localid with partner $remoteid")
-    end
-    nothing # We have all data here ;)
-end
-
-function minmax(data, localid, remoteid, keepsmall)
-    # Keep min-max on local data
-    temp = copy(data[localid+1, :])
-    if (keepsmall)
-        view(data, localid+1, :)  .= min.(temp, data[remoteid+1, :])
-        view(data, remoteid+1, :) .= max.(temp, data[remoteid+1, :])
-    else
-        view(data, localid+1, :)  .= max.(temp, data[remoteid+1, :])
-        view(data, remoteid+1, :) .= min.(temp, data[remoteid+1, :])
-    end
-end
-
-"""
-    distbubletonic!(p, data)
-
-distributed bitonic v0.5 sort using a "bubble sort"-like functionality to propagate large and small
-items between nodes.
-
-p:    The number of processes
-data: (p, N/p) array
-"""
-function distbubletonic!(p, data) 
-  
-    pid = 0:p-1
-    ascending = mod.(pid,2) .== 0
-    if verbose
-        println("ascending: $ascending")
-    end
-    # local full sort here (run all MPI nodes)
-    for i in 1:p
-        sort!(view(data, i, :), rev = !ascending[i])
-    end
-    for step in 0:p-2
-        direction = [true for x = 1:p]
-        partnerid = partner.(pid, step)
-        activeids = active.(partnerid, p)
-        keepsmall = pid .< partnerid
-        if verbose
-            println("step: $step | active ids: $activeids | partner: $partnerid | keepsmall: $keepsmall")
-        end
-        # exchange with partner and keep small or large (run all MPI nodes)
-        for i in 0:p-1
-            l_idx = i+1
-            r_idx = partnerid[i+1]+1
-            if activeids[l_idx] && i < partnerid[l_idx]
-                exchange(i, partnerid[l_idx])
-                minmax(data, i, partnerid[l_idx], keepsmall[l_idx])
-                sort!(view(data, l_idx, :), rev = !ascending[l_idx])    # elbow sort here
-                sort!(view(data, r_idx, :), rev = !ascending[r_idx])    # elbow sort here
-            end
-        end
-    end
-    # [optional] reverse the odd positions (run all MPI nodes)
-    for i in 1:p
-        if !ascending[i]
-            sort!(view(data, i, :))
-        end
-    end
-    nothing
-end
-
-
-
-#
-# Homework setup
-# ---------------------------------------------
-#
-p::Int8 = 3  # The order of number of "processors"
-q::Int8 = 8  # The data size order (power of 2) of each "processor"
-verbose = false;
-
-
-# Run Script
-# ---------------------------------------------
-P::Int = 2^p
-Q::Int = 2^q
-N::Int = 2^(q+p)
-
-println("Distributed Bubbletonic (v0.5) test")
-println("p: $p -> Number of processors: $P")
-println("q: $q -> Data length for each node: $Q, Total: $(P*Q)")
-
-println("Create an $P x $Q array")
-Data = rand(Int8, P, Q)
-
-println("Sort array with $P (MPI) nodes")
-@time distbubletonic!(P, Data)
-
-# Test
-if issorted(vec(permutedims(Data)))
-    println("Test: Passed")
-else
-    println("Test: Failed")
-end

homework_2/julia/bitonic_v1.jl

@@ -1,113 +0,0 @@
-#
-# ---------------------------------------------
-# Bitonic v0.5 functionality
-#
-
-function exchange(localid,  remoteid)
-    if verbose
-        println("Exchange local data from $localid with partner $remoteid")
-    end
-    nothing # We have all data here ;)
-end
-
-function minmax(data, localid, remoteid, keepsmall)
-    # Keep min-max on local data
-    temp = copy(data[localid+1, :])
-    if (keepsmall)
-        view(data, localid+1, :)  .= min.(temp, data[remoteid+1, :])
-        view(data, remoteid+1, :) .= max.(temp, data[remoteid+1, :])
-    else
-        view(data, localid+1, :)  .= max.(temp, data[remoteid+1, :])
-        view(data, remoteid+1, :) .= min.(temp, data[remoteid+1, :])
-    end
-end
-
-
-
-function sort_network!(data, n, depth)
-    nodes = 0:n-1
-    for step = depth-1:-1:0
-        partnerid = nodes .⊻ (1 << step)
-        direction = (nodes .& (1 << depth)) .== 0 .& (nodes .< partnerid)
-        keepsmall = ((nodes .< partnerid) .& direction) .| ((nodes .> partnerid) .& .!direction)
-        if verbose
-            println("depth: $depth | step: $step | partner: $partnerid | keepsmall: $keepsmall")
-        end
-        # exchange with partner and keep small or large (run all MPI nodes)
-        for i in 0:n-1
-            if (i < partnerid[i+1])
-                exchange(i, partnerid[i+1])
-                minmax(data, i, partnerid[i+1], keepsmall[i+1])
-            end
-        end
-    end
-end
-
-
-
-"""
-    distbitonic!(p, data)
-
-distributed bitonic sort v1 using elbow merge locally except for the first step
-p:    The number of processes
-data: (p, N/p) array
-"""
-function distbitonic!(p, data) 
-
-    q = Int(log2(p))    # CPU order
-  
-    pid = 0:p-1
-    ascending = mod.(pid,2) .== 0
-    if verbose
-        println("ascending: $ascending")
-    end
-    # local full sort here (run all MPI nodes)
-    for i in 1:p
-        sort!(view(data, i, :), rev = !ascending[i])
-    end
-    for depth = 1:q
-        sort_network!(data, p, depth)
-        ascending = (pid .& (1 << depth)) .== 0
-        if verbose
-            println("ascending: $ascending")
-        end
-        # local elbowmerge here (run all MPI nodes)
-        for i in 1:p
-            sort!(view(data, i, :), rev = !ascending[i])
-        end
-    end
-  
-    nothing
-end
-
-#
-# Homework setup
-# ---------------------------------------------
-#
-p::Int8 = 3  # The order of number of "processors"
-q::Int8 = 8  # The data size order (power of 2) of each "processor"
-verbose = false;
-
-
-# Run Script
-# ---------------------------------------------
-P::Int = 2^p
-Q::Int = 2^q
-N::Int = 2^(q+p)
-
-println("Distributed bitonic (v1) test")
-println("p: $p -> Number of processors: $P")
-println("q: $q -> Data length for each node: $Q, Total: $(P*Q)")
-
-println("Create an $P x $Q array")
-Data = rand(Int8, P, Q)
-
-println("Sort array with $P (MPI) nodes")
-@time distbitonic!(P, Data)
-
-# Test
-if issorted(vec(permutedims(Data)))
-    println("Test: Passed")
-else
-    println("Test: Failed")
-end

homework_2/julia/distbitonic.jl

@@ -1,27 +0,0 @@
-# distributed bitonic sort using elbow merge locally except for the first step
-function distbitonic!(p)
-
-    q = Int(log2(p))
-  
-    pid = 0:p-1
-    ascending = mod.(pid,2) .== 0
-    println("ascending: $ascending")
-    # local full sort here
-    for k = 1:q
-      kk = 1 << k
-      for j = k-1:-1:0
-        jj = 1 << j
-        partnerid = pid .⊻ jj
-        direction = (pid .& kk) .== 0 .& (pid .< partnerid)
-        keepsmall = ((pid .< partnerid) .& direction) .| ((pid .> partnerid) .& .!direction)
-        println("k: $k | j: $j | partner: $partnerid | keepsmall: $keepsmall")
-        # exchange with partner and keep small or large
-      end
-      ascending = (pid .& kk) .== 0
-      println("ascending: $ascending")
-      # local elbowmerge here
-    end
-  
-    nothing
-  end
- 

homework_2/julia/elbowmerge.jl

@@ -1,23 +0,0 @@
-
-# given a bitonic sequence b, merge it into a sorted sequence s
-@inbounds function elbowmerge!(s, b)
-    n = length(b)
-    l = argmin(b)
-    r = l == n ? 1 : l + 1
-  
-    i = 1
-    while i <= n
-      if b[l] < b[r]
-        s[i] = b[l]
-        l = l == 1 ? n : l - 1
-      else
-        s[i] = b[r]
-        r = r == n ? 1 : r + 1
-      end
-      i += 1
-    end
-    nothing
-  end
-  
-
-  

homework_2/src/main.cpp

@@ -1,78 +0,0 @@
-/*!
- * \file    main.cpp
- * \brief   Main application file for PDS HW2 (MPI)
- *
- * \author
- *    Christos Choutouridis AEM:8997
- *    <cchoutou@ece.auth.gr>
- */
-#include <exception>
-#include <iostream>
-
-#include <mpi.h>
-
-#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<argc ; ++i) {
-        std::string arg(argv[i]);     // get current argument
-
-        if (arg == "-x" || arg == "--xxxxx") {
-            if (i+2 < argc) {
-//                session.corpusMtxFile = std::string(argv[++i]);
-//                session.corpusDataSet = std::string(argv[++i]);
-            }
-            else
-                status = false;
-        }
-
-        else if (arg == "-v" || arg == "--verbose")
-            session.verbose = true;
-        else if (arg == "-h" || arg == "--help") {
-            std::cout << "distbitonic - A distributed bitonic sort\n\n";
-            std::cout << "distbitonic -x <> [-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 <xxxxx>  \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);
-}