Init commit: a bare implementation before parallelism

.gitignore

@@ -0,0 +1,9 @@
+# project
+resources/
+bin/
+
+# eclipse
+.project
+.cproject
+.settings/
+

Makefile

@@ -0,0 +1,167 @@
+#
+# PDS excercise 1 Makefile
+#
+# Copyright (C) 2019-2020 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         := exercise_1
+# Excecutable's name
+TARGET          := triangCounting
+# 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    := src inc
+# 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=c++14
+REL_CFLAGS      := -Wall -Wextra -O2 -std=c++14
+# Pre-defines
+# PRE_DEFS := MYCAB=1729 SUPER_MODE
+
+# ============== 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  := $${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)
+CXX             := g++
+
+
+#
+# =========== 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: %.cpp
+	@mkdir -p $(@D)
+	@$(DOCKER) $(CXX) -E $(CFLAGS) $(INC) $(DEF) -MM -MT $(OBJ_DIR)/$(<:.cpp=.o) -MF $@ $<
+
+# 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 $(CFLAGS) $(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) $(LDFLAGS) $(MAP_FLAG) -o $(@D)/$(TARGET) '$$(OBJ)'
+	@$(DOCKER) $(CXX) $(LDFLAGS) $(MAP_FLAG) -o $(@D)/$(TARGET) $(OBJ)
+	@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 =================
+# examples:
+# make debug
+
+.PHONY: debug
+debug: CFLAGS := $(DEB_CFLAGS)
+debug: $(BUILD_DIR)/$(TARGET)
+
+.PHONY: release
+release: CFLAGS := $(REL_CFLAGS)
+release: $(BUILD_DIR)/$(TARGET)
+
+.PHONY: all
+all: release
+
+#
+# ================ Docker based rules ================
+# examples:
+# make IMAGE="gcc:8.3" dock
+#
+.PHONY: dock
+dock: DOCKER := $(DOCKER_CMD)
+dock: CFLAGS := $(REL_CFLAGS)
+dock: $(BUILD_DIR)/$(TARGET)
+
+

inc/config.h

@@ -0,0 +1,36 @@
+/*!
+ * \file    config,h
+ * \brief   Build configuration file.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+
+#ifndef CONFIG_H_
+#define CONFIG_H_
+
+#include <v12.h>
+#include <v3.h>
+#include <v4.h>
+
+#define  V12   2
+#define  V3    3
+#define  V4    4
+
+#if !defined CODE_VERSION
+ #define CODE_VERSION   V4
+#endif
+
+#if   CODE_VERSION == V12
+   using namespace v12;
+   using matrix = v12::matrix;
+#elif CODE_VERSION == V3
+   using namespace v3;
+   using matrix = v3::matrix;
+#elif CODE_VERSION == V4
+   using namespace v4;
+   using matrix = v4::matrix;
+#endif
+
+#endif /* CONFIG_H_ */

inc/impl.hpp

