diff --git a/src/benchmarks/Sofa.Type/Matrix.cpp b/src/benchmarks/Sofa.Type/Matrix.cpp
index 0c8e1e9..655ce73 100644
--- a/src/benchmarks/Sofa.Type/Matrix.cpp
+++ b/src/benchmarks/Sofa.Type/Matrix.cpp
@@ -7,188 +7,258 @@
 #include <Eigen/LU> // for inverse(?)
 #include <array>
 
-static void BM_Matrix_typemat3x3f_construct(benchmark::State& state);
-static void BM_Matrix_typemat3x3f_construct_noinit(benchmark::State& state);
-static void BM_Matrix_eigenmat33_construct(benchmark::State& state);
-static void BM_Matrix_typemat3x3f_transpose(benchmark::State& state);
-static void BM_Matrix_eigenmat33_transpose(benchmark::State& state);
-static void BM_Matrix_typemat3x3f_invert(benchmark::State& state);
-static void BM_Matrix_eigenmat33_invert(benchmark::State& state);
-static void BM_Matrix_typemat3x3f_determinant(benchmark::State& state);
-static void BM_Matrix_eigenmat33_determinant(benchmark::State& state);
-template<int N>
-static void BM_Matrix_typematf_matmult(benchmark::State& state);
-template<int N>
-static void BM_Matrix_eigenmatf_matmult(benchmark::State& state);
-template<int N>
-static void BM_Matrix_typematf_vecmult(benchmark::State& state);
-template<int N>
-static void BM_Matrix_eigenmatf_vecmult(benchmark::State& state);
-static void BM_Matrix_typemat3x3f_assign(benchmark::State& state);
-static void BM_Matrix_eigenmat33_assign(benchmark::State& state);
-static void BM_Matrix_typemat3x3f_multTranspose(benchmark::State& state);
-static void BM_Matrix_eigenmat3x3f_multTranspose(benchmark::State& state);
-
-constexpr int64_t minSubIterations = 8 << 4;
-constexpr int64_t maxSubIterations = 8 << 6;
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_construct(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_construct_noinit(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_transpose(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_invert(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_determinant(benchmark::State& state);
+template<typename ScalarType, int N>
+static void BM_Matrix_typemat_matmult(benchmark::State& state);
+template<typename ScalarType, int N>
+static void BM_Matrix_typemat_vecmult(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_assign(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_typemat3x3_multTranspose(benchmark::State& state);
+
+// eigen version
+template<typename ScalarType>
+static void BM_Matrix_eigenmat3x3_construct(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_eigenmat3x3_transpose(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_eigenmat3x3_invert(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_eigenmat3x3_determinant(benchmark::State& state);
+template<typename ScalarType, int N>
+static void BM_Matrix_eigenmat_matmult(benchmark::State& state);
+template<typename ScalarType, int N>
+static void BM_Matrix_eigenmat_vecmult(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_eigenmat3x3_assign(benchmark::State& state);
+template<typename ScalarType>
+static void BM_Matrix_eigenmat3x3_multTranspose(benchmark::State& state);
+
+constexpr int64_t minSubIterations = 8 << 6;
+constexpr int64_t maxSubIterations = 8 << 8;
 
 #define BMARGS ->RangeMultiplier(2)->Ranges({ {minSubIterations, maxSubIterations} })->Unit(benchmark::kMicrosecond)
 
-BENCHMARK(BM_Matrix_typemat3x3f_construct) BMARGS;
-BENCHMARK(BM_Matrix_typemat3x3f_construct_noinit) BMARGS;
-BENCHMARK(BM_Matrix_eigenmat33_construct) BMARGS;
-BENCHMARK(BM_Matrix_typemat3x3f_transpose) BMARGS;
-BENCHMARK(BM_Matrix_eigenmat33_transpose) BMARGS;
-BENCHMARK(BM_Matrix_typemat3x3f_invert) BMARGS;
-BENCHMARK(BM_Matrix_eigenmat33_invert) BMARGS;
-BENCHMARK(BM_Matrix_typemat3x3f_determinant) BMARGS;
-BENCHMARK(BM_Matrix_eigenmat33_determinant) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_typematf_matmult, 3) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_typematf_matmult, 24) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_eigenmatf_matmult, 3) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_eigenmatf_matmult, 24) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_typematf_vecmult, 3) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_typematf_vecmult, 24) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_eigenmatf_vecmult, 3) BMARGS;
-BENCHMARK_TEMPLATE(BM_Matrix_eigenmatf_vecmult, 24) BMARGS;
-BENCHMARK(BM_Matrix_typemat3x3f_assign) BMARGS;
-BENCHMARK(BM_Matrix_eigenmat33_assign) BMARGS;
-BENCHMARK(BM_Matrix_typemat3x3f_multTranspose) BMARGS;
-BENCHMARK(BM_Matrix_eigenmat3x3f_multTranspose) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_construct, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_construct, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_construct_noinit, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_construct_noinit, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_transpose, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_transpose, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_invert, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_invert, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_determinant,float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_determinant,double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_matmult, float, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_matmult, double, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_matmult, float, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_matmult, double, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_matmult, float, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_matmult, double, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_vecmult, float, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_vecmult, double, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_vecmult, float, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_vecmult, double, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_vecmult, float, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat_vecmult, double, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_assign, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_assign, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_multTranspose, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_typemat3x3_multTranspose, double) BMARGS;
+
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_construct, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_construct, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_transpose, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_transpose, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_invert, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_invert, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_determinant, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_determinant, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_matmult, float, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_matmult, double, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_matmult, float, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_matmult, double, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_matmult, float, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_matmult, double, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_vecmult, float, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_vecmult, double, 3) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_vecmult, float, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_vecmult, double, 6) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_vecmult, float, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat_vecmult, double, 24) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_assign, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_assign, double) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_multTranspose, float) BMARGS;
+BENCHMARK_TEMPLATE(BM_Matrix_eigenmat3x3_multTranspose, double) BMARGS;
 
 #undef BMARGS
 
