Refactor progress

Next step is ramdom and equal. I also need to have a look at add and the serial naming of the functions.

include/numerics/detail/abs_serial.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <cmath> // std::abs
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_abs_scalar_serial(T& c) {
+		c = static_cast<T>(std::abs(c));
+	}
+		
+	template <typename T>
+	void inplace_abs_elementwise_serial(utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = static_cast<T>(std::abs(A(i,j)));
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_abs_elementwise_serial(utils::Vector<T>& v) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = static_cast<T>(std::abs(v[i]));
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/add_serial.h

@@ -0,0 +1,85 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Scalar ----------------
+	template <typename T>
+	void inplace_add_scalar_serial(utils::Matrix<T>& A, const T c) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) += c;
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_add_scalar_serial(utils::Vector<T>& v, const T c) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] += c;
+		}
+	}
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_add_elementwise_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_add_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) += B(i,j);
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_add_elementwise_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    if (v.size() != p.size()) {
+	        throw std::runtime_error("inplace_add_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] += p[i];
+		}
+	}
+
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	void inplace_add_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (cols != v.size()) {
+	        throw std::runtime_error("inplace_add_rowwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) += v[j];
+			}
+		}
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	void inplace_add_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != v.size()) {
+	        throw std::runtime_error("inplace_add_colwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+		    const T vi = v[i];
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) += vi;
+			}
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/argmax_serial.h

@@ -0,0 +1,83 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+//#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	utils::Vector<uint64_t> argmax_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    utils::Vector<uint64_t> idx(2, 0);
+		T max = A(0,0);
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (max < A(i,j)){
+					max = A(i,j);
+					idx[0] = i;
+					idx[1] = j;
+				}
+			}
+		}
+		return idx;
+	}
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	uint64_t argmax_serial(const utils::Vector<T>& v) {
+		const uint64_t N = v.size();
+		uint64_t idx = 0;
+		T max = v[0];
+		for (uint64_t i = 1; i < N; ++i){
+			if (max < v[i]){
+				max = v[i];
+				idx = i;
+			}
+		}
+		return idx;
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	utils::Vector<uint64_t> argmax_rowwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    utils::Vector<T> max(rows, T{0});
+		utils::Vector<uint64_t> idx(rows, 0);
+		for (uint64_t i = 0; i < rows; ++i){
+			max[i] = A(i,0);
+			for (uint64_t j = 1; j < cols; ++j){
+				if (max[i] < A(i,j)){
+					max[i] = A(i,j);
+					idx[i] = j;
+
+				}
+			}
+		}
+		return idx;
+	}
+
+	template <typename T>
+	utils::Vector<uint64_t> argmax_colwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> max(cols, T{0});
+		utils::Vector<uint64_t> idx(cols, 0);
+		for (uint64_t j = 0; j < cols; ++j){
+			max[j] = A(0, j);
+			for (uint64_t i = 1; i < rows; ++i){
+				if (max[j] < A(i,j)){
+					max[j] = A(i,j);
+					idx[j] = i;
+				}
+			}
+		}
+		return idx;
+	}
+} // namespace numerics
+

include/numerics/detail/argmin_serial.h

@@ -0,0 +1,83 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+//#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	utils::Vector<uint64_t> argmin_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    utils::Vector<uint64_t> idx(2, 0);
+		T min = A(0,0);
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (min > A(i,j)){
+					min = A(i,j);
+					idx[0] = i;
+					idx[1] = j;
+				}
+			}
+		}
+		return idx;
+	}
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	uint64_t argmin_serial(const utils::Vector<T>& v) {
+		const uint64_t N = v.size();
+		uint64_t idx = 0;
+		T min = v[0];
+		for (uint64_t i = 1; i < N; ++i){
+			if (min > v[i]){
+				min = v[i];
+				idx = i;
+			}
+		}
+		return idx;
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	utils::Vector<uint64_t> argmin_rowwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    utils::Vector<T> min(rows, T{0});
+		utils::Vector<uint64_t> idx(rows, 0);
+		for (uint64_t i = 0; i < rows; ++i){
+			min[i] = A(i,0);
+			for (uint64_t j = 1; j < cols; ++j){
+				if (min[i] > A(i,j)){
+					min[i] = A(i,j);
+					idx[i] = j;
+
+				}
+			}
+		}
+		return idx;
+	}
+
+	template <typename T>
+	utils::Vector<uint64_t> argmin_colwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> min(cols, T{0});
+		utils::Vector<uint64_t> idx(cols, 0);
+		for (uint64_t j = 0; j < cols; ++j){
+			min[j] = A(0, j);
+			for (uint64_t i = 1; i < rows; ++i){
+				if (min[j] > A(i,j)){
+					min[j] = A(i,j);
+					idx[j] = i;
+				}
+			}
+		}
+		return idx;
+	}
+} // namespace numerics
+