@@ -0,0 +1,519 @@
+/**
+ * \file    impl.hpp
+ * \brief   Implementation common header file
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+#ifndef IMPL_HPP_
+#define IMPL_HPP_
+
+#include <type_traits>
+#include <utility>
+#include <algorithm>
+#include <vector>
+#include <tuple>
+#include <cstddef>
+#include <iostream>
+#include <fstream>
+
+using std::size_t;
+
+/*
+ * Small helper to strip types
+ */
+template<typename T>
+struct remove_cvref {
+    typedef std::remove_cv_t<std::remove_reference_t<T>> type;
+};
+template<typename T>
+using remove_cvref_t = typename remove_cvref<T>::type;
+
+/*!
+ * Enumerator to denote the storage type of the array to use.
+ */
+enum class MatrixType {
+   FULL,          /*!< Matrix is asymmetric */
+   SYMMETRIC,     /*!< Matrix is symmetric */
+};
+
+/*
+ * Forward type declerations
+ */
+template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC> struct Matrix;
+template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC> struct SpMat;
+template<typename DataType, typename IndexType> struct SpMatCol;
+template<typename DataType, typename IndexType> struct SpMatRow;
+template<typename DataType, typename IndexType> struct SpMatVal;
+
+ /*!
+ * 2D-array wrapper for v1 and v2 part of the exercise t use as RAII
+ * and copy-prevention.
+ *
+ * This is a very thin abstraction layer over a native array.
+ * This is tested using compiler explorer and our template produce
+ * almost identical assembly.
+ *
+ * The penalty hit we have is due to the fact that we use a one dimension array
+ * and we have to calculate the actual position from an (i,j) pair.
+ * The use of 1D array was our intention from the beginning, so the penalty
+ * was pretty much unavoidable.
+ *
+ * \tparam DataType  The underling data type of the array
+ * \tparam Type      The storage type of the array
+ *    \arg  FULL              For full matrix
+ *    \arg  SYMMETRIC         For symmetric matrix (we use only the lower part)
+ */
+template<typename DataType, typename IndexType, MatrixType Type>
+struct Matrix {
+
+   using dataType    = DataType;                   //!< meta:export of underling data type
+   using indexType   = IndexType;                  //!< meta:export of underling index type
+   static constexpr MatrixType matrixType = Type;  //!< export of array type
+
+   /*!
+    * \name Obj lifetime
+    */
+   //! @{
+
+   //! Constructor using data type and size
+   Matrix(DataType* t, IndexType s) noexcept : m_(t), size_(s) { }
+
+   //! RAII deleter
+   ~Matrix() { delete m_; }
+   //! move ctor
+   Matrix(Matrix&& a) noexcept { a.swap(*this); }
+   //! move
+   Matrix& operator=(Matrix&& a) noexcept { a.swap(*this); return *this; }
+   Matrix(const Matrix& a)             = delete;  //!< No copy ctor
+   Matrix& operator=(const Matrix& a)  = delete;  //!< No copy
+   //! @}
+
+   //! \name Data exposure
+   //! @{
+
+   //! memory capacity of the matrix with diagonal data
+   template<MatrixType AT=Type> std::enable_if_t<AT==MatrixType::SYMMETRIC, IndexType>
+      static constexpr capacity(IndexType N) noexcept { return (N+1)*N/2; }
+   //! memory capacity for full matrix
+   template<MatrixType AT=Type> std::enable_if_t<AT==MatrixType::FULL, IndexType>
+      static constexpr capacity(IndexType N) noexcept { return N*N; }
+
+   //! Get the size of each dimension
+   IndexType size() noexcept { return size_; }
+
+   /*
+    * virtual 2D accessors
+    */
+   template<MatrixType AT=Type>
+   std::enable_if_t <AT==MatrixType::SYMMETRIC, DataType>
+   get (IndexType i, IndexType j) {
+      if (i<j)    return m_[j*(j+1)/2 + i];  // Upper, use opposite index
+      else        return m_[i*(i+1)/2 + j];  // Lower, use our notation
+   }
+   template<MatrixType AT=Type>
+   std::enable_if_t <AT==MatrixType::FULL, DataType>
+   get(IndexType i, IndexType j) {
+      return m_[i*size_ + j];
+   }
+   template<MatrixType AT=Type>
+   std::enable_if_t <AT==MatrixType::SYMMETRIC, DataType>
+   set(DataType v, IndexType i, IndexType j) {
+      if (i<j)    return m_[j*(j+1)/2 + i] =v;  // Upper, use opposite index
+      else        return m_[i*(i+1)/2 + j] =v;  // Lower, use our notation
+   }
+   template<MatrixType AT=Type>
+   std::enable_if_t <AT==MatrixType::FULL, DataType>
+   set(DataType v, IndexType i, IndexType j) {
+      return m_[i*size_ + j] =v;
+   }
+   DataType operator()(IndexType i, IndexType j) { return get(i, j); }
+
+   // a basic serial iterator support
+   DataType* begin() noexcept { return m_; }
+   DataType* end()   noexcept { return m_ + Matrix::capacity(size_); }
+   //! @}
+
+   /*!
+    * \name Safe iteration API
+    *
+    * This api automates the iteration over the array based on
+    * MatrixType
+    */
+   //! @{
+   template<typename F, typename... Args>
+   void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) {
+      for (IndexType it=begin ; it<end ; ++it) {
+         std::forward<F>(lambda)(std::forward<Args>(args)..., it);
+      }
+   }
+   //! @}
+
+   // move helper
+   void swap(Matrix& src) noexcept {
+      std::swap(m_, src.m_);
+      std::swap(size_, src.size_);
+   }
+private:
+   DataType*   m_;      //!< Pointer to actual data.
+   IndexType   size_;   //!< the virtual size of each dimension.
+};
+
+ /*!
+ * RAII allocation helper, the smart_ptr way.
+ */
+template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC>
+Matrix<DataType, IndexType, Type> make_Matrix(IndexType s) {
+   return Matrix<DataType, IndexType, Type>(new DataType[Matrix<DataType, IndexType, Type>::capacity(s)], s);
+}
+
+template<typename DataType, typename IndexType>
+using CooVal = std::tuple<DataType, IndexType, IndexType>;
+
+template<typename DataType, typename IndexType, MatrixType Type>
+struct SpMat {
+   using dataType    = DataType;                   //!< meta:export of underling data type
+   using indexType   = IndexType;                  //!< meta:export of underling index type
+   static constexpr MatrixType matrixType = Type;  //!< export of array type
+
+   friend class SpMatCol<DataType, IndexType>;
+   friend class SpMatRow<DataType, IndexType>;
+   friend class SpMatVal<DataType, IndexType>;
+
+   /*!
+    * \name Obj lifetime
+    */
+   //! @{
+
+   //! allocation with init value ctor
+   SpMat(IndexType n=IndexType{}, IndexType nnz=IndexType{}) :
+      values(nnz, DataType{}),
+      rows(nnz, IndexType{}),
+      col_ptr((n)? n+1:2, IndexType{}),
+      N(n),
+      NNZ(nnz) { }
+
+   SpMat(IndexType n, IndexType nnz, const IndexType* row, const IndexType* col) :
+      values(nnz, 1),
+      rows(row, row+nnz),
+      col_ptr(col, col+n+1),
+      N(n),
+      NNZ(nnz) { }
+
+   SpMat(IndexType n, IndexType nnz, const DataType* v, const IndexType* row, const IndexType* col) :
+      values(v, v+nnz),
+      rows(row, row+nnz),
+      col_ptr(col, col+n+1),
+      N(n),
+      NNZ(nnz) { }
+
+   SpMat(IndexType n, IndexType nnz, const DataType v, const std::vector<IndexType>& row, const std::vector<IndexType>& col) :
+      values(nnz, v),
+      rows (row),
+      col_ptr(col),
+      N(n),
+      NNZ(nnz) { }
+
+   //! move ctor
+   SpMat(SpMat&& a) noexcept { moves(std::move(a)); }
+   //! move
+   SpMat& operator=(SpMat&& a) noexcept { moves(std::move(a)); return *this; }
+   SpMat(const SpMat& a)             = delete;  //!< No copy ctor
+   SpMat& operator=(const SpMat& a)  = delete;  //!< No copy assignment
+   //! @}
+
+   //! \name Data exposure
+   //! @{
+   IndexType size()     noexcept { return N; }
+   IndexType size(IndexType n) {
+      col_ptr.resize(n+1);
+      return N = n;
+   }
+   IndexType capacity() noexcept { return NNZ; }
+   IndexType capacity(IndexType nnz) {
+      values.reserve(nnz);
+      rows.reserve(nnz);
+      return NNZ;
+   }
+   std::vector<IndexType>& getRows() noexcept { return rows; }
+   std::vector<IndexType>& getCols() noexcept { return col_ptr; }
+
+   SpMatVal<DataType, IndexType> operator()(IndexType i, IndexType j) {
+      return SpMatVal<DataType, IndexType>(this, get(i, j), i, j);
+   }
+   DataType get(IndexType i, IndexType j) {
+      IndexType idx;
+      bool found;
+      std::tie(idx, found) =find_idx(rows, col_ptr[j], col_ptr[j+1], i);
+      return (found) ? values[idx] : 0;
+   }
+   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;
+      }
+   }
+
+
+   SpMatCol<DataType, IndexType> getCol(IndexType j) {
+      return SpMatCol<DataType, IndexType>(this, col_ptr[j], col_ptr[j+1]);
+   }
+   template<MatrixType AT= Type>
+   std::enable_if_t<AT==MatrixType::SYMMETRIC, SpMatCol<DataType, IndexType>>
+   getRow(IndexType i) {
+      return getCol(i);
+   }
+   template<MatrixType AT= Type>
+   std::enable_if_t<AT==MatrixType::FULL, SpMatCol<DataType, IndexType>>
+   getRow(IndexType i) {
+      return SpMatRow<DataType, IndexType>(this, i);
+   }
+
+   // iterator support
+   DataType* begin() noexcept { return values.begin(); }
+   DataType* end()   noexcept { return values.end(); }
+   //! @}
+
+   /*!
+    * \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);
+      }
+   }
+
+   //! @}
+
+   // operations
+   template<typename D, typename I, MatrixType T> friend void print(SpMat<D, I, T>& mat);
+   template<typename D, typename I, MatrixType T> friend void print_dense(SpMat<D, I, T>& mat);
+
+private:
+   // index-find helper
+   std::pair<IndexType, bool> find_idx(const std::vector<IndexType>& v, IndexType begin, IndexType end, IndexType match) {
+      for ( ; begin < end ; ++begin) {
+         if       (match == v[begin])   return std::make_pair(begin, true);
+         else if  (match <  v[begin])   return std::make_pair(begin, false);
+      }
+      return std::make_pair(end, false);
+   }
+   // move helper
+   void moves(SpMat&& src) noexcept {
+      values  = std::move(src.values);
+      rows    = std::move(src.rows);
+      col_ptr = std::move(src.col_ptr);
+      N       = std::move(src.N);   // redundant for primitives
+      NNZ     = std::move(src.NNZ); //
+   }
+   //! \name Data
+   //! @{
+   std::vector<DataType>   values {};
+   std::vector<IndexType>  rows{};
+   std::vector<IndexType>  col_ptr{1,0};
+   IndexType   N{0};
+   IndexType   NNZ{0};
+   //! @}
+};
+
+template<typename DataType, typename IndexType>
+struct SpMatCol {
+   using owner_t = SpMat<DataType, IndexType>;
+
+   SpMatCol(owner_t* own, const IndexType begin, const IndexType end) noexcept :
+      owner_(own), index_(begin), begin_(begin), end_(end) {
+      vindex_ = vIndexCalc(index_);
+   }
+   SpMatCol()                          = default;
+   SpMatCol(const SpMatCol&)           = delete;
+   SpMatCol& operator=(const SpMatCol&)= delete;
+   SpMatCol(SpMatCol&&)                = default;
+   SpMatCol& operator=(SpMatCol&&)     = default;
+
+   DataType operator* () {
+      return get();
+   }
+   SpMatCol& operator++ ()    { advance(); return *this; }
+   SpMatCol& operator++ (int) { SpMatCol& p = *this; advance(); return p; }
+
+   DataType operator()(IndexType x) {
+      return (x == index())? get() : DataType{};
+   }
+   DataType operator= (DataType v) { return owner_->values[index_] = v; }
+   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; }
+
+   template <typename C>
+   DataType operator* (C&& c) {
+      static_assert(std::is_same<remove_cvref_t<C>, SpMatCol<DataType, IndexType>>(), "");
+      DataType v{};
+      while (index() != end() && c.index() != c.end()) {
+         if      (index() < c.index())  advance();
+         else if (index() > c.index())  ++c;
+         else { //index() == c.index()
+            v += get() * *c;
+            ++c;
+            advance();
+         }
+      }
+      return v;
+   }
+
+private:
+   void advance() noexcept {
+      ++index_;
+      vindex_ = vIndexCalc(index_);
+   }
+   IndexType vIndexCalc(IndexType idx) {
+      return (idx < end_) ? owner_->rows[idx] : end();
+   }
+   DataType get() { return owner_->values[index_]; }
+   owner_t*    owner_{nullptr};
+   IndexType   vindex_  {IndexType{}};
+   IndexType   index_{IndexType{}};
+   IndexType   begin_{IndexType{}};
+   IndexType   end_{IndexType{}};
+};
+
+template<typename DataType, typename IndexType>
+struct SpMatRow {
+   using owner_t = SpMat<DataType, IndexType>;
+
+   SpMatRow(owner_t* own, const IndexType row) noexcept :
+      owner_(own), vindex_(IndexType{}), row_(row), index_(IndexType{}),
+      begin_(IndexType{}), end_(owner_->NNZ) {
+      // place begin
+      while(begin_ != end_ && owner_->rows[begin_] != row_)
+         ++begin_;
+      // place index_ and vindex_
+      if (owner_->rows[index_] != row_)
+         advance();
+   }
+   SpMatRow()                          = default;
+   SpMatRow(const SpMatRow&)           = delete;
+   SpMatRow& operator=(const SpMatRow&)= delete;
+   SpMatRow(SpMatRow&&)                = default;
+   SpMatRow& operator=(SpMatRow&&)     = default;
+
+   DataType operator* () {
+      return get();
+   }
+   SpMatRow& operator++ ()    { advance(); return *this; }
+   SpMatRow& operator++ (int) { SpMatRow& p = *this; advance(); return p; }
+   DataType operator()(IndexType x) {
+      return (x == index())? get() : DataType{};
+   }
+   DataType operator= (DataType v) { return owner_->values[index_] = v; }
+   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; }
+
+   template <typename C>
+   DataType operator* (C&& c) {
+      static_assert(std::is_same<remove_cvref_t<C>, SpMatCol<DataType, IndexType>>(), "");
+      DataType v{};
+      while (index() != end() && c.index() != c.end()) {
+         if      (index() < c.index())  advance();
+         else if (index() > c.index())  ++c;
+         else { //index() == c.index()
+            v += get()* *c;
+            ++c;
+            advance();
+         }
+      }
+      return v;
+   }
+private:
+   void advance() noexcept {
+      do
+         ++index_;
+      while(index_ != end_ && owner_->rows[index_] != row_);
+      vindex_ = vIndexCalc(index_);
+   }
+   IndexType vIndexCalc(IndexType idx) {
+      for(IndexType i =0 ; i<(owner_->N+1) ; ++i)
+         if (idx < owner_->col_ptr[i])
+            return i-1;
+      return end();
+   }
+   DataType get() { return owner_->values[index_]; }
+   owner_t*    owner_   {nullptr};
+   IndexType   vindex_  {IndexType{}};
+   IndexType   row_     {IndexType{}};
+   IndexType   index_   {IndexType{}};
+   IndexType   begin_   {IndexType{}};
+   IndexType   end_     {IndexType{}};
+};
+
+template<typename DataType, typename IndexType>
+struct SpMatVal {
+   using owner_t = SpMat<DataType, IndexType>;
+
+   SpMatVal(owner_t* own, DataType v, IndexType i, IndexType j) :
+      owner_(own), v_(v), i_(i), j_(j) { }
+   SpMatVal()                           = default;
+   SpMatVal(const SpMatVal&)            = delete;
+   SpMatVal& operator=(const SpMatVal&) = delete;
+   SpMatVal(SpMatVal&&)                 = default;
+   SpMatVal& operator=(SpMatVal&&)      = default;
+
+   operator DataType() { return v_; }
+   SpMatVal& operator=(DataType v) {
+      v_ = v;
+      owner_->set(v_, i_, j_);
+      return *this;
+   }
+private:
+   owner_t*    owner_{nullptr};;
+   DataType    v_{DataType{}};
+   IndexType   i_{IndexType{}};
+   IndexType   j_{IndexType{}};
+};
+
+
+enum class InputMatrix{
+   GENERATE,
+   MTX
+};
+
+/*!
+ * Session option for each invocation of the executable
+ */
+struct session_t {
+   std::size_t    size           {0};
+   double         probability    {0};
+   InputMatrix    inputMatrix    {InputMatrix::GENERATE};
+   std::ifstream  mtxFile        {};
+   std::size_t    print_size     {80};
+   bool           timing         {false};
+   bool           print          {false};
+   bool           makeSymmetric  {false};
+};
+
+extern session_t session;
+
+
+#endif /* IMPL_HPP_ */