-void BM_Matrix_typemat3x3f_construct(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_construct(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using Line = typename Matrix::Line;
+    
     for (auto _ : state)
     {
         for (unsigned int i = 0; i < state.range(0); ++i)
         {
-            benchmark::DoNotOptimize(sofa::type::Mat3x3f{
-                sofa::type::Vec<3, float>{values[i], values[i + 1], values[i + 2]},
-                sofa::type::Vec<3, float>{values[i + 3], values[i + 4], values[i + 5]},
-                sofa::type::Vec<3, float>{values[i + 6], values[i + 7], values[i + 8] }
+            benchmark::DoNotOptimize(Matrix{
+                Line{values[i], values[i + 1], values[i + 2]},
+                Line{values[i + 3], values[i + 4], values[i + 5]},
+                Line{values[i + 6], values[i + 7], values[i + 8] }
             });
         }
     }
 }
 
-void BM_Matrix_typemat3x3f_construct_noinit(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_construct_noinit(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using Line = typename Matrix::LineNoInit;
+    
     for (auto _ : state)
     {
         for (unsigned int i = 0; i < state.range(0); ++i)
         {
-            benchmark::DoNotOptimize(sofa::type::Mat3x3f{
-                sofa::type::Mat3x3f::LineNoInit{values[i], values[i + 1], values[i + 2]},
-                sofa::type::Mat3x3f::LineNoInit{values[i + 3], values[i + 4], values[i + 5]},
-                sofa::type::Mat3x3f::LineNoInit{values[i + 6], values[i + 7], values[i + 8] }
+            benchmark::DoNotOptimize(Matrix{
+                Line{values[i], values[i + 1], values[i + 2]},
+                Line{values[i + 3], values[i + 4], values[i + 5]},
+                Line{values[i + 6], values[i + 7], values[i + 8] }
             });
         }
     }
 }
 
-void BM_Matrix_eigenmat33_construct(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_eigenmat3x3_construct(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
+    using Matrix = Eigen::Matrix<ScalarType, 3, 3>;
+    
     for (auto _ : state)
     {
         for (unsigned int i = 0; i < state.range(0); ++i)
         {
             //only eigen3.4
-            //Eigen::Matrix<float, 3,3>(
+            //Eigen::Matrix<ScalarType, 3,3>(
             //    { vc[i], vc[i + 1], vc[i + 2] },
             //    { vc[i + 3], vc[i + 4], vc[i + 5] },
             //    { vc[i + 6], vc[i + 7], vc[i + 8] }
             //);
 
-            Eigen::Matrix<float, 3, 3> m;
+            Matrix m;
             m << values[i], values[i + 1], values[i + 2], values[i + 3], values[i + 4], values[i + 5], values[i + 6], values[i + 7], values[i + 8];
             benchmark::DoNotOptimize(m);
         }
     }
 }
 
-void BM_Matrix_typemat3x3f_assign(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_assign(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using MatrixNoInit = sofa::type::MatNoInit<3,3,ScalarType>;
+    using Line = typename Matrix::LineNoInit;
+    
     for (auto _ : state)
     {
         for (unsigned int i = 0; i < state.range(0); ++i)
         {
-            sofa::type::Mat3x3f m1{
-                sofa::type::Mat3x3f::LineNoInit{values[i], values[i + 1], values[i + 2]},
-                sofa::type::Mat3x3f::LineNoInit{values[i + 3], values[i + 4], values[i + 5]},
-                sofa::type::Mat3x3f::LineNoInit{values[i + 6], values[i + 7], values[i + 8] }
+            Matrix m1{
+                Line{values[i], values[i + 1], values[i + 2]},
+                Line{values[i + 3], values[i + 4], values[i + 5]},
+                Line{values[i + 6], values[i + 7], values[i + 8] }
             };
-            sofa::type::Mat3x3f m2;
+            MatrixNoInit m2{};
             m2 = m1;
             benchmark::DoNotOptimize(m2);
         }
     }
 }
 
-void BM_Matrix_eigenmat33_assign(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_eigenmat3x3_assign(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
+    using Matrix = Eigen::Matrix<ScalarType, 3, 3>;
+    
     for (auto _ : state)
     {
         for (unsigned int i = 0; i < state.range(0); ++i)
         {
-            Eigen::Matrix<float, 3, 3> m1;
+            Matrix m1;
             m1 << values[i], values[i + 1], values[i + 2], values[i + 3], values[i + 4], values[i + 5], values[i + 6], values[i + 7], values[i + 8];
-            Eigen::Matrix<float, 3, 3> m2;
+            Matrix m2;
             m2 = m1;
             benchmark::DoNotOptimize(m2);
         }
     }
 }
 
-void BM_Matrix_typemat3x3f_transpose(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_transpose(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<sofa::type::Mat3x3f> vc;
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using Line = typename Matrix::LineNoInit;
+    
+    std::vector<Matrix> vc;
     vc.reserve(state.range(0));
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        vc.push_back(sofa::type::Mat3x3f{
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 0], values[i * 9 + 1], values[i * 9 + 2]},
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 3], values[i * 9 + 4], values[i * 9 + 5]},
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 6], values[i * 9 + 7], values[i * 9 + 8]}
-            });
+        vc.emplace_back(
+            Line{values[i * 9 + 0], values[i * 9 + 1], values[i * 9 + 2]},
+            Line{values[i * 9 + 3], values[i * 9 + 4], values[i * 9 + 5]},
+            Line{values[i * 9 + 6], values[i * 9 + 7], values[i * 9 + 8]}
+        );
     }
 
     for (auto _ : state)
     {
         for (auto& mat : vc)
         {
-            auto t = mat.transposed();
-            benchmark::DoNotOptimize(t);
+            benchmark::DoNotOptimize(mat.transposed());
         }
     }
 }
 
-void BM_Matrix_eigenmat33_transpose(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_eigenmat3x3_transpose(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<Eigen::Matrix3f> vc;
+    using Matrix = Eigen::Matrix<ScalarType, 3, 3>;
+    
+    std::vector<Matrix> vc;
     vc.resize(state.range(0));
 
     for (unsigned int i = 0; i < state.range(0); i++)
@@ -202,20 +272,23 @@ void BM_Matrix_eigenmat33_transpose(benchmark::State& state)
     {
         for (auto& mat : vc)
         {
-            Eigen::Matrix3f t = mat.transpose();
-            benchmark::DoNotOptimize(t);
+            Matrix res = mat.transpose();
+            benchmark::DoNotOptimize(res);
         }
     }
 }
 
-template<int N>
-void BM_Matrix_typematf_matmult(benchmark::State& state)
+template<typename ScalarType, int N>
+void BM_Matrix_typemat_matmult(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * N*N * 2;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<sofa::type::Mat<N, N, float > > vc1;
-    std::vector<sofa::type::Mat<N, N, float > > vc2;
+    using Matrix = sofa::type::Mat<N,N,ScalarType>;
+    using MatrixNoInit = sofa::type::MatNoInit<N,N,ScalarType>;
+    
+    std::vector<Matrix> vc1;
+    std::vector<Matrix> vc2;
     vc1.reserve(state.range(0));
     vc2.reserve(state.range(0));
 
@@ -223,7 +296,7 @@ void BM_Matrix_typematf_matmult(benchmark::State& state)
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        sofa::type::Mat<N, N, float > mat1, mat2;
+        MatrixNoInit mat1, mat2;
         for (int a = 0; a < N; ++a)
         {
             for (int b = 0; b < N; ++b)
@@ -245,14 +318,16 @@ void BM_Matrix_typematf_matmult(benchmark::State& state)
     }
 }
 
-template<int N>
-void BM_Matrix_eigenmatf_matmult(benchmark::State& state)
+template<typename ScalarType, int N>
+void BM_Matrix_eigenmat_matmult(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * N*N * 2;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<Eigen::Matrix<float, N, N > > vc1;
-    std::vector<Eigen::Matrix<float, N, N > > vc2;
+    using Matrix = Eigen::Matrix<ScalarType, N, N>;
+    
+    std::vector<Matrix> vc1;
+    std::vector<Matrix> vc2;
     vc1.reserve(state.range(0));
     vc2.reserve(state.range(0));
     benchmark::DoNotOptimize(vc1);
@@ -262,7 +337,7 @@ void BM_Matrix_eigenmatf_matmult(benchmark::State& state)
 
     for (auto i = 0; i < state.range(0); i++)
     {
-        Eigen::Matrix<float, N, N > mat1, mat2;
+        Matrix mat1, mat2;
         for (int a = 0; a < N; ++a)
         {
             for (int b = 0; b < N; ++b)
@@ -279,20 +354,25 @@ void BM_Matrix_eigenmatf_matmult(benchmark::State& state)
     {
         for (auto i = 0; i < state.range(0); i++)
         {
-            Eigen::Matrix<float, N, N > product = vc1[i] * vc2[i];
-            benchmark::DoNotOptimize(product);
+            Matrix res = vc1[i] * vc2[i];
+            benchmark::DoNotOptimize(res);
         }
     }
 }
 
-template <int N>
-void BM_Matrix_typematf_vecmult(benchmark::State& state)
+template<typename ScalarType, int N>
+void BM_Matrix_typemat_vecmult(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * (N*N + N);
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
-
-    std::vector<sofa::type::Mat<N, N, float > > vc1;
-    std::vector<sofa::type::Vec<N, float > > vc2;
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
+
+    using Matrix = sofa::type::Mat<N,N,ScalarType>;
+    using Vec = sofa::type::Vec<N,ScalarType>;
+    using MatrixNoInit = sofa::type::MatNoInit<N,N,ScalarType>;
+    using VecNoInit = sofa::type::VecNoInit<N,ScalarType>;
+    
+    std::vector<Matrix> vc1;
+    std::vector<Vec> vc2;
     vc1.reserve(state.range(0));
     vc2.reserve(state.range(0));
 
@@ -300,8 +380,8 @@ void BM_Matrix_typematf_vecmult(benchmark::State& state)
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        sofa::type::Mat<N, N, float > mat;
-        sofa::type::Vec<N, float > vec;
+        MatrixNoInit mat;
+        VecNoInit vec;
         for (int a = 0; a < N; ++a)
         {
             for (int b = 0; b < N; ++b)
@@ -323,14 +403,17 @@ void BM_Matrix_typematf_vecmult(benchmark::State& state)
     }
 }
 
-template <int N>
-void BM_Matrix_eigenmatf_vecmult(benchmark::State& state)
+template<typename ScalarType, int N>
+void BM_Matrix_eigenmat_vecmult(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * (N*N + N);
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<Eigen::Matrix<float, N, N > > vc1;
-    std::vector<Eigen::Matrix<float, N, 1 > > vc2;
+    using Matrix = Eigen::Matrix<ScalarType, N, N>;
+    using Vec = Eigen::Matrix<ScalarType, N, 1>;
+    
+    std::vector<Matrix> vc1;
+    std::vector<Vec> vc2;
     vc1.reserve(state.range(0));
     vc2.reserve(state.range(0));
     benchmark::DoNotOptimize(vc1);
@@ -340,8 +423,8 @@ void BM_Matrix_eigenmatf_vecmult(benchmark::State& state)
 
     for (auto i = 0; i < state.range(0); i++)
     {
-        Eigen::Matrix<float, N, N > mat1;
-        Eigen::Matrix<float, N, 1 > mat2;
+        Matrix mat1;
+        Vec mat2;
         for (int a = 0; a < N; ++a)
         {
             for (int b = 0; b < N; ++b)
@@ -358,28 +441,31 @@ void BM_Matrix_eigenmatf_vecmult(benchmark::State& state)
     {
         for (auto i = 0; i < state.range(0); i++)
         {
-            Eigen::Matrix<float, N, 1 > product = vc1[i] * vc2[i];
-            benchmark::DoNotOptimize(product);
+            Vec res = vc1[i] * vc2[i];
+            benchmark::DoNotOptimize(res);
         }
     }
 }
 
-
-void BM_Matrix_typemat3x3f_determinant(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_determinant(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<sofa::type::Mat3x3f> vc;
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using Line = typename Matrix::LineNoInit;
+    
+    std::vector<Matrix> vc;
     vc.reserve(state.range(0));
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        vc.push_back(sofa::type::Mat3x3f{
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 0], values[i * 9 + 1], values[i * 9 + 2]},
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 3], values[i * 9 + 4], values[i * 9 + 5]},
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 6], values[i * 9 + 7], values[i * 9 + 8]}
-            });
+        vc.emplace_back(
+            Line{values[i * 9 + 0], values[i * 9 + 1], values[i * 9 + 2]},
+            Line{values[i * 9 + 3], values[i * 9 + 4], values[i * 9 + 5]},
+            Line{values[i * 9 + 6], values[i * 9 + 7], values[i * 9 + 8]}
+        );
     }
 
     for (auto _ : state)
@@ -391,12 +477,15 @@ void BM_Matrix_typemat3x3f_determinant(benchmark::State& state)
     }
 }
 
-void BM_Matrix_eigenmat33_determinant(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_eigenmat3x3_determinant(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<Eigen::Matrix3f> vc;
+    using Matrix = Eigen::Matrix<ScalarType, 3, 3>;
+    
+    std::vector<Matrix> vc;
     vc.resize(state.range(0));
 
     for (unsigned int i = 0; i < state.range(0); i++)
@@ -415,22 +504,25 @@ void BM_Matrix_eigenmat33_determinant(benchmark::State& state)
     }
 }
 
-
-void BM_Matrix_typemat3x3f_invert(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_invert(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<sofa::type::Mat3x3f> vc;
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using Line = typename Matrix::LineNoInit;
+    
+    std::vector<Matrix> vc;
     vc.reserve(state.range(0));
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        vc.push_back(sofa::type::Mat3x3f{
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 0], values[i * 9 + 1], values[i * 9 + 2]},
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 3], values[i * 9 + 4], values[i * 9 + 5]},
-                sofa::type::Mat3x3f::LineNoInit{values[i * 9 + 6], values[i * 9 + 7], values[i * 9 + 8]}
-            });
+        vc.emplace_back(
+            Line{values[i * 9 + 0], values[i * 9 + 1], values[i * 9 + 2]},
+            Line{values[i * 9 + 3], values[i * 9 + 4], values[i * 9 + 5]},
+            Line{values[i * 9 + 6], values[i * 9 + 7], values[i * 9 + 8]}
+        );
     }
 
     for (auto _ : state)