include/numerics/detail/clip_serial.h

@@ -0,0 +1,136 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ------------- Helper Functions -----------
+	template <typename T>
+	inline T clip_low_value(T x, const T low) {
+	    if (x < low) {
+	        return low;
+	    }
+	    return x;
+	}
+
+	template <typename T>
+	inline T clip_high_value(T x, const  T high) {
+	    if (x > high) {
+	        return high;
+	    }
+	    return x;
+	}
+
+	template <typename T>
+	inline T clip_value(T x, const  T low, const  T high) {
+	    if (x < low) {
+	        return low;
+	    }
+	    if (x > high) {
+	        return high;
+	    }
+	    return x;
+	}
+
+	//---------------------------------------------
+	// -----------------	Clip 	---------------
+	//---------------------------------------------
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	inline void inplace_clip_scalar_serial(T& c, const T low, const T high) {
+
+	    if (low > high) {
+	        throw std::runtime_error("inplace_clip_scalar_serial: lower clip  > higher clip");
+	    }
+	    c = detail::clip_value(c, low, high);
+
+	}
+		
+	template <typename T>
+	inline void inplace_clip_elementwise_serial(utils::Matrix<T>& A, const T low, const T high) {
+	    if (low > high) {
+	        throw std::runtime_error("inplace_clip_elementwise_serial: lower clip  > higher clip");
+	    }
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = detail::clip_value(A(i,j), low, high);
+			}
+		}
+	}
+
+	template <typename T>
+	inline void inplace_clip_elementwise_serial(utils::Vector<T>& v, const T low, const T high) {
+	    if (low > high) {
+	        throw std::runtime_error("inplace_clip_elementwise_serial: lower clip  > higher clip");
+	    }
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = detail::clip_value(v[i], low, high);
+		}
+	}
+
+	//---------------------------------------------
+	// -----------------  Clip Low  ---------------
+	//---------------------------------------------
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	inline void inplace_clip_low_scalar_serial(T& c, const T low) {
+	    c = detail::clip_low_value(c, low);
+
+	}
+		
+	template <typename T>
+	inline void inplace_clip_low_elementwise_serial(utils::Matrix<T>& A, const T low) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = detail::clip_low_value(A(i,j), low);
+			}
+		}
+	}
+
+	template <typename T>
+	inline void inplace_clip_low_elementwise_serial(utils::Vector<T>& v, const T low) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = detail::clip_low_value(v[i], low);
+		}
+	}
+
+	//---------------------------------------------
+	// ----------------- Clip High  ---------------
+	//---------------------------------------------
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	inline void inplace_clip_high_scalar_serial(T& c, const T high) {
+	    c = detail::clip_high_value(c, high);
+
+	}
+		
+	template <typename T>
+	inline void inplace_clip_high_elementwise_serial(utils::Matrix<T>& A, const T high) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = detail::clip_high_value(A(i,j), high);
+			}
+		}
+	}
+
+	template <typename T>
+	inline void inplace_clip_high_elementwise_serial(utils::Vector<T>& v, const T high) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = detail::clip_high_value(v[i], high);
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/div_serial.h

