diff --git a/Fastor/tensor/TensorFunctions.h b/Fastor/tensor/TensorFunctions.h
index da75710..235d562 100644
--- a/Fastor/tensor/TensorFunctions.h
+++ b/Fastor/tensor/TensorFunctions.h
@@ -124,8 +124,10 @@ FASTOR_INLINE Tensor<T,Rest...> torowmajor(const TensorType<T,Rest...> &a) {
     }
 }
 
-/* Squeeze - removes dimenions of 1 from a tensor and returns a map
-    Note that you cannot use tensor expressions on squeeze as squeeze is
+/* squeeze - removes dimenions of 1 from a tensor and returns a TensorMap
+    A TensorMap is a view in to an existing tensor so modifying the squeezed
+    tensor will modify the original tensor and vice versa
+    Note that you cannot call this function on tensor expressions as this is
     a view in to a concrete tensor type holding storage
 */
 template<template <typename,size_t...> class TensorType, typename T, size_t ... Rest>
@@ -135,7 +137,34 @@ squeeze(const TensorType<T,Rest...> &a) {
     return index_to_tensor_map_t<T,filter_t<1,Rest...>>(a.data());
 }
 
+/* reshape - reshapes a tensor to a tensor of different shape and returns a TensorMap
+    A TensorMap is a view in to an existing tensor so modifying the reshaped
+    tensor will modify the original tensor and vice versa
+    Note that you cannot call this function on tensor expressions as this is
+    a view in to a concrete tensor type holding storage
 
+    example:
+        auto b = reshape<shapes...>(a);
+*/
+template<size_t ... shapes,typename T, size_t ... Rest>
+FASTOR_INLINE TensorMap<T,shapes...> reshape(const Tensor<T,Rest...> &a) {
+    static_assert(pack_prod<shapes...>::value==pack_prod<Rest...>::value, "SIZE OF TENSOR SHOULD REMAIN THE SAME DURING RESHAPE");
+    return TensorMap<T,shapes...>(a.data());
+}
+
+/* flatten - creates a flattened 1D view of a tensor and returns a TensorMap
+    A TensorMap is a view in to an existing tensor so modifying the reshaped
+    tensor will modify the original tensor and vice versa
+    Note that you cannot call this function on tensor expressions as this is
+    a view in to a concrete tensor type holding storage
+
+    example:
+        auto b = flatten(a);
+*/
+template<typename T, size_t ... Rest>
+FASTOR_INLINE TensorMap<T,pack_prod<Rest...>::value> flatten(const Tensor<T,Rest...> &a) {
+    return TensorMap<T,pack_prod<Rest...>::value>(a.data());
+}
 
 
 
diff --git a/Fastor/tensor_algebra/einsum.h b/Fastor/tensor_algebra/einsum.h
index cfef5ef..910cbc7 100644
--- a/Fastor/tensor_algebra/einsum.h
+++ b/Fastor/tensor_algebra/einsum.h
@@ -7,7 +7,6 @@
 #include "Fastor/tensor_algebra/indicial.h"
 #include "Fastor/backend/voigt.h"
 
-#include "Fastor/tensor_algebra/reshape.h"
 #include "Fastor/tensor_algebra/permutation.h"
 #include "Fastor/tensor_algebra/permute.h"
 #include "Fastor/tensor_algebra/innerproduct.h"
diff --git a/Fastor/tensor_algebra/reshape.h b/Fastor/tensor_algebra/reshape.h
deleted file mode 100644
index 8efcfd6..0000000
--- a/Fastor/tensor_algebra/reshape.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef RESHAPE_H
-#define RESHAPE_H
-
-#include "Fastor/tensor/Tensor.h"
-#include "Fastor/tensor_algebra/indicial.h"
-#include "Fastor/tensor_algebra/permutation.h"
-
-namespace Fastor {
-
-
-template<size_t ... shapes,typename T, size_t ... Rest>
-FASTOR_INLINE Tensor<T,shapes...> reshape(const Tensor<T,Rest...> &a) {
-    // Call it as reshape<shapes...>(a)
-    static_assert(pack_prod<shapes...>::value==pack_prod<Rest...>::value, "SIZE OF TENSOR SHOULD REMAIN THE SAME DURING RESHAPE");
-    Tensor<T,shapes...> out;
-    //! Reshape involves deep copy, as in place permutation is not possible since Tensor is static
-    std::copy(a.data(),a.data()+a.Size,out.data());
-    return out;
-  }
-
-
-// flatten
-template<typename T, size_t ... Rest>
-FASTOR_INLINE Tensor<T,pack_prod<Rest...>::value> flatten(const Tensor<T,Rest...> &a) {
-    Tensor<T,pack_prod<Rest...>::value> out;
-    std::copy(a.data(),a.data()+a.Size,out.data());
-    return out;
-}
-
-}
-#endif // RESHAPE_H
-
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index f0501d1..633edab 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -22,6 +22,10 @@ add_subdirectory(test_simd_vectors)
 
 add_subdirectory(test_tensor_basics)
 
+add_subdirectory(test_tensormap)
+
+add_subdirectory(test_numerics)
+
 add_subdirectory(test_math_functions)
 
 add_subdirectory(test_booleans)