inc/utils.h

@@ -0,0 +1,198 @@
+/*!
+ * \file    utils.h
+ * \brief   Utilities to handle matrix files, chrono, etc...
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+#ifndef UTILS_H_
+#define UTILS_H_
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <chrono>
+#include <random>
+#include <impl.hpp>
+#include <config.h>
+
+template <typename T>
+struct buffer_t {
+   buffer_t(size_t s) { p = new T[s]; }
+   ~buffer_t()        { delete[] p; }
+   buffer_t()                           = default;
+   buffer_t(buffer_t&&)                 = default;
+   buffer_t& operator=(buffer_t&&)      = default;
+   buffer_t(const buffer_t&)            = delete;
+   buffer_t& operator=(const buffer_t&) = delete;
+
+   T* allocate(size_t s)   { return p = new T[s]; }
+   T* operator() ()        { return p; }
+   T& operator[] (size_t i){ return p[i]; }
+private:
+   T* p{nullptr};
+};
+
+struct Mtx {
+
+   template<typename I>
+   static void coo2csc(I *row, I *col, I const* row_coo, I const* col_coo, I nnz, I n, I isOneBased) {
+      // ----- cannot assume that input is already 0!
+      for (I l = 0; l < n+1; l++) col[l] = 0;
+
+      // ----- find the correct column sizes
+      for (I l = 0; l < nnz; l++)
+        col[col_coo[l] - isOneBased]++;
+
+      // ----- cumulative sum
+      for (I i = 0, cumsum = 0; i < n; i++) {
+        I temp = col[i];
+        col[i] = cumsum;
+        cumsum += temp;
+      }
+      col[n] = nnz;
+      // ----- copy the row indices to the correct place
+      for (I l = 0; l < nnz; l++) {
+         I col_l;
+        col_l = col_coo[l] - isOneBased;
+
+        I dst = col[col_l];
+        row[dst] = row_coo[l] - isOneBased;
+
+        col[col_l]++;
+      }
+      // ----- revert the column pointers
+      for (I i = 0, last = 0; i < n; i++) {
+         I temp = col[i];
+        col[i] = last;
+        last = temp;
+      }
+   }
+
+   template<typename I>
+   static bool is_triangular (std::ifstream& file) {
+      std::string line, token;
+      enum state_en {HEADER, SIZE, DATA} state = HEADER;
+      enum LU_t {Z, LOWER, UPPER} LU = Z;
+
+      while (std::getline (file, line, '\n')) {
+         std::stringstream ss(line);
+         switch (state) {
+         case HEADER:
+            ss >> token;
+            if (token != "%%MatrixMarket")   return false;
+            else                             state = SIZE;
+            break;
+         case SIZE:
+            if (line[0] == '%')  continue;
+            else                 state = DATA;
+            break;
+         case DATA:
+            if (line[0] == '%')  continue;
+            I i, j;
+            ss >> i >> j;
+            switch (LU) {
+               case Z:     LU = (i<j) ? UPPER: LOWER; break;
+               case LOWER: if (i<=j)   return false;  break;
+               case UPPER: if (j<=i)   return false;  break;
+            }
+            break;
+         }
+      }
+      file.clear();  // rewind
+      file.seekg(0);
+      return true;
+   }
+
+   template<typename DataT, typename IndexT, MatrixType MatrixT>
+   static bool load (SpMat<DataT, IndexT, MatrixT>& M, std::ifstream& file) {
+      std::string line, token;
+      enum state_en {HEADER, SIZE, DATA} state = HEADER;
+      enum LU_t {Z, LOWER, UPPER} LU = Z;
+      IndexT n1, n2, nnz;
+      buffer_t<IndexT> col{}, row{}, coo_col{}, coo_row{};
+
+      IndexT cnt{};
+      while (std::getline (file, line, '\n')) {
+         std::stringstream ss(line);
+         switch (state) {
+         case HEADER:
+            ss >> token;
+            if (token != "%%MatrixMarket")   return false;
+            else                             state = SIZE;
+            break;
+         case SIZE:
+            if (line[0] == '%')  continue;
+            else {
+               ss >> n1 >> n2 >> nnz;
+               if (session.makeSymmetric)
+                  nnz *= 2;
+               col.allocate(nnz);
+               row.allocate(nnz);
+               coo_col.allocate(nnz);
+               coo_row.allocate(nnz);
+               state = DATA;
+            }
+            break;
+         case DATA:
+            if (line[0] == '%')  continue;
+            IndexT i, j;
+            ss >> i >> j;
+            if (session.makeSymmetric) {
+               if (LU == Z) {
+                  LU = (i<j) ? UPPER: LOWER;
+               }
+               if ((LU==LOWER && j<i) || (LU==UPPER && i<j)) {
+                  coo_row[cnt]   = i;
+                  coo_col[cnt++] = j;
+                  coo_row[cnt]   = j;
+                  coo_col[cnt++] = i;
+               }
+            }
+            else {
+               coo_row[cnt]   = i;
+               coo_col[cnt++] = j;
+            }
+            break;
+         }
+      }
+      coo2csc(&row[0], &col[0], &coo_row[0], &coo_col[0], cnt, n1, 1);
+      M = SpMat<DataT, IndexT, MatrixT>(n1, cnt, &row[0], &col[0]);
+      return true;
+   }
+};
+
+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;
+
+
+   Tpoint start () noexcept { return start_ = std::chrono::steady_clock::now(); }
+   Tpoint stop ()  noexcept { return stop_ = std::chrono::steady_clock::now(); }
+
+   auto dt () noexcept {
+      return std::chrono::duration_cast<std::chrono::microseconds>(stop_ - start_).count();
+   }
+   void print_dt () noexcept {
+      auto t = stop_ - start_;
+      if       (std::chrono::duration_cast<microseconds>(t).count() < 10000)
+         std::cout << "time: " << std::to_string(std::chrono::duration_cast<microseconds>(t).count()) << " [usec]\n";
+      else if (std::chrono::duration_cast<milliseconds>(t).count() < 10000)
+         std::cout << "time: " << std::to_string(std::chrono::duration_cast<milliseconds>(t).count()) << " [msec]\n";
+      else
+         std::cout << "time: " << std::to_string(std::chrono::duration_cast<seconds>(t).count()) << " [sec]\n";
+
+   }
+private:
+   Tpoint start_;
+   Tpoint stop_;
+};
+
+
+void init_ER_graph (matrix& A, double p);
+void print_ER_graph (matrix& A);
+
+#endif /* UTILS_H_ */