@@ -0,0 +1,85 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Scalar ----------------
+	template <typename T>
+	void inplace_div_scalar_serial(utils::Matrix<T>& A, const T c) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) /= c;
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_div_scalar_serial(utils::Vector<T>& v, const T c) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] /= c;
+		}
+	}
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_div_elementwise_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_div_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) /= B(i,j);
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_div_elementwise_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    if (v.size() != p.size()) {
+	        throw std::runtime_error("inplace_div_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] /= p[i];
+		}
+	}
+
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	void inplace_div_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (cols != v.size()) {
+	        throw std::runtime_error("inplace_div_rowwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) /= v[j];
+			}
+		}
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	void inplace_div_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != v.size()) {
+	        throw std::runtime_error("inplace_div_colwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+		    const T vi = v[i];
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) /= vi;
+			}
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/exp_serial.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <cmath> // std::exp
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_exp_scalar_serial(T& c) {
+		c = static_cast<T>(std::exp(c));
+	}
+		
+	template <typename T>
+	void inplace_exp_elementwise_serial(utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = static_cast<T>(std::exp(A(i,j)));
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_exp_elementwise_serial(utils::Vector<T>& v) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = static_cast<T>(std::exp(v[i]));
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/log_serial.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <cmath> // std::log
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_log_scalar_serial(T& c) {
+		c = static_cast<T>(std::log(c));
+	}
+		
+	template <typename T>
+	void inplace_log_elementwise_serial(utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = static_cast<T>(std::log(A(i,j)));
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_log_elementwise_serial(utils::Vector<T>& v) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = static_cast<T>(std::log(v[i]));
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/max_serial.h

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+//#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	T max_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		T max = A(0,0);
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (max < A(i,j)){
+					max = A(i,j);
+				}
+			}
+		}
+		return max;
+	}
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	T max_serial(const utils::Vector<T>& v) {
+		const uint64_t N = v.size();
+		T max = v[0];
+		for (uint64_t i = 0; i < N; ++i){
+			if (max < v[i]){
+				max = v[i];
+			}
+		}
+		return max;
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	utils::Vector<T> max_rowwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> max(rows, T{0});
+		for (uint64_t i = 0; i < rows; ++i){
+			max[i] = A(i,0);
+			for (uint64_t j = 1; j < cols; ++j){
+				if (max[i] < A(i,j)){
+					max[i] = A(i,j);
+				}
+			}
+		}
+		return max;
+	}
+
+	template <typename T>
+	utils::Vector<T> max_colwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> max(cols, T{0});
+		for (uint64_t j = 0; j < cols; ++j){
+			max[j] = A(0, j);
+			for (uint64_t i = 1; i < rows; ++i){
+				if (max[j] < A(i,j)){
+					max[j] = A(i,j);
+				}
+			}
+		}
+		return max;
+	}
+} // namespace numerics
+

include/numerics/detail/mean_serial.h

@@ -0,0 +1,55 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+//#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+#include "sum_serial.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	T mean_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		T sum = detail::sum_serial(A);
+		sum /= static_cast<T>(rows*cols);
+		return sum;
+	}
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	T mean_serial(const utils::Vector<T>& v) {
+		const uint64_t N = v.size();
+		T sum = detail::sum_serial(v);
+		sum /= static_cast<T>(N);
+		return sum;
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	utils::Vector<T> mean_rowwise_serial(const utils::Matrix<T>& A) {
+		const T cols = static_cast<T>(A.cols());
+		utils::Vector<T> sum = detail::sum_rowwise_serial(A);
+		const uint64_t N = sum.size();
+		for (uint64_t i = 0; i < N; ++i){
+			sum[i] /= cols;
+		}
+		return sum;
+	}
+
+	template <typename T>
+	utils::Vector<T> mean_colwise_serial(const utils::Matrix<T>& A) {
+		const T rows = static_cast<T>(A.rows());
+		utils::Vector<T> sum = detail::sum_colwise_serial(A);
+		const uint64_t N = sum.size();
+		for (uint64_t i = 0; i < N; ++i){
+			sum[i] /= rows;
+		}
+		return sum;
+	}
+} // namespace numerics
+

include/numerics/detail/min_serial.h

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+//#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	T min_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		T min = A(0,0);
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (min > A(i,j)){
+					min = A(i,j);
+				}
+			}
+		}
+		return min;
+	}
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	T min_serial(const utils::Vector<T>& v) {
+		const uint64_t N = v.size();
+		T min = v[0];
+		for (uint64_t i = 0; i < N; ++i){
+			if (min > v[i]){
+				min = v[i];
+			}
+		}
+		return min;
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	utils::Vector<T> min_rowwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> min(rows, T{0});
+		for (uint64_t i = 0; i < rows; ++i){
+			min[i] = A(i,0);
+			for (uint64_t j = 1; j < cols; ++j){
+				if (min[i] > A(i,j)){
+					min[i] = A(i,j);
+				}
+			}
+		}
+		return min;
+	}
+
+	template <typename T>
+	utils::Vector<T> min_colwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> min(cols, T{0});
+		for (uint64_t j = 0; j < cols; ++j){
+			min[j] = A(0, j);
+			for (uint64_t i = 1; i < rows; ++i){
+				if (min[j] > A(i,j)){
+					min[j] = A(i,j);
+				}
+			}
+		}
+		return min;
+	}
+} // namespace numerics
+

include/numerics/detail/mul_serial.h