@@ -55,8 +59,6 @@ add_subdirectory(test_mixed_views)
 
 add_subdirectory(test_complex_expressions)
 
-add_subdirectory(test_tensormap)
-
 add_subdirectory(test_unary_bool_ops)
 
 add_subdirectory(test_binary_cmp_ops)
\ No newline at end of file
diff --git a/tests/test_booleans/CMakeLists.txt b/tests/test_booleans/CMakeLists.txt
index 4b69534..4d0788e 100644
--- a/tests/test_booleans/CMakeLists.txt
+++ b/tests/test_booleans/CMakeLists.txt
@@ -7,7 +7,7 @@ add_executable(test_booleans test_booleans.cpp)
 add_test(test_booleans test_booleans)
 
 if(MSVC)
-    add_compile_options(test_booleans PRIVATE "/W4" "$<$<CONFIG:RELEASE>:/O2>")
+    add_compile_options(test_booleans PRIVATE "/W2" "$<$<CONFIG:RELEASE>:/O2>")
 else()
     add_compile_options(test_booleans PRIVATE "$<$<CONFIG:RELEASE>:-O3>" "$<$<CONFIG:RELEASE>:-march=native>")
 endif()
diff --git a/tests/test_numerics/CMakeLists.txt b/tests/test_numerics/CMakeLists.txt
new file mode 100644
index 0000000..0005cb7
--- /dev/null
+++ b/tests/test_numerics/CMakeLists.txt
@@ -0,0 +1,16 @@
+cmake_minimum_required(VERSION 3.1)
+project(test_numerics)
+
+set(CMAKE_CXX_STANDARD 14)
+
+add_executable(test_numerics test_numerics.cpp)
+add_test(test_numerics test_numerics)
+
+if(MSVC)
+    add_compile_options(test_numerics PRIVATE "/W2" "$<$<CONFIG:RELEASE>:/O2>")
+else()
+    add_compile_options(test_numerics PRIVATE "$<$<CONFIG:RELEASE>:-O3>" "$<$<CONFIG:RELEASE>:-march=native>")
+endif()
+
+target_include_directories(test_numerics PRIVATE ${FASTOR_INCLUDE_DIR})
+target_include_directories(test_numerics PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../)
\ No newline at end of file
diff --git a/tests/test_numerics/test_numerics.cpp b/tests/test_numerics/test_numerics.cpp
new file mode 100644
index 0000000..4535e89
--- /dev/null
+++ b/tests/test_numerics/test_numerics.cpp
@@ -0,0 +1,71 @@
+#include <Fastor/Fastor.h>
+
+using namespace Fastor;
+
+
+#define Tol 1e-12
+#define BigTol 1e-5
+
+
+template<typename T>
+void test_numerics() {
+
+    // The test assumes that TensorMap is already
+    // tested elsewhere
+
+    // reshape
+    {
+        Tensor<T,3,4> a; a.iota(1);
+
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - reshape<4,3>(a).sum()) < Tol);
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - reshape<2,2,3>(a).sum()) < Tol);
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - reshape<1,2,2,3>(a).sum()) < Tol);
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - reshape<2,3,2>(a).sum()) < Tol);
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - reshape<3,2,2>(a).sum()) < Tol);
+
+        FASTOR_EXIT_ASSERT(reshape<4,3>(a).dimension(0) - 4 < Tol);
+        FASTOR_EXIT_ASSERT(reshape<4,3>(a).dimension(1) - 3 < Tol);
+
+        FASTOR_EXIT_ASSERT(reshape<2,2,3>(a).dimension(0) - 2 < Tol);
+        FASTOR_EXIT_ASSERT(reshape<2,2,3>(a).dimension(1) - 2 < Tol);
+        FASTOR_EXIT_ASSERT(reshape<2,2,3>(a).dimension(2) - 3 < Tol);
+    }
+
+    // flatten
+    {
+        Tensor<T,3,4> a; a.iota(1);
+
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - flatten(a).sum()) < Tol);
+
+        FASTOR_EXIT_ASSERT(flatten(a).dimension(0) - 12 < Tol);
+    }
+
+    // squeeze
+    {
+        Tensor<T,1,3,4> a; a.iota(1);
+
+        FASTOR_EXIT_ASSERT(std::abs(a.sum() - squeeze(a).sum()) < Tol);
+
+        FASTOR_EXIT_ASSERT(squeeze(a).dimension(0) - 3 < Tol);
+        FASTOR_EXIT_ASSERT(squeeze(a).dimension(1) - 4 < Tol);
+
+        FASTOR_EXIT_ASSERT(squeeze(Tensor<T,1,2,1>{}).dimension(0) - 2 < Tol);
+        FASTOR_EXIT_ASSERT(is_same_v_<decltype(squeeze(Tensor<T,1,1>{})),TensorMap<T>> == true);
+        FASTOR_EXIT_ASSERT(is_same_v_<decltype(squeeze(Tensor<T,1,1,3>{})),TensorMap<T,3>> == true);
+    }
+
+    print(FGRN(BOLD("All tests passed successfully")));
+
+}
+
+
+int main() {
+
+    print(FBLU(BOLD("Testing numerics with single precision")));
+    test_numerics<float>();
+    print(FBLU(BOLD("Testing numerics with double precision")));
+    test_numerics<double>();
+
+    return 0;
+}
+