inc/v12.h

@@ -0,0 +1,21 @@
+/*!
+ * \file    v12.h
+ * \brief   v1 and v2 part of the exercise header file.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+#ifndef V12_H_
+#define V12_H_
+
+#include <impl.hpp>
+
+namespace v12 {
+
+using matrix = Matrix<int, int>;
+
+int triang_count (matrix& A) ;
+
+};
+#endif /* V12_H_ */

inc/v3.h

@@ -0,0 +1,21 @@
+/*!
+ * \file    v3.h
+ * \brief   v3 part of the exercise header file.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+#ifndef V3_H_
+#define V3_H_
+
+#include <impl.hpp>
+
+namespace v3 {
+
+using matrix = SpMat<int, int>;
+
+int triang_count (matrix& A);
+
+};
+#endif /* V3_H_ */

inc/v4.h

@@ -0,0 +1,21 @@
+/*!
+ * \file    v4.h
+ * \brief   v4 part of the exercise header file.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+#ifndef V4_H_
+#define V4_H_
+
+#include <impl.hpp>
+
+namespace v4 {
+
+using matrix = SpMat<int, int>;
+
+int triang_count (matrix& A);
+
+};
+#endif /* V4_H_ */

src/main.cpp

@@ -0,0 +1,112 @@
+/*!
+ * \file    main.cpp
+ * \brief   Main application file
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+
+#include <iostream>
+#include <string>
+#include <exception>
+
+#include <utils.h>
+#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 == "-i" || arg == "--input") {
+         session.inputMatrix = InputMatrix::MTX;
+         if (i+1 < argc)
+            session.mtxFile = std::ifstream(argv[++i]);
+         else
+            status = false;
+      }
+      else if (arg == "-g" || arg == "--generate")
+         session.inputMatrix = InputMatrix::GENERATE;
+      else if (arg == "-s" || arg == "--size")
+         session.size = (i+1 < argc) ? std::atoi(argv[++i]) : session.size;
+      else if (arg == "-p" || arg == "--probability")
+         session.probability = (i+1 < argc) ? std::atof(argv[++i]) : session.probability;
+      else if (arg == "--print") {
+         session.print = true;
+         session.print_size = (i+1 < argc) ? std::atoi(argv[++i]) : session.print_size;
+      }
+      else if (arg == "--make_symmetric")
+         session.makeSymmetric = true;
+      else if (arg == "-t" || arg == "--timing")
+         session.timing = true;
+      else if (arg == "-h" || arg == "--help") {
+         std::cout << "Help message\n";
+         exit(0);
+      }
+      else {
+         std::cout << "Error message\n";
+         status = false;
+      }
+   }
+   return status;
+}
+
+
+int main(int argc, char* argv[]) try {
+   Timing timer;
+   matrix A;
+
+   // try to read command line
+   if (!get_options(argc, argv))
+      exit(1);
+
+   // get or generate matrix
+   if (session.inputMatrix == InputMatrix::GENERATE) {
+      std::cout << "Initialize matrix with size: " << session.size << " and probability: " << session.probability << '\n';
+      timer.start();
+      A.size(session.size);
+      init_ER_graph(A, session.probability);
+      timer.stop();
+      if (session.timing)  timer.print_dt();
+   }
+   else {
+      std::cout << "Read matrix from file\n";
+      timer.start();
+      if (session.makeSymmetric && !Mtx::is_triangular<matrix::indexType> (session.mtxFile))
+         throw std::runtime_error("Error: Matrix is not strictly upper or lower");
+      if (!Mtx::load (A, session.mtxFile)) {
+         throw std::runtime_error("Error: fail to load matrix");
+      }
+      timer.stop();
+      std::cout << "Matrix size: " << A.size() << " and capacity: " << A.capacity() <<'\n';
+      if (session.timing)  timer.print_dt();
+   }
+
+   if (session.print) {
+      std::cout << "Array A:\n";
+      print_ER_graph (A);
+   }
+
+   std::cout << "count triangles\n";
+   timer.start();
+   std::cout << "There are " << triang_count(A) << " triangles\n";
+   timer.stop();
+   if (session.timing)  timer.print_dt();
+
+	return 0;
+}
+catch (std::exception& e) {
+    //we probably pollute the user's screen. Comment `cerr << ...` if you don't like it.
+    std::cerr << e.what() << '\n';
+   exit(1);
+}

src/utils.cpp

@@ -0,0 +1,53 @@
+/*!
+ * \file    utils.cpp
+ * \brief   Utilities to handle matrix files, chrono, etc...
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+
+#include <utils.h>
+#include <algorithm>
+
+/*!
+ * Initialize the matrix as Erdős-Rényi graph
+ * \param A    The matrix to initialize
+ * \param p    The probability of each edge
+ */
+void init_ER_graph (matrix& A, double p) {
+   std::random_device rd;
+   std::mt19937 gen(rd());
+   std::binomial_distribution<> d(1, p);
+
+#if CODE_VERSION == V12
+   std::transform (A.begin(), A.end(), A.begin(),
+         [&] (int x) { std::ignore =x; return d(gen); }
+#else
+   A.for_each_in(0, A.size(), [&](auto i) {
+      A.for_each_in(i+1, A.size(), [&](auto j){
+         matrix::dataType edge = d(gen);
+         if (edge) {
+            A.set(edge, i, j);
+            A.set(edge, j, i);
+         }
+      });
+   });
+#endif
+
+}
+
+/*!
+ * Utility to print the graph to sdtout
+ */
+void print_ER_graph (matrix& A) {
+   matrix::indexType N = (A.size() < (matrix::indexType)session.print_size) ? A.size() : session.print_size;
+
+   A.for_each_in(0, N, [&](auto i){
+      A.for_each_in(0, N, [&](auto j) {
+         std::cout << A(i, j) << ' ';
+      });
+      std::cout << '\n';
+   });
+}
+

src/v12.cpp

@@ -0,0 +1,35 @@
+/*!
+ * \file    v12.cpp
+ * \brief   v1 and v2 part of the exercise.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+
+#include <iostream>
+#include <random>
+#include <v12.h>
+
+namespace v12 {
+
+/*!
+ * A naive triangle counting algorithm
+ * \param A    The adjacency matrix
+ * \return     The number of triangles
+ */
+int triang_count (matrix& A) {
+   int count =0;
+
+   // We use a symmetric matrix so we iterate using the constrain i<j<k
+   A.for_each_in(0, A.size(), [&](auto i) {
+      A.for_each_in(i+1, A.size(), [&](auto j) {
+         A.for_each_in(j+1, A.size(), [&](auto k){
+            count += (A(i,j) && A(i,k) && A(j,k)) ? 1:0;
+         });
+      });
+   });
+   return count;
+}
+
+}