@@ -0,0 +1,85 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Scalar ----------------
+	template <typename T>
+	void inplace_mul_scalar_serial(utils::Matrix<T>& A, const T c) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) *= c;
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_mul_scalar_serial(utils::Vector<T>& v, const T c) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] *= c;
+		}
+	}
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_mul_elementwise_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_mul_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) *= B(i,j);
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_mul_elementwise_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    if (v.size() != p.size()) {
+	        throw std::runtime_error("inplace_mul_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] *= p[i];
+		}
+	}
+
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	void inplace_mul_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (cols != v.size()) {
+	        throw std::runtime_error("inplace_mul_rowwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) *= v[j];
+			}
+		}
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	void inplace_mul_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != v.size()) {
+	        throw std::runtime_error("inplace_mul_colwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+		    const T vi = v[i];
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) *= vi;
+			}
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/neg_serial.h

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_neg_scalar_serial(T& c) {
+		c = -c;
+	}
+		
+	template <typename T>
+	void inplace_neg_elementwise_serial(utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = -A(i,j);
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_neg_elementwise_serial(utils::Vector<T>& v) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = v[i];
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/pow_serial.h

@@ -0,0 +1,55 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <stdexcept> // std::runtime_error
+#include <cmath> // std::pow
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Scalar ----------------
+	template <typename T>
+	void inplace_pow_scalar_serial(utils::Matrix<T>& A, const T c) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = static_cast<T>(std::pow(A(i,j),c));
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_pow_scalar_serial(utils::Vector<T>& v, const T c) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = static_cast<T>(std::pow(v[i], c));
+		}
+	}
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_pow_elementwise_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_pow_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = static_cast<T>(std::pow(A(i,j), B(i,j)));
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_pow_elementwise_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    if (v.size() != p.size()) {
+	        throw std::runtime_error("inplace_pow_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = static_cast<T>(std::pow(v[i], p[i]));
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/sqrt_serial.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <cmath> // std::sqrt
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_sqrt_scalar_serial(T& c) {
+		c = static_cast<T>(std::sqrt(c));
+	}
+		
+	template <typename T>
+	void inplace_sqrt_elementwise_serial(utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) = static_cast<T>(std::sqrt(A(i,j)));
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_sqrt_elementwise_serial(utils::Vector<T>& v) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] = static_cast<T>(std::sqrt(v[i]));
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/sub_serial.h

@@ -0,0 +1,85 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Scalar ----------------
+	template <typename T>
+	void inplace_sub_scalar_serial(utils::Matrix<T>& A, const T c) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) -= c;
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_add_scalar_serial(utils::Vector<T>& v, const T c) {
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] -= c;
+		}
+	}
+	// ---------------- Elemenwise ----------------
+	template <typename T>
+	void inplace_sub_elementwise_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_sub_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) -= B(i,j);
+			}
+		}
+	}
+
+	template <typename T>
+	void inplace_sub_elementwise_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    if (v.size() != p.size()) {
+	        throw std::runtime_error("inplace_sub_elementwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < v.size(); ++i){
+			v[i] -= p[i];
+		}
+	}
+
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	void inplace_sub_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (cols != v.size()) {
+	        throw std::runtime_error("inplace_sub_rowwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) -= v[j];
+			}
+		}
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	void inplace_sub_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    if (rows != v.size()) {
+	        throw std::runtime_error("inplace_sub_colwise_serial: dimension mismatch");
+	    }
+		for (uint64_t i = 0; i < rows; ++i){
+		    const T vi = v[i];
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) -= vi;
+			}
+		}
+	}
+} // namespace numerics
+

include/numerics/detail/sum_serial.h

@@ -0,0 +1,63 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+//#include <stdexcept> // std::runtime_error
+
+#include "../utils/vector.h"
+#include "../utils/matrix.h"
+
+namespace numerics::detail{
+	
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	T sum_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		T sum = T{0};
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				sum += A(i,j);
+			}
+		}
+		return sum;
+	}
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	T sum_serial(const utils::Vector<T>& v) {
+		const uint64_t N = v.size();
+		T sum = T{0};
+		for (uint64_t i = 0; i < N; ++i){
+			sum += v[i];
+		}
+		return sum;
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	utils::Vector<T> sum_rowwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> sum(cols, T{0});
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				sum[j] += A(i,j);
+			}
+		}
+		return sum;
+	}
+
+	template <typename T>
+	utils::Vector<T> sum_colwise_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+		utils::Vector<T> sum(rows, T{0});
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				sum[i] += A(i,j);
+			}
+		}
+		return sum;
+	}
+} // namespace numerics
+