@@ -442,12 +534,15 @@ void BM_Matrix_typemat3x3f_invert(benchmark::State& state)
     }
 }
 
-void BM_Matrix_eigenmat33_invert(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_eigenmat3x3_invert(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 9;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<Eigen::Matrix3f> vc;
+    using Matrix = Eigen::Matrix<ScalarType, 3, 3>;
+    
+    std::vector<Matrix> vc;
     vc.resize(state.range(0));
 
     for (unsigned int i = 0; i < state.range(0); i++)
@@ -461,19 +556,23 @@ void BM_Matrix_eigenmat33_invert(benchmark::State& state)
     {
         for (auto& mat : vc)
         {
-            Eigen::Matrix3f inv = mat.inverse();
-            benchmark::DoNotOptimize(inv);
+            Matrix res = mat.inverse();
+            benchmark::DoNotOptimize(res);
         }
     }
 }
 
-void BM_Matrix_typemat3x3f_multTranspose(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_typemat3x3_multTranspose(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 3*3 * 2;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<sofa::type::Mat<3, 3, float > > vc1;
-    std::vector<sofa::type::Mat<3, 3, float > > vc2;
+    using Matrix = sofa::type::Mat<3,3,ScalarType>;
+    using MatrixNoInit = sofa::type::MatNoInit<3,3,ScalarType>;
+    
+    std::vector<Matrix> vc1;
+    std::vector<Matrix> vc2;
     vc1.reserve(state.range(0));
     vc2.reserve(state.range(0));
 
@@ -481,7 +580,7 @@ void BM_Matrix_typemat3x3f_multTranspose(benchmark::State& state)
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        sofa::type::Mat<3, 3, float > mat1, mat2;
+        MatrixNoInit mat1, mat2;
         for (int a = 0; a < 3; ++a)
         {
             for (int b = 0; b < 3; ++b)
@@ -503,13 +602,16 @@ void BM_Matrix_typemat3x3f_multTranspose(benchmark::State& state)
     }
 }
 
-void BM_Matrix_eigenmat3x3f_multTranspose(benchmark::State& state)
+template<typename ScalarType>
+void BM_Matrix_eigenmat3x3_multTranspose(benchmark::State& state)
 {
     constexpr auto totalsize = maxSubIterations * 3*3 * 2;
-    const std::array<float, totalsize>& values = RandomValuePool<float, totalsize>::get();
+    const std::array<ScalarType, totalsize>& values = RandomValuePool<ScalarType, totalsize>::get();
 
-    std::vector<Eigen::Matrix<float, 3, 3> > vc1;
-    std::vector<Eigen::Matrix<float, 3, 3> > vc2;
+    using Matrix = Eigen::Matrix<ScalarType, 3, 3>;
+    
+    std::vector<Matrix> vc1;
+    std::vector<Matrix> vc2;
     vc1.reserve(state.range(0));
     vc2.reserve(state.range(0));
 
@@ -517,7 +619,7 @@ void BM_Matrix_eigenmat3x3f_multTranspose(benchmark::State& state)
 
     for (unsigned int i = 0; i < state.range(0); i++)
     {
-        Eigen::Matrix<float, 3, 3> mat1, mat2;
+        Matrix mat1, mat2;
         for (int a = 0; a < 3; ++a)
         {
             for (int b = 0; b < 3; ++b)
@@ -534,8 +636,8 @@ void BM_Matrix_eigenmat3x3f_multTranspose(benchmark::State& state)
     {
         for (auto i = 0 ; i < state.range(0); i++)
         {
-            Eigen::Matrix<float, 3, 3 > product = vc1[i].transpose() * vc2[i];
-            benchmark::DoNotOptimize(product);
+            Matrix res = vc1[i].transpose() * vc2[i];
+            benchmark::DoNotOptimize(res);
         }
     }
 }