diff --git a/.gitignore b/.gitignore
index 6976e6147..75e5d999c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,4 @@ test/Testing/Temporary/CTestCostData.txt
 build/
 .vscode
 .DS_Store
+.cache
\ No newline at end of file
diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h
index bfd506a53..e30aa25bb 100644
--- a/include/convex_bodies/hpolytope.h
+++ b/include/convex_bodies/hpolytope.h
@@ -96,7 +96,7 @@ class HPolytope {
         for (unsigned int i = 1; i < Pin.size(); i++) {
             b(i - 1) = Pin[i][0];
             for (unsigned int j = 1; j < _d + 1; j++) {
-                A(i - 1, j - 1) = -Pin[i][j];
+                A.coeffRef(i - 1, j - 1) = -Pin[i][j];
             }
         }
         has_ball = false;
@@ -645,12 +645,12 @@ class HPolytope {
 
         int m = num_of_hyperplanes();
 
-        lamdas.noalias() += A.col(rand_coord_prev)
-                         * (r_prev[rand_coord_prev] - r[rand_coord_prev]);
+        lamdas.noalias() += (DenseMT)(A.col(rand_coord_prev)
+                         * (r_prev[rand_coord_prev] - r[rand_coord_prev]));
         NT* data = lamdas.data();
 
         for (int i = 0; i < m; i++) {
-            NT a = A(i, rand_coord);
+            NT a = A.coeff(i, rand_coord);
 
             if (a == NT(0)) {
                 //std::cout<<"div0"<<std::endl;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 61e499a9e..587966638 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -263,6 +263,7 @@ add_test(NAME volume_cg_hpolytope_birkhoff COMMAND volume_cg_hpolytope -tc=birk)
 add_test(NAME volume_cg_hpolytope_prod_simplex COMMAND volume_cg_hpolytope -tc=prod_simplex)
 add_test(NAME volume_cg_hpolytope_simplex COMMAND volume_cg_hpolytope -tc=simplex)
 add_test(NAME volume_cg_hpolytope_skinny_cube COMMAND volume_cg_hpolytope -tc=skinny_cube)
+add_test(NAME volume_cg_hpolytope_sparse_simplex COMMAND volume_cg_hpolytope -tc=sparse_simplex)
 
 add_executable (volume_cg_vpolytope volume_cg_vpolytope.cpp $<TARGET_OBJECTS:test_main>)
 add_test(NAME volume_cg_vpolytope_cube COMMAND volume_cg_vpolytope -tc=cube)
@@ -305,6 +306,8 @@ add_test(NAME round_volumetric_barrier_test
           COMMAND rounding_test -tc=round_volumetric_barrier_test)
 add_test(NAME round_vaidya_barrier_test
           COMMAND rounding_test -tc=round_vaidya_barrier_test)
+add_test(NAME round_sparse_test
+          COMMAND rounding_test -tc=round_sparse)
 
 
 
diff --git a/test/rounding_test.cpp b/test/rounding_test.cpp
index 5beae0415..b4739b47e 100644
--- a/test/rounding_test.cpp
+++ b/test/rounding_test.cpp
@@ -57,7 +57,7 @@ void rounding_min_ellipsoid_test(Polytope &HP,
 {
     typedef typename Polytope::PointType Point;
     typedef typename Point::FT NT;
-    typedef typename Polytope::MT MT;
+    typedef typename Polytope::DenseMT MT;
     typedef typename Polytope::VT VT;
 
     int d = HP.dimension();
@@ -102,7 +102,7 @@ void rounding_max_ellipsoid_test(Polytope &HP,
 {
     typedef typename Polytope::PointType Point;
     typedef typename Point::FT NT;
-    typedef typename Polytope::MT MT;
+    typedef typename Polytope::DenseMT MT;
     typedef typename Polytope::VT VT;
 
     int d = HP.dimension();
@@ -178,7 +178,7 @@ void rounding_log_barrier_test(Polytope &HP,
 {
     typedef typename Polytope::PointType Point;
     typedef typename Point::FT NT;
-    typedef typename Polytope::MT MT;
+    typedef typename Polytope::DenseMT MT;
     typedef typename Polytope::VT VT;
 
     int d = HP.dimension();
@@ -208,7 +208,7 @@ void rounding_volumetric_barrier_test(Polytope &HP,
 {
     typedef typename Polytope::PointType Point;
     typedef typename Point::FT NT;
-    typedef typename Polytope::MT MT;
+    typedef typename Polytope::DenseMT MT;
     typedef typename Polytope::VT VT;
 
     int d = HP.dimension();
@@ -238,7 +238,7 @@ void rounding_vaidya_barrier_test(Polytope &HP,
 {
     typedef typename Polytope::PointType Point;
     typedef typename Point::FT NT;
-    typedef typename Polytope::MT MT;
+    typedef typename Polytope::DenseMT MT;
     typedef typename Polytope::VT VT;
 
     int d = HP.dimension();
@@ -381,6 +381,39 @@ void call_test_svd() {
     rounding_svd_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0);
 }
 
+template <typename NT>
+void call_test_sparse() {
+    typedef Cartesian <NT> Kernel;
+    typedef typename Kernel::Point Point;
+    typedef HPolytope <Point, Eigen::SparseMatrix<NT>> Hpolytope;
+    Hpolytope P;
+
+    //std::cout << "\n--- Testing SVD rounding of sparse H-skinny_cube5" << std::endl;
+    //P = generate_skinny_cube<Hpolytope>(5);
+    //rounding_svd_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0);
+
+    std::cout << "\n--- Testing log-barrier rounding of sparse H-skinny_cube5" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(5);
+    rounding_log_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0);
+
+    std::cout << "\n--- Testing volumetric barrier rounding of sparse H-skinny_cube5" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(5);
+    rounding_volumetric_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0);
+
+    std::cout << "\n--- Testing vaidya barrier rounding of sparse H-skinny_cube5" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(5);
+    rounding_vaidya_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0);
+
+    std::cout << "\n--- Testing min ellipsoid rounding of sparse H-skinny_cube10" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(10);
+    rounding_min_ellipsoid_test(P, 0, 122550, 108426, 105003.0, 102400.0);
+
+    std::cout << "\n--- Testing max ellipsoid rounding of sparse H-skinny_cube5" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(5);
+    rounding_max_ellipsoid_test(P, 0, 3070.64, 3188.25, 3262.61, 3200.0);
+}
+
+
 
 TEST_CASE("round_min_ellipsoid") {
     call_test_min_ellipsoid<double>();
@@ -410,3 +443,6 @@ TEST_CASE("round_svd") {
     call_test_svd<double>();
 }
 
+TEST_CASE("round_sparse") {
+    call_test_sparse<double>();
+}
diff --git a/test/volume_cg_hpolytope.cpp b/test/volume_cg_hpolytope.cpp
index 2e9060e45..e924c4ba3 100644
--- a/test/volume_cg_hpolytope.cpp
+++ b/test/volume_cg_hpolytope.cpp
@@ -251,6 +251,32 @@ void call_test_skinny_cube() {
     //test_volume(P, 104857600, 104857600.0);
 }
 
+template <typename NT>
+void call_test_sparse_simplex() {
+    typedef Cartesian<NT>    Kernel;
+    typedef typename Kernel::Point    Point;
+
+    typedef HPolytope<Point, Eigen::SparseMatrix<NT>> Hpolytope;
+    Hpolytope SP;
+
+    std::cout << "--- Testing volume of sparse H-simplex10" << std::endl;
+    SP = generate_simplex<Hpolytope>(10, false);
+    test_volume(SP,
+                2.14048 * std::pow(10,-7),
+                2.70598 * std::pow(10,-7),
+                2.53893e-07,
+                1.0 / factorial(10.0));
+
+    std::cout << "--- Testing volume of sparse H-simplex20" << std::endl;
+    SP = generate_simplex<Hpolytope>(20, false);
+    test_volume(SP,
+                2.00646 * std::pow(10,-21),
+                4.16845 * std::pow(10,-19),
+                5.0348e-19,
+                1.0 / factorial(20.0));
+}
+
+
 
 TEST_CASE("cube") {
     call_test_cube<double>();
@@ -276,3 +302,7 @@ TEST_CASE("simplex") {
 TEST_CASE("skinny_cube") {
     call_test_skinny_cube<double>();
 }
+
+TEST_CASE("sparse_simplex") {
+    call_test_sparse_simplex<double>();
+}