src/v3.cpp

@@ -0,0 +1,48 @@
+/*!
+ * \file    v3.cpp
+ * \brief   vv3 part of the exercise.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+
+#include <iostream>
+#include <random>
+#include <v3.h>
+
+namespace v3 {
+
+using index_t  = typename matrix::indexType;
+using value_t  = typename matrix::dataType;
+
+/*!
+ * A naive triangle counting algorithm
+ * \param A    The adjacency matrix
+ * \return     The number of triangles
+ */
+std::vector<value_t> triang_v(matrix& A) {
+   std::vector<value_t> c(A.size());
+
+   for (int i=0 ; i<A.size() ; ++i) {
+      for (auto j = A.getCol(i); j.index() != j.end() ; ++j)               // j list all the edges with i
+         for (auto k = A.getCol(j.index()); k.index() != k.end() ; ++k)    // k list all the edges with j
+            for (auto ii = A.getCol(i) ; ii.index() <= k.index() ; ++ii)   // search for i-k edge (this could be binary search)
+               if (ii.index() == k.index())  ++c[i];
+   }
+   return c;
+}
+
+value_t sum (std::vector<value_t>& v) {
+   value_t s =0;
+   for (auto& it : v)
+      s += it;
+   return s;
+}
+
+value_t triang_count (matrix& A) {
+   auto v = triang_v(A);
+   return sum(v);
+}
+
+}

src/v4.cpp

@@ -0,0 +1,43 @@
+/*!
+ * \file    v4.cpp
+ * \brief   vv3 part of the exercise.
+ *
+ * \author
+ *    Christos Choutouridis AEM:8997
+ *    <cchoutou@ece.auth.gr>
+ */
+
+#include <iostream>
+#include <random>
+#include <v4.h>
+
+namespace v4 {
+
+using index_t  = typename matrix::indexType;
+using value_t  = typename matrix::dataType;
+
+
+std::vector<value_t> mmacc_v(matrix& A, matrix& B) {
+   std::vector<value_t> c(A.size());
+   for (int i=0 ; i<A.size() ; ++i) {
+      for (auto j = A.getRow(i); j.index() != j.end() ; ++j){
+         c[i] += A.getRow(i)*B.getCol(j.index());
+      }
+   }
+   return c;
+}
+
+
+value_t sum (std::vector<value_t>& v) {
+   value_t s =0;
+   for (auto& it : v)
+      s += it;
+   return s;
+}
+
+value_t triang_count (matrix& A) {
+   auto v = mmacc_v(A, A);
+   return (session.makeSymmetric) ? sum(v)/6 : sum(v);
+}
+
+}