diff --git a/include/numerics/abs.h b/include/numerics/abs.h
index f2af0a8..ebc8d5c 100644
--- a/include/numerics/abs.h
+++ b/include/numerics/abs.h
@@ -1,23 +1,46 @@
 #pragma once
 
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+#include "detail/abs_serial.h"
 
 
 namespace numerics{
-	
-	template <typename T>
-	T abs(const T a){
 
-		if(a < 0){
-			return -a;
-		}else{
-			return a;
-		}
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_abs(T& c) {
+		detail::inplace_abs_scalar_serial(c);
+	}
+	template <typename T>
+	inline T abs(const T c) {
+		T out = c;
+		inplace_abs(out);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_abs(utils::Matrix<T>& A) {
+		detail::inplace_abs_elementwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> abs(const utils::Matrix<T>& A) {
+	    utils::Matrix<T> out = A;
+	    inplace_abs(out);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_abs(utils::Vector<T>& v) {
+		detail::inplace_abs_elementwise_serial(v);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> abs(const utils::Vector<T>& v) {
+	    utils::Vector<T> out = v;
+	    inplace_abs(out);
+	    return out;
 	}
 
 
-
-} // namespace numerics
-
+}
\ No newline at end of file
diff --git a/include/numerics/add.h b/include/numerics/add.h
new file mode 100644
index 0000000..fda479f
--- /dev/null
+++ b/include/numerics/add.h
@@ -0,0 +1,83 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/add_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Scalar ----------------
+	template <typename T>
+	inline void inplace_add(utils::Matrix<T>& A, const T b) {
+		detail::inplace_add_scalar_serial(A,b);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> add(const utils::Matrix<T>& A, const T b) {
+	    utils::Matrix<T> out = A;
+	    inplace_add(out, b);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_add(utils::Vector<T>& v, const T b) {
+		detail::inplace_add_scalar_serial(v,b);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> add(const utils::Vector<T>& v, const T b) {
+	    utils::Vector<T> out = v;
+	    inplace_add(out, b);
+	    return out;
+	}
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_add(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		detail::inplace_add_elementwise_serial(A,B);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> add(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    utils::Matrix<T> out = A;
+	    inplace_add(out, B);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_add(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_add_elementwise_serial(v,p);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> add(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    utils::Vector<T> out = v;
+	    inplace_add(out, p);
+	    return out;
+	}
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	inline void inplace_add_rowwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_add_rowwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> add_rowwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_add_rowwise(out, v);
+	    return out;
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	inline void inplace_add_colwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_add_colwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> add_colwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_add_colwise(out, v);
+	    return out;
+	}
+
+}
\ No newline at end of file
diff --git a/include/numerics/argmax.h b/include/numerics/argmax.h
new file mode 100644
index 0000000..785917d
--- /dev/null
+++ b/include/numerics/argmax.h
@@ -0,0 +1,33 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/argmax_serial.h"
+
+
+
+namespace numerics{
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	inline uint64_t argmax(const utils::Vector<T>& v) {
+		return detail::argmax_serial(v);
+	}
+
+
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	inline utils::Vector<uint64_t> argmax(const utils::Matrix<T>& A) {
+		return detail::argmax_serial(A);
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	inline utils::Vector<uint64_t> argmax_rowwise(const utils::Matrix<T>& A) {
+		return detail::argmax_rowwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Vector<uint64_t> argmax_colwise(const utils::Matrix<T>& A) {
+		return detail::argmax_colwise_serial(A);
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/argmin.h b/include/numerics/argmin.h
new file mode 100644
index 0000000..52f9e3b
--- /dev/null
+++ b/include/numerics/argmin.h
@@ -0,0 +1,33 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/argmin_serial.h"
+
+
+
+namespace numerics{
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	inline uint64_t argmin(const utils::Vector<T>& v) {
+		return detail::argmin_serial(v);
+	}
+
+
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	inline utils::Vector<uint64_t> argmin(const utils::Matrix<T>& A) {
+		return detail::argmin_serial(A);
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	inline utils::Vector<uint64_t> argmin_rowwise(const utils::Matrix<T>& A) {
+		return detail::argmin_rowwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Vector<uint64_t> argmin_colwise(const utils::Matrix<T>& A) {
+		return detail::argmin_colwise_serial(A);
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/clip.h b/include/numerics/clip.h
new file mode 100644
index 0000000..2de4640
--- /dev/null
+++ b/include/numerics/clip.h
@@ -0,0 +1,128 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/clip_serial.h"
+
+
+namespace numerics{
+
+	//---------------------------------------------
+	// -----------------	Clip 	---------------
+	//---------------------------------------------
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_clip(T& c, const T low, const T high) {
+		detail::inplace_clip_scalar_serial(c, low, high);
+	}
+	template <typename T>
+	inline T clip(const T c, const  T low, const  T high) {
+		T out = c;
+		inplace_clip(out, low, high);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_clip(utils::Matrix<T>& A, const T low, const T high) {
+		detail::inplace_clip_elementwise_serial(A, low, high);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> clip(const utils::Matrix<T>& A, const T low, const T high) {
+	    utils::Matrix<T> out = A;
+	    inplace_clip(out, low, high);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_clip(utils::Vector<T>& v, const T low, const T high) {
+		detail::inplace_clip_elementwise_serial(v, low, high);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> clip(const utils::Vector<T>& v, const T low, const T high) {
+	    utils::Vector<T> out = v;
+	    inplace_clip(out, low, high);
+	    return out;
+	}
+
+	//---------------------------------------------
+	// -----------------  Clip Low  ---------------
+	//---------------------------------------------
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_clip_low(T& c, const T low) {
+		detail::inplace_clip_low_scalar_serial(c, low);
+	}
+	template <typename T>
+	inline T clip_low(const T c, const  T low) {
+		T out = c;
+		inplace_clip_low(out, low);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_clip_low(utils::Matrix<T>& A, const T low) {
+		detail::inplace_clip_low_elementwise_serial(A, low);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> clip_low(const utils::Matrix<T>& A, const T low) {
+	    utils::Matrix<T> out = A;
+	    inplace_clip_low(out, low);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_clip_low(utils::Vector<T>& v, const T low) {
+		detail::inplace_clip_low_elementwise_serial(v, low);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> clip_low(const utils::Vector<T>& v, const T low) {
+	    utils::Vector<T> out = v;
+	    inplace_clip_low(out, low);
+	    return out;
+	}
+
+	//---------------------------------------------
+	// -----------------  Clip High ---------------
+	//---------------------------------------------
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_clip_high(T& c, const T high) {
+		detail::inplace_clip_high_scalar_serial(c, high);
+	}
+	template <typename T>
+	inline T clip_high(const T c, const  T high) {
+		T out = c;
+		inplace_clip_high(out, high);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_clip_high(utils::Matrix<T>& A, const T high) {
+		detail::inplace_clip_high_elementwise_serial(A, high);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> clip_high(const utils::Matrix<T>& A, const T high) {
+	    utils::Matrix<T> out = A;
+	    inplace_clip_high(out, high);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_clip_high(utils::Vector<T>& v, const T high) {
+		detail::inplace_clip_high_elementwise_serial(v, high);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> clip_high(const utils::Vector<T>& v, const T high) {
+	    utils::Vector<T> out = v;
+	    inplace_clip_high(out, high);
+	    return out;
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/detail/abs_serial.h b/include/numerics/detail/abs_serial.h
new file mode 100644
index 0000000..f16e6e5
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/add_serial.h b/include/numerics/detail/add_serial.h
new file mode 100644
index 0000000..1f1dc86
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/argmax_serial.h b/include/numerics/detail/argmax_serial.h
new file mode 100644
index 0000000..6e18b39
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/argmin_serial.h b/include/numerics/detail/argmin_serial.h
new file mode 100644
index 0000000..42ee54c
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/clip_serial.h b/include/numerics/detail/clip_serial.h
new file mode 100644
index 0000000..326f521
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/div_serial.h b/include/numerics/detail/div_serial.h
new file mode 100644
index 0000000..48f8a9d
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/exp_serial.h b/include/numerics/detail/exp_serial.h
new file mode 100644
index 0000000..2111ccc
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/log_serial.h b/include/numerics/detail/log_serial.h
new file mode 100644
index 0000000..cead6e4
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/max_serial.h b/include/numerics/detail/max_serial.h
new file mode 100644
index 0000000..9f8b357
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/mean_serial.h b/include/numerics/detail/mean_serial.h
new file mode 100644
index 0000000..33d3fda
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/min_serial.h b/include/numerics/detail/min_serial.h
new file mode 100644
index 0000000..3a99710
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/mul_serial.h b/include/numerics/detail/mul_serial.h
new file mode 100644
index 0000000..063f834
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/neg_serial.h b/include/numerics/detail/neg_serial.h
new file mode 100644
index 0000000..365f69c
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/pow_serial.h b/include/numerics/detail/pow_serial.h
new file mode 100644
index 0000000..c54a753
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/sqrt_serial.h b/include/numerics/detail/sqrt_serial.h
new file mode 100644
index 0000000..2e24a63
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/sub_serial.h b/include/numerics/detail/sub_serial.h
new file mode 100644
index 0000000..b480257
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/detail/sum_serial.h b/include/numerics/detail/sum_serial.h
new file mode 100644
index 0000000..3531c50
--- /dev/null
+++ b/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
+
diff --git a/include/numerics/div.h b/include/numerics/div.h
new file mode 100644
index 0000000..042d655
--- /dev/null
+++ b/include/numerics/div.h
@@ -0,0 +1,83 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/div_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Scalar ----------------
+	template <typename T>
+	inline void inplace_div(utils::Matrix<T>& A, const T b) {
+		detail::inplace_div_scalar_serial(A,b);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> div(const utils::Matrix<T>& A, const T b) {
+	    utils::Matrix<T> out = A;
+	    inplace_div(out, b);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_div(utils::Vector<T>& v, const T b) {
+		detail::inplace_div_scalar_serial(v,b);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> div(const utils::Vector<T>& v, const T b) {
+	    utils::Vector<T> out = v;
+	    inplace_div(out, b);
+	    return out;
+	}
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_div(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		detail::inplace_div_elementwise_serial(A,B);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> div(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    utils::Matrix<T> out = A;
+	    inplace_div(out, B);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_div(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_div_elementwise_serial(v,p);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> div(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    utils::Vector<T> out = v;
+	    inplace_div(out, p);
+	    return out;
+	}
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	inline void inplace_div_rowwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_div_rowwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> div_rowwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_div_rowwise(out, v);
+	    return out;
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	inline void inplace_div_colwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_div_colwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> dub_colwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_div_colwise(out, v);
+	    return out;
+	}
+
+}
\ No newline at end of file
diff --git a/include/numerics/exp.h b/include/numerics/exp.h
index 238e988..ca39e95 100644
--- a/include/numerics/exp.h
+++ b/include/numerics/exp.h
@@ -1,18 +1,46 @@
 #pragma once
 
-#include <cmath>
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+#include "detail/exp_serial.h"
 
 
 namespace numerics{
-	
+
+	// ---------------- Elementwise ----------------
 	template <typename T>
-	T exp(const T a){
-		return std::exp(a);
+	inline void inplace_exp(T& c) {
+		detail::inplace_exp_scalar_serial(c);
+	}
+	template <typename T>
+	inline T exp(const T c) {
+		T out = c;
+		inplace_exp(out);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_exp(utils::Matrix<T>& A) {
+		detail::inplace_exp_elementwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> exp(const utils::Matrix<T>& A) {
+	    utils::Matrix<T> out = A;
+	    inplace_exp(out);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_exp(utils::Vector<T>& v) {
+		detail::inplace_exp_elementwise_serial(v);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> exp(const utils::Vector<T>& v) {
+	    utils::Vector<T> out = v;
+	    inplace_exp(out);
+	    return out;
 	}
 
 
-} // namespace numerics
-
+}
\ No newline at end of file
diff --git a/include/numerics/log.h b/include/numerics/log.h
index d178580..6b36feb 100644
--- a/include/numerics/log.h
+++ b/include/numerics/log.h
@@ -1,25 +1,46 @@
 #pragma once
 
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+#include "detail/log_serial.h"
 
 
 namespace numerics{
-	
 
+	// ---------------- Elementwise ----------------
 	template <typename T>
-	void inplace_log(T a){
-		a = std::log(a);
+	inline void inplace_log(T& c) {
+		detail::inplace_log_scalar_serial(c);
+	}
+	template <typename T>
+	inline T log(const T c) {
+		T out = c;
+		inplace_log(out);
+		return out;
 	}
 
 	template <typename T>
-	T log(const T a){
-		T b = a;
-		inplace_log(b);
-		return b;
+	inline void inplace_log(utils::Matrix<T>& A) {
+		detail::inplace_log_elementwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> log(const utils::Matrix<T>& A) {
+	    utils::Matrix<T> out = A;
+	    inplace_log(out);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_log(utils::Vector<T>& v) {
+		detail::inplace_log_elementwise_serial(v);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> log(const utils::Vector<T>& v) {
+	    utils::Vector<T> out = v;
+	    inplace_log(out);
+	    return out;
 	}
 
 
-} // namespace numerics
-
+}
\ No newline at end of file
diff --git a/include/numerics/matabs.h b/include/numerics/matabs.h
deleted file mode 100644
index b9b73b5..0000000
--- a/include/numerics/matabs.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#pragma once
-
-#include "./numerics/abs.h"
-#include "./utils/matrix.h"
-
-namespace numerics{
-	
-	template <typename T>
-	void inplace_matabs(utils::Matrix<T>& A){
-	
-		for (uint64_t i = 0; i < A.rows(); ++i){
-			for (uint64_t j = 0; j < A.cols(); ++j){
-				A(i,j) = numerics::abs(A(i,j));
-			}
-		}
-	}
-
-
-	template <typename T>
-	utils::Matrix<T> matabs(const utils::Matrix<T>& A){
-		utils::Matrix<T> B = A;
-		inplace_matabs(B);
-		return B;
-	}
-
-
-
-} // namespace numerics
-
diff --git a/include/numerics/matadd.h b/include/numerics/matadd.h
deleted file mode 100644
index 624e3c3..0000000
--- a/include/numerics/matadd.h
+++ /dev/null
@@ -1,132 +0,0 @@
-#ifndef _matadd_n_
-#define _matadd_n_
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-#include "./core/omp_config.h"
-
-namespace numerics{
-	
-	template <typename T>
-	void inplace_matadd_mat(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_matadd: dimension mismatch");
-	    }
-
-	    for (uint64_t i = 0; i < cols; ++i) {
-	        for (uint64_t j = 0; j < rows; ++j) {
-	            A(j, i) +=  B(j, i);
-	        }
-	    }
-	}
-
-	template <typename T>
-	utils::Matrix<T> matadd_mat(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
-	    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-	    utils::Matrix<T> C = A;
-
-	    if (rows != B.rows() || cols != B.cols()) {
-	        throw std::runtime_error("inplace_matadd: dimension mismatch");
-	    }
-
-	    inplace_matadd_mat(C, B);
-	    return C;
-	}
-
-	template <typename T>
-	void inplace_matadd_colvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-	    if (rows != x.size()) {
-	        throw std::runtime_error("inplace_matadd_colvec: dimension mismatch");
-	    }
-
-	    for (uint64_t i = 0; i < cols; ++i) {
-	        for (uint64_t j = 0; j < rows; ++j) {
-	            A(j, i) +=  x[j];
-	        }
-	    }
-	}
-
-	template <typename T>
-	void inplace_matadd_rowvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-	    if (cols != x.size()) {
-	        throw std::runtime_error("inplace_matadd_rowvec: dimension mismatch");
-	    }
-
-	    for (uint64_t i = 0; i < cols; ++i) {
-	        for (uint64_t j = 0; j < rows; ++j) {
-	            A(j, i) +=  x[i];
-	        }
-	    }
-	}
-
-	template <typename T>
-	utils::Matrix<T> matadd_colvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    //const uint64_t rows = A.rows();
-	    //const uint64_t cols = A.cols();
-
-	    utils::Matrix<T> B = A;
-
-	    inplace_matadd_colvec(B, x);
-
-	    return B;
-	}
-
-	template <typename T>
-	utils::Matrix<T> matadd_rowvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    //const uint64_t rows = A.rows();
-	    //const uint64_t cols = A.cols();
-
-	    utils::Matrix<T> B = A;
-
-	    inplace_matadd_rowvec(B, x);
-
-	    return B;
-	}
-
-	template <typename T>
-	utils::Matrix<T> matadd(const utils::Matrix<T>& A, const utils::Vector<T>& x, std::string method = "auto"){
-
-		const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-	    const uint64_t N = x.size();
-
-	    if (method=="auto"){
-	    	
-			if (rows==cols){
-				throw std::runtime_error("matadd: too many options for dimensions");
-			} else if (rows == N){
-				return matadd_rowvec(A, x);
-			} else if (cols == N){
-				return matadd_colvec(A, x);
-			}else{
-				throw std::runtime_error("matadd: undefined fault - auto");
-			}
-	    }else if(method=="row"){
-				return matadd_rowvec(A, x);
-		} else if (method=="col"){
-			return matadd_colvec(A, x);
-		}else{
-			throw std::runtime_error("matadd: undefined fault - defined method");
-		}
-    }
-
-} // namespace numerics
-
-#endif // _matadd_n_
\ No newline at end of file
diff --git a/include/numerics/matargmax.h b/include/numerics/matargmax.h
deleted file mode 100644
index 7fd411f..0000000
--- a/include/numerics/matargmax.h
+++ /dev/null
@@ -1,62 +0,0 @@
-#pragma once
-
-#include "./utils/matrix.h"
-
-namespace numerics{
-	
-	template <typename Ti, typename Td>
-	void inplace_matargmax_row(const utils::Matrix<Td>& A, utils::Vector<Ti>& b){
-
-		if (b.size() != A.rows()){
-			b.resize(A.rows(), Ti{0});
-		}
-		Td value;
-		
-		for (uint64_t i = 0; i < A.rows(); ++i){
-			value = Td{0};
-			for (uint64_t j = 0; j < A.cols(); ++j){
-				if (value < A(i,j)){
-					value = A(i,j);
-					b[i] = j;
-				}
-			}
-		}
-	}
-
-	template <typename Ti, typename Td>
-	void inplace_matargmax_col(const utils::Matrix<Td>& A, utils::Vector<Ti>& b){
-
-		if (b.size() != A.cols()){
-			b(A.cols(), Ti{0});
-		}
-		Td value;
-		
-		for (uint64_t j = 0; j < A.cols(); ++j){
-			value = Td{0};
-			for (uint64_t i = 0; i < A.cols(); ++i){
-				if (value < A(i,j)){
-					value = A(i,j);
-					b[j] = i;
-				}
-			}
-		}
-	}
-
-	template <typename Ti, typename Td>
-	utils::Vector<Ti> matargmax_row(const utils::Matrix<Td>& A){
-		utils::Vector<Ti> b(A.rows(), Ti{0});
-		inplace_matargmax_row(A, b);
-		return b;
-
-	}
-
-	template <typename Ti, typename Td>
-	utils::Vector<Ti> matargmax_col(const utils::Matrix<Td>& A){
-		utils::Vector<Ti> b(A.rows(), Ti{0});
-		inplace_matargmax_col(A, b);
-		return b;
-
-	}
-
-} // namespace numerics
-
diff --git a/include/numerics/matclip.h b/include/numerics/matclip.h
deleted file mode 100644
index 3ddc46c..0000000
--- a/include/numerics/matclip.h
+++ /dev/null
@@ -1,76 +0,0 @@
-#pragma once
-
-#include "./utils/matrix.h"
-
-
-namespace numerics{
-	
-	template <typename T>
-	void inplace_matclip_high(utils::Matrix<T>& A, T high){
-		uint64_t rows = A.rows();
-		uint64_t cols = A.cols();
-
-		for (uint64_t i = 0; i < rows; ++i){
-			for (uint64_t j = 0; j < cols; ++j){
-				if (A(i,j) > high){
-					A(i,j) = high;
-				}
-			}
-		}
-	}
-	template <typename T>
-	void inplace_matclip_low(utils::Matrix<T>& A, T low){
-		uint64_t rows = A.rows();
-		uint64_t cols = A.cols();
-
-		for (uint64_t i = 0; i < rows; ++i){
-			for (uint64_t j = 0; j < cols; ++j){
-				if (A(i,j) < low){
-					A(i,j) = low;
-				}
-			}
-		}
-	}
-	template <typename T>
-	void inplace_matclip(utils::Matrix<T>& A, T low, T high){
-		uint64_t rows = A.rows();
-		uint64_t cols = A.cols();
-
-		for (uint64_t i = 0; i < rows; ++i){
-			for (uint64_t j = 0; j < cols; ++j){
-				if (A(i,j) > high){
-					A(i,j) = high;
-				}else if (A(i,j) < low){
-					A(i,j) = low;
-				}
-			}
-		}
-	}
-	template <typename Td, typename Ti>
-	utils::Matrix<Td> matclip_high(const utils::Matrix<Ti>& A, Td high){
-
-		utils::Matrix<Td> B = A;
-		inplace_matclip_high(B, high);
-
-		return B;
-	}
-	template <typename Td, typename Ti>
-	utils::Matrix<Td> matclip_low(const utils::Matrix<Ti>& A, Td low){
-
-		utils::Matrix<Td> B = A;
-		inplace_matclip_low(B, low);
-
-		return B;
-	}
-
-	template <typename Td, typename Ti>
-	utils::Matrix<Td> matclip(const utils::Matrix<Ti>& A, Td low, Td high){
-
-		utils::Matrix<Td> B = A;
-		inplace_matclip(B, low, high);
-
-		return B;
-	}
-
-} // namespace numerics
-
diff --git a/include/numerics/matdiv.h b/include/numerics/matdiv.h
deleted file mode 100644
index c4c175e..0000000
--- a/include/numerics/matdiv.h
+++ /dev/null
@@ -1,91 +0,0 @@
-#ifndef _matdiv_n_
-#define _matdiv_n_
-
-
-#include "./utils/matrix.h"
-#include "./core/omp_config.h"
-
-
-namespace numerics{
-
-	template <typename T>
-	utils::Matrix<T> matdiv(const utils::Matrix<T>& A, const utils::Vector<T>& b, std::string method){
-        
-        utils::Matrix<T> C = A;
-
-        if (method == "row"){
-            for (uint64_t i = 0; i < A.rows(); ++i){
-                for (uint64_t j = 0; j < A.cols(); ++j){
-                    C(i,j) /= b[j];
-                }
-            }
-        }else if (method == "col"){
-            for (uint64_t i = 0; i < A.rows(); ++i){
-                for (uint64_t j = 0; j < A.cols(); ++j){
-                    C(i,j) /= b[i];
-                }
-            }
-        }else{
-            throw std::runtime_error("matdiv: choose div by: 'row' or 'col'");
-        }
-        return C;
-	}
-
-
-
-
-    template <typename T>
-    void inplace_matdiv(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_matdiv: rows and cols are not the same'");
-        }
-
-        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>
-    utils::Matrix<T> matdiv(const 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("matdiv: choose div by: 'row' or 'col'");
-        }
-
-        utils::Matrix<T> C = A;
-
-        inplace_matdiv(C, B);
-
-        return C;
-    }
-
-    template <typename T>
-    void inplace_matdiv(utils::Matrix<T>& A, const T b){
-        
-        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) /= b;
-            }
-        }
-    }
-
-
-
-
-} // namespace numerics
-
-#endif // _matdiv_n_
\ No newline at end of file
diff --git a/include/numerics/matexp.h b/include/numerics/matexp.h
deleted file mode 100644
index dd7f6be..0000000
--- a/include/numerics/matexp.h
+++ /dev/null
@@ -1,22 +0,0 @@
-#pragma once
-
-#include "./utils/matrix.h"
-
-namespace numerics{
-	
-
-	template <typename T>
-	utils::Matrix<T>  matexp(const utils::Matrix<T>& A){
-		utils::Matrix<T> B = A;
-		for (uint64_t i = 0; i < A.rows(); ++i){
-			for (uint64_t j = 0; j < A.cols(); ++j){
-				B(i,j) = numerics::exp(A(i,j));
-			}
-		}
-		return B;
-	}
-
-
-
-} // namespace numerics
-
diff --git a/include/numerics/matlog.h b/include/numerics/matlog.h
deleted file mode 100644
index 00c3469..0000000
--- a/include/numerics/matlog.h
+++ /dev/null
@@ -1,34 +0,0 @@
-#pragma once
-
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-
-
-namespace numerics{
-	
-
-	template <typename T>
-	void inplace_matlog(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){
-				numerics::inplace_log(A(i,j));
-			}
-		}
-	}
-
-	template <typename T>
-	utils::Matrix<T> log(const utils::Matrix<T>& A){
-
-		utils::Matrix<T> B = A;
-		inplace_matlog(B);
-		return B;
-	}
-
-
-} // namespace numerics
-
diff --git a/include/numerics/matmax.h b/include/numerics/matmax.h
deleted file mode 100644
index 74b2102..0000000
--- a/include/numerics/matmax.h
+++ /dev/null
@@ -1,52 +0,0 @@
-#pragma once
-
-#include "./utils/matrix.h"
-
-namespace numerics{
-	
-	template <typename T>
-	T matmax(const utils::Matrix<T>& A){
-
-		const uint64_t rows = A.rows();
-		const uint64_t cols = A.cols();
-
-		T max_value(T{0});
-
-		for (uint64_t i = 0; i < rows; ++i){
-			for (uint64_t j = 0; j < cols; ++j){
-				max_value = numerics::max(max_value, A(i,j));
-			}
-		}
-		return max_value;
-	}
-
-	template <typename T>
-	utils::Vector<T> matmax(const utils::Matrix<T>& A, std::string method){
-
-		utils::Vector<T> b;
-
-		if (method == "cols"){
-			b.resize(A.cols(), T{0});
-			for (uint64_t i = 0; i < A.cols(); ++i){
-				for (uint64_t j = 0; j < A.rows(); ++j){
-					b[i] = numerics::max(A(j, i), b[i]);
-				}
-			}
-		}else if (method == "rows"){
-			b.resize(A.rows(), T{0});
-			for (uint64_t i = 0; i < A.rows(); ++i){
-				for (uint64_t j = 0; j < A.cols(); ++j){
-					b[i] = numerics::max(A(i, j), b[i]);
-				}
-			}
-		}else{
-			throw std::runtime_error("max: choose 'rows or 'cols'");
-		}
-		return b;
-
-	}
-
-
-
-} // namespace numerics
-
diff --git a/include/numerics/matmean.h b/include/numerics/matmean.h
deleted file mode 100644
index 4a98152..0000000
--- a/include/numerics/matmean.h
+++ /dev/null
@@ -1,85 +0,0 @@
-#pragma once
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-#include "./core/omp_config.h"
-
-namespace numerics{
-	
-	template <typename T>
-	T matmean(const utils::Matrix<T>& A) {
-
-		T mean(T{0});
-	    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-
-	    for (uint64_t i = 0; i < cols; ++i) {
-	        for (uint64_t j = 0; j < rows; ++j) {
-	            mean += A(j, i);
-	        }
-	    }
-	    mean /= (static_cast<T>(rows)* static_cast<T>(cols));
-	    return mean;
-	}
-
-
-	template <typename T>
-	void inplace_matmean_row(const utils::Matrix<T>& A, utils::Vector<T>& b) {
-    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-		if (b.size() != cols){
-			b.resize(cols, T{0});
-		}
-
-    	for(uint64_t j = 0; j < cols; ++j){
-    		for (uint64_t i = 0; i < rows; ++i){
-    			b[j] += A(i, j);
-    		}
-    		b[j] /= static_cast<T>(rows);
-    	}
-	}
-
-	template <typename T>
-	void inplace_matmean_cols(const utils::Matrix<T>& A, utils::Vector<T>& b) {
-
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-		if (b.size() != rows){
-			b.resize(rows, T{0});
-		}
-		
-    	for(uint64_t i = 0; i < rows; ++i){
-    		for (uint64_t j = 0; j < cols; ++j){
-    			b[i] += A(i, j);
-    		}
-    		b[i] /= static_cast<T>(cols);
-    	}
-	}
-
-
-	template <typename T>
-	utils::Vector<T> matmean_row(const utils::Matrix<T>& A) {
-    
-	    utils::Vector<T> b(A.rows(), T{0});
-
-	    inplace_matmean_row(A, b);
-	    
-	    return b;
-	}
-
-	template <typename T>
-	utils::Vector<T> matmean_col(const utils::Matrix<T>& A) {
-
-	    utils::Vector<T> b(A.cols(), T{0});
-
-	    inplace_matmean_cols(A, b);
-	    
-	    return b;
-	}
-
-} // namespace numerics
diff --git a/include/numerics/matscalar.h b/include/numerics/matscalar.h
deleted file mode 100644
index a4c2dfb..0000000
--- a/include/numerics/matscalar.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#pragma once
-
-#include "./utils/matrix.h"
-
-namespace numerics{
-	
-	template <typename T>
-	void inplace_matscalar(utils::Matrix<T>& A, const T scalar){
-
-		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) *= scalar;
-			}
-		}
-	}
-
-	template <typename T>
-	utils::Matrix<T> matscalar(const utils::Matrix<T>& A, T scalar){
-
-		utils::Matrix<T> B = A;
-		inplace_matscalar(B, scalar);
-		return B;
-
-	}
-
-
-
-} // namespace numerics
-
diff --git a/include/numerics/matsubtract.h b/include/numerics/matsubtract.h
deleted file mode 100644
index b58e881..0000000
--- a/include/numerics/matsubtract.h
+++ /dev/null
@@ -1,102 +0,0 @@
-#ifndef _matsubtract_n_
-#define _matsubtract_n_
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-#include "./core/omp_config.h"
-
-namespace numerics{
-	
-	template <typename T>
-	void inplace_matsubtract_colvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-	    if (rows != x.size()) {
-	        throw std::runtime_error("inplace_matsubtract_colvec: dimension mismatch");
-	    }
-
-	    for (uint64_t i = 0; i < cols; ++i) {
-	        for (uint64_t j = 0; j < rows; ++j) {
-	            A(j, i) -=  x[j];
-	        }
-	    }
-	}
-
-	template <typename T>
-	void inplace_matsubtract_rowvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-
-	    if (cols != x.size()) {
-	        throw std::runtime_error("inplace_matsubtract_rowvec: dimension mismatch");
-	    }
-
-	    for (uint64_t i = 0; i < cols; ++i) {
-	        for (uint64_t j = 0; j < rows; ++j) {
-	            A(j, i) -=  x[i];
-	        }
-	    }
-	}
-
-	template <typename T>
-	utils::Matrix<T> matsubtract_colvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    //const uint64_t rows = A.rows();
-	    //const uint64_t cols = A.cols();
-
-	    utils::Matrix<T> B = A;
-
-	    inplace_matsubtract_colvec(B, x);
-
-	    return B;
-	}
-
-	template <typename T>
-	utils::Matrix<T> matsubtract_rowvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    
-	    //const uint64_t rows = A.rows();
-	    //const uint64_t cols = A.cols();
-
-	    utils::Matrix<T> B = A;
-
-	    inplace_matsubtract_rowvec(B, x);
-
-	    return B;
-	}
-
-	template <typename T>
-	utils::Matrix<T> matsubtract(const utils::Matrix<T>& A, const utils::Vector<T>& x, std::string method = "auto"){
-
-		const uint64_t rows = A.rows();
-	    const uint64_t cols = A.cols();
-	    const uint64_t N = x.size();
-
-	    if (method=="auto"){
-	    	
-			if (rows==cols){
-				throw std::runtime_error("matsubtract: too many options for dimensions");
-			} else if (rows == N){
-				return matsubtract_rowvec(A, x);
-			} else if (cols == N){
-				return matsubtract_colvec(A, x);
-			}else{
-				throw std::runtime_error("matsubtract: undefined fault - auto");
-			}
-	    }else if(method=="row"){
-				return matsubtract_rowvec(A, x);
-		} else if (method=="col"){
-			return matsubtract_colvec(A, x);
-		}else{
-			throw std::runtime_error("matsubtract: undefined fault - defined method");
-		}
-    }
-
-
-
-
-} // namespace numerics
-
-#endif // _matsubtract_n_
\ No newline at end of file
diff --git a/include/numerics/matsum.h b/include/numerics/matsum.h
deleted file mode 100644
index 1dd1f12..0000000
--- a/include/numerics/matsum.h
+++ /dev/null
@@ -1,53 +0,0 @@
-#ifndef _matsum_n_
-#define _matsum_n_
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-#include "./core/omp_config.h"
-
-namespace numerics{
-	
-	template <typename T>
-	T matsum_coeff(const utils::Matrix<T>& A) {
-
-		T b;
-
-		for (uint64_t i = 0; i < A.cols(); ++i){
-			for (uint64_t j = 0; j < A.rows(); ++j){
-				b += A(i, j); 
-			}
-		}
-		return b;
-	}
-
-
-	template <typename T>
-	utils::Vector<T> matsum(const utils::Matrix<T>& A, std::string method) {
-
-		utils::Vector<T> b;
-
-		if (method == "row"){
-			b.resize(A.cols(), T{0});
-			for (uint64_t i = 0; i < A.cols(); ++i){
-				for (uint64_t j = 0; j < A.rows(); ++j){
-					b[i] += A(j, i); 
-				}
-			}
-		}else if (method == "col"){
-			b.resize(A.rows(), T{0});
-
-			for (uint64_t i = 0; i < A.cols(); ++i){
-				for (uint64_t j = 0; j < A.rows(); ++j){
-					b[j] += A(j, i); 
-				}
-			}
-		}else{
-			throw std::runtime_error("matsum: choose sum by: 'row' or 'col'");
-		}
-		return b;
-	}
-
-
-} // namespace numerics
-
-#endif // _matadd_n_
\ No newline at end of file
diff --git a/include/numerics/max.h b/include/numerics/max.h
index 1740b2b..5756925 100644
--- a/include/numerics/max.h
+++ b/include/numerics/max.h
@@ -1,16 +1,33 @@
 #pragma once
 
-namespace numerics{
-	
-	template <typename T>
-	T max(const T a, const T b){
+#include "./core/omp_config.h"
+#include "detail/max_serial.h"
 
-		if(a < b){
-			return b;
-		}else{
-			return a;
-		}
+
+
+namespace numerics{
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	inline T max(const utils::Vector<T>& v) {
+		return detail::max_serial(v);
 	}
 
-} // namespace numerics
 
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	inline T max(const utils::Matrix<T>& A) {
+		return detail::max_serial(A);
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	inline utils::Vector<T> max_rowwise(const utils::Matrix<T>& A) {
+		return detail::max_rowwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> max_colwise(const utils::Matrix<T>& A) {
+		return detail::max_colwise_serial(A);
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/mean.h b/include/numerics/mean.h
new file mode 100644
index 0000000..7273b31
--- /dev/null
+++ b/include/numerics/mean.h
@@ -0,0 +1,33 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/mean_serial.h"
+
+
+
+namespace numerics{
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	inline T mean(const utils::Vector<T>& v) {
+		return detail::mean_serial(v);
+	}
+
+
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	inline T mean(const utils::Matrix<T>& A) {
+		return detail::mean_serial(A);
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	inline utils::Vector<T> mean_rowwise(const utils::Matrix<T>& A) {
+		return detail::mean_rowwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> mean_colwise(const utils::Matrix<T>& A) {
+		return detail::mean_colwise_serial(A);
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/min.h b/include/numerics/min.h
index f6b436a..a04d92c 100644
--- a/include/numerics/min.h
+++ b/include/numerics/min.h
@@ -1,21 +1,33 @@
 #pragma once
 
+#include "./core/omp_config.h"
+#include "detail/min_serial.h"
 
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
 
 
 namespace numerics{
-	
-	template <typename T>
-	T min(const T a, const T b){
 
-		if(a < b){
-			return a;
-		}else{
-			return b;
-		}
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	inline T min(const utils::Vector<T>& v) {
+		return detail::min_serial(v);
 	}
 
-} // namespace numerics
 
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	inline T min(const utils::Matrix<T>& A) {
+		return detail::min_serial(A);
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	inline utils::Vector<T> min_rowwise(const utils::Matrix<T>& A) {
+		return detail::min_rowwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> min_colwise(const utils::Matrix<T>& A) {
+		return detail::min_colwise_serial(A);
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/mul.h b/include/numerics/mul.h
new file mode 100644
index 0000000..7b7fcec
--- /dev/null
+++ b/include/numerics/mul.h
@@ -0,0 +1,83 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/mul_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Scalar ----------------
+	template <typename T>
+	inline void inplace_mul(utils::Matrix<T>& A, const T b) {
+		detail::inplace_mul_scalar_serial(A,b);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> mul(const utils::Matrix<T>& A, const T b) {
+	    utils::Matrix<T> out = A;
+	    inplace_mul(out, b);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_mul(utils::Vector<T>& v, const T b) {
+		detail::inplace_mul_scalar_serial(v,b);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> mul(const utils::Vector<T>& v, const T b) {
+	    utils::Vector<T> out = v;
+	    inplace_mul(out, b);
+	    return out;
+	}
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_mul(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		detail::inplace_mul_elementwise_serial(A,B);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> mul(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    utils::Matrix<T> out = A;
+	    inplace_mul(out, B);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_mul(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_mul_elementwise_serial(v,p);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> mul(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    utils::Vector<T> out = v;
+	    inplace_mul(out, p);
+	    return out;
+	}
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	inline void inplace_mul_rowwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_mul_rowwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> mul_rowwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_mul_rowwise(out, v);
+	    return out;
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	inline void inplace_mul_colwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_mul_colwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> mul_colwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_mul_colwise(out, v);
+	    return out;
+	}
+
+}
\ No newline at end of file
diff --git a/include/numerics/neg.h b/include/numerics/neg.h
new file mode 100644
index 0000000..ba6038f
--- /dev/null
+++ b/include/numerics/neg.h
@@ -0,0 +1,46 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/neg_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_neg(T& c) {
+		detail::inplace_neg_scalar_serial(c);
+	}
+	template <typename T>
+	inline T neg(const T c) {
+		T out = c;
+		inplace_neg(out);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_neg(utils::Matrix<T>& A) {
+		detail::inplace_neg_elementwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> neg(const utils::Matrix<T>& A) {
+	    utils::Matrix<T> out = A;
+	    inplace_neg(out);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_neg(utils::Vector<T>& v) {
+		detail::inplace_neg_elementwise_serial(v);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> neg(const utils::Vector<T>& v) {
+	    utils::Vector<T> out = v;
+	    inplace_neg(out);
+	    return out;
+	}
+
+
+}
\ No newline at end of file
diff --git a/include/numerics/pow.h b/include/numerics/pow.h
new file mode 100644
index 0000000..0bf1476
--- /dev/null
+++ b/include/numerics/pow.h
@@ -0,0 +1,57 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/pow_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Scalar ----------------
+	template <typename T>
+	inline void inplace_pow(utils::Matrix<T>& A, const T b) {
+		detail::inplace_pow_scalar_serial(A,b);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> pow(const utils::Matrix<T>& A, const T b) {
+	    utils::Matrix<T> out = A;
+	    inplace_pow(out, b);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_pow(utils::Vector<T>& v, const T b) {
+		detail::inplace_pow_scalar_serial(v,b);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> pow(const utils::Vector<T>& v, const T b) {
+	    utils::Vector<T> out = v;
+	    inplace_pow(out, b);
+	    return out;
+	}
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_pow(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		detail::inplace_pow_elementwise_serial(A,B);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> pow(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    utils::Matrix<T> out = A;
+	    inplace_pow(out, B);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_pow(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_pow_elementwise_serial(v,p);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> pow(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    utils::Vector<T> out = v;
+	    inplace_pow(out, p);
+	    return out;
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/sqrt.h b/include/numerics/sqrt.h
new file mode 100644
index 0000000..e03b50f
--- /dev/null
+++ b/include/numerics/sqrt.h
@@ -0,0 +1,46 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/sqrt_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_sqrt(T& c) {
+		detail::inplace_sqrt_scalar_serial(c);
+	}
+	template <typename T>
+	inline T sqrt(const T c) {
+		T out = c;
+		inplace_sqrt(out);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_sqrt(utils::Matrix<T>& A) {
+		detail::inplace_sqrt_elementwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> sqrt(const utils::Matrix<T>& A) {
+	    utils::Matrix<T> out = A;
+	    inplace_sqrt(out);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_sqrt(utils::Vector<T>& v) {
+		detail::inplace_sqrt_elementwise_serial(v);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> sqrt(const utils::Vector<T>& v) {
+	    utils::Vector<T> out = v;
+	    inplace_sqrt(out);
+	    return out;
+	}
+
+
+}
\ No newline at end of file
diff --git a/include/numerics/sub.h b/include/numerics/sub.h
new file mode 100644
index 0000000..942e46f
--- /dev/null
+++ b/include/numerics/sub.h
@@ -0,0 +1,83 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/sub_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Scalar ----------------
+	template <typename T>
+	inline void inplace_sub(utils::Matrix<T>& A, const T b) {
+		detail::inplace_sub_scalar_serial(A,b);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> sub(const utils::Matrix<T>& A, const T b) {
+	    utils::Matrix<T> out = A;
+	    inplace_sub(out, b);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_sub(utils::Vector<T>& v, const T b) {
+		detail::inplace_sub_scalar_serial(v,b);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> sub(const utils::Vector<T>& v, const T b) {
+	    utils::Vector<T> out = v;
+	    inplace_sub(out, b);
+	    return out;
+	}
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_sub(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		detail::inplace_sub_elementwise_serial(A,B);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> sub(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    utils::Matrix<T> out = A;
+	    inplace_sub(out, B);
+	    return out;
+	}
+
+	template <typename T>
+	inline void inplace_sub(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_sub_elementwise_serial(v,p);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> sub(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+	    utils::Vector<T> out = v;
+	    inplace_sub(out, p);
+	    return out;
+	}
+	// ---------------- Rowwise ----------------
+	template <typename T>
+	inline void inplace_sub_rowwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_sub_rowwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> sub_rowwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_sub_rowwise(out, v);
+	    return out;
+	}
+
+	// ---------------- Colwise ----------------
+	template <typename T>
+	inline void inplace_sub_colwise(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		detail::inplace_sub_colwise_serial(A,v);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> sub_colwise(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
+	    utils::Matrix<T> out = A;
+	    inplace_sub_colwise(out, v);
+	    return out;
+	}
+
+}
\ No newline at end of file
diff --git a/include/numerics/sum.h b/include/numerics/sum.h
new file mode 100644
index 0000000..a80ce35
--- /dev/null
+++ b/include/numerics/sum.h
@@ -0,0 +1,33 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/sum_serial.h"
+
+
+
+namespace numerics{
+
+	// ---------------- Vector -> Scalar ----------------
+	template <typename T>
+	inline T sum(const utils::Vector<T>& v) {
+		return detail::sum_serial(v);
+	}
+
+
+	// ---------------- Matrix -> Scalar ----------------
+	template <typename T>
+	inline T sum(const utils::Matrix<T>& A) {
+		return detail::sum_serial(A);
+	}
+
+	// ---------------- Matrix -> Vector ----------------
+	template <typename T>
+	inline utils::Vector<T> sum_rowwise(const utils::Matrix<T>& A) {
+		return detail::sum_rowwise_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> sum_colwise(const utils::Matrix<T>& A) {
+		return detail::sum_colwise_serial(A);
+	}
+}
\ No newline at end of file
diff --git a/include/numerics/vecargmax.h b/include/numerics/vecargmax.h
deleted file mode 100644
index f61b99d..0000000
--- a/include/numerics/vecargmax.h
+++ /dev/null
@@ -1,26 +0,0 @@
-#pragma once
-
-#include "./utils/vector.h"
-
-namespace numerics{
-	
-	template <typename T>
-	T vecarmax(const utils::Vector<T>& a){
-
-		const uint64_t N = a.size();
-		uint64_t idx = 0;
-		T value;
-
-		for (uint64_t i = 0; i < N; ++i){
-			if (value > a[i]){
-				value = a[i];
-				idx = i;
-
-			}
-		}
-		return idx;
-	}
-
-
-} // namespace numerics
-
diff --git a/include/numerics/vecclip.h b/include/numerics/vecclip.h
deleted file mode 100644
index e259c66..0000000
--- a/include/numerics/vecclip.h
+++ /dev/null
@@ -1,69 +0,0 @@
-#pragma once
-
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-
-
-namespace numerics{
-
-	template <typename T>
-	void inplace_vecclip_high(utils::Vector<T>& a, T high){
-		uint64_t N = a.size();
-
-		for (uint64_t i = 0; i < N; ++i){
-			if (a[i] > high){
-				a[i] = high;
-			}
-		}
-	}
-	template <typename T>
-	void inplace_vecclip_low(utils::Vector<T>& a, T low){
-		uint64_t N = a.size();
-
-		for (uint64_t i = 0; i < N; ++i){
-			if (a[i] < low){
-				a[i] = low;
-			}
-		}
-	}
-	template <typename T>
-	void inplace_vecclip(utils::Vector<T>& a, T low, T high){
-		uint64_t N = a.size();
-
-		for (uint64_t i = 0; i < N; ++i){
-			if (a[i] > high){
-				a[i] = high;
-			}else if (a[i] < low){
-				a[i] = low;
-			}
-		}
-	}
-	template <typename Td, typename Ti>
-	utils::Vector<Td> vecclip_high(const utils::Vector<Ti>& a, Td high){
-
-		utils::Vector<Td> b = a;
-		inplace_vecclip_high(b, high);
-
-		return b;
-	}
-	template <typename Td, typename Ti>
-	utils::Vector<Td> vecclip_low(const utils::Vector<Ti>& a, Td low){
-
-		utils::Vector<Td> b = a;
-		inplace_vecclip_low(b, low);
-
-		return b;
-	}
-	template <typename Td, typename Ti>
-	utils::Vector<Td> vecclip(const utils::Vector<Ti>& a, Td low, Td high){
-
-		utils::Vector<Td> b = a;
-		inplace_vecclip(b, low, high);
-
-		return b;
-	}
-
-
-} // namespace numerics
-
diff --git a/include/numerics/vecexp.h b/include/numerics/vecexp.h
deleted file mode 100644
index 73a9b56..0000000
--- a/include/numerics/vecexp.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#pragma once
-
-#include "./utils/vector.h"
-
-namespace numerics{
-	
-
-	template <typename T>
-	utils::Vector<T>  vecexp(const utils::Vector<T>& a){
-		utils::Vector<T> b = a;
-		for (uint64_t i = 0; i < a.size(); ++i){
-			b[i] =  numerics::exp(a[i]);
-		}
-		return b;
-	}
-
-
-} // namespace numerics
-
diff --git a/include/numerics/veclog.h b/include/numerics/veclog.h
deleted file mode 100644
index 53560f5..0000000
--- a/include/numerics/veclog.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#pragma once
-
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-
-
-namespace numerics{
-	
-
-	template <typename T>
-	void inplace_veclog(utils::Vector<T>& a){
-
-		const uint64_t N = a.size();
-
-		for (uint64_t i = 0; i < N; ++i){
-			numerics::inplace_log(a[i]);
-		}
-	}
-
-	template <typename T>
-	utils::Vector<T> veclog(const utils::Vector<T>& a){
-
-		utils::Vector<T> b = a;
-		inplace_veclog(b);
-		return b;
-	}
-
-} // namespace numerics
-
diff --git a/include/numerics/vecmax.h b/include/numerics/vecmax.h
deleted file mode 100644
index e22f42d..0000000
--- a/include/numerics/vecmax.h
+++ /dev/null
@@ -1,22 +0,0 @@
-#pragma once
-
-
-#include "./utils/vector.h"
-
-namespace numerics{
-	
-
-
-	template <typename T>
-	T vecmax(const utils::Vector<T>& a){
-		T max_value(T{0});
-		uint64_t N = a.size();
-
-		for (uint64_t i = 0; i < N; ++i){
-			max_value = numerics::max(max_value, a[i]);
-		}
-		return max_value;
-	}
-
-} // namespace numerics
-
diff --git a/include/numerics/vecmean.h b/include/numerics/vecmean.h
deleted file mode 100644
index 65b0b7a..0000000
--- a/include/numerics/vecmean.h
+++ /dev/null
@@ -1,66 +0,0 @@
-#ifndef _mean_n_
-#define _mean_n_
-
-#include "./utils/vector.h"
-#include "./utils/matrix.h"
-#include "./core/omp_config.h"
-
-namespace numerics{
-	
-	template <typename T>
-	T vecmean(utils::Vector<T>& a) {
-
-		T mean(T{0});
-
-		const uint64_t N = a.size();
-
-
-	    for (uint64_t i = 0; i < N; ++i) {
-            mean += a[i];
-	    }
-	    mean /= (static_cast<T>(N));
-	    return mean;
-	}
-
-	template <typename T>
-	T vecmean_isequal(utils::Vector<T>& a, utils::Vector<T>& b, T tol=1e-12) {
-
-		uint64_t count = T{0};
-
-		const uint64_t N = a.size();
-
-		if (a.size() != b.size()){
-			throw std::runtime_error("vecmean_equal: dimension mismatch");
-		}
-
-
-	    for (uint64_t i = 0; i < N; ++i) {
-	    	if ((a[i]-b[i]) < tol){
-	    		count += 1;
-	    	}
-	    }
-	    return static_cast<T>(count)/static_cast<T>(N);
-	}
-	template <typename Td, typename Ti>
-	Td vecmean_equal(utils::Vector<Ti>& a, utils::Vector<Ti>& b) {
-
-		Ti count = Ti{0};
-
-		const uint64_t N = a.size();
-
-		if (a.size() != b.size()){
-			throw std::runtime_error("vecmean_equal: dimension mismatch");
-		}
-
-
-	    for (uint64_t i = 0; i < N; ++i) {
-	    	if (a[i]==b[i]){
-	    		count += Ti{1};
-	    	}
-	    }
-	    return static_cast<Td>(count)/static_cast<Td>(N);
-	}
-
-} // namespace numerics
-
-#endif // _mean_n_
\ No newline at end of file