diff --git a/test/SparseGenMatProd.cpp b/test/SparseGenMatProd.cpp
index 332a5bdd..1941dd26 100644
--- a/test/SparseGenMatProd.cpp
+++ b/test/SparseGenMatProd.cpp
@@ -36,7 +36,31 @@ Eigen::SparseMatrix<T> generate_random_sparse(Index rows, Index cols)
     return mat;
 }
 
-TEMPLATE_TEST_CASE("matrix operations", "[SparseGenMatProd]", float, double)
+template <typename T>
+Eigen::SparseMatrix<T> generate_random_sparse_entries(Index rows, Index cols, int entries)
+{
+    std::default_random_engine gen;
+    std::uniform_real_distribution<T> dist(0.0, 1.0);
+    std::uniform_real_distribution<double> dist_x(0, rows);
+    std::uniform_real_distribution<double> dist_y(0, cols);
+
+    std::vector<Eigen::Triplet<T>> tripletVector;
+    for (int n = 0; n < entries ; n++) {
+        auto v_ij = dist(gen);
+        int x_c = std::round(dist_x(gen));
+        int y_c = std::round(dist_y(gen));
+        //Dont care about duplicates
+        tripletVector.push_back(Eigen::Triplet<T>(x_c, y_c, v_ij));
+    }
+    
+    Eigen::SparseMatrix<T> mat(rows, cols);
+    //create the matrix
+    mat.setFromTriplets(tripletVector.begin(), tripletVector.end());
+
+    return mat;
+}
+
+TEMPLATE_TEST_CASE("matrix operations [100x100]", "[SparseGenMatProd]", float, double)
 {
     std::srand(123);
     constexpr Index n = 100;
@@ -53,3 +77,26 @@ TEMPLATE_TEST_CASE("matrix operations", "[SparseGenMatProd]", float, double)
     INFO("The accesor operator must produce the same element as in eigen")
     REQUIRE(mat1.coeff(45, 22) == sparse1(45, 22));
 }
+
+TEMPLATE_TEST_CASE("matrix operations [100000000x100000000]", "[SparseGenMatProd]", float, double)
+{
+    std::srand(123);
+    constexpr Index n = 100000000; //adding another 0 will explode it, so perhaps it can be more efficient
+
+    Eigen::SparseMatrix<TestType> mat1 = generate_random_sparse_entries<TestType>(n, n, 1000);
+    Eigen::SparseMatrix<TestType> mat2 = generate_random_sparse_entries<TestType>(n, n, 1000);
+
+    SparseGenMatProd<TestType> sparse1(mat1);
+    INFO("The matrix-matrix product must not explode.")
+    Eigen::SparseMatrix<TestType> xs = sparse1 * mat2;
+    if (xs.coeff(50000,3042) == 0) {
+        xs.coeffRef(50000,3242) = 1; // Dont let eigen optimalize it away
+    } else
+    {
+        xs.coeffRef(50000,3242) = 2; // Dont let eigen optimalize it away
+    }
+    
+
+    INFO("The accessor operator must not explode")
+    REQUIRE(mat1.coeff(45, 22) == sparse1(45, 22));
+}
diff --git a/test/SparseSymMatProd.cpp b/test/SparseSymMatProd.cpp
index 0b3d36c9..5425ac3c 100644
--- a/test/SparseSymMatProd.cpp
+++ b/test/SparseSymMatProd.cpp
@@ -34,6 +34,29 @@ Eigen::SparseMatrix<T> generate_random_sparse(Index rows, Index cols)
     return mat;
 }
 
+template <typename T>
+Eigen::SparseMatrix<T> generate_random_sparse_entries(Index rows, Index cols, int entries)
+{
+    std::default_random_engine gen;
+    std::uniform_real_distribution<T> dist(0.0, 1.0);
+    std::uniform_real_distribution<double> dist_x(0, rows);
+    std::uniform_real_distribution<double> dist_y(0, cols);
+
+    std::vector<Eigen::Triplet<T>> tripletVector;
+    for (int n = 0; n < entries ; n++) {
+        auto v_ij = dist(gen);
+        int x_c = std::round(dist_x(gen));
+        int y_c = std::round(dist_y(gen));
+        //Dont care about duplicates
+        tripletVector.push_back(Eigen::Triplet<T>(x_c, y_c, v_ij));
+    }
+    
+    Eigen::SparseMatrix<T> mat(rows, cols);
+    //create the matrix
+    mat.setFromTriplets(tripletVector.begin(), tripletVector.end());
+
+    return mat;
+}
 TEMPLATE_TEST_CASE("matrix operations", "[SparseSymMatProd]", float, double)
 {
     std::srand(123);
@@ -52,3 +75,27 @@ TEMPLATE_TEST_CASE("matrix operations", "[SparseSymMatProd]", float, double)
     INFO("The accesor operator must produce the same element as in eigen")
     REQUIRE(mat1.coeff(45, 22) == sparse1(45, 22));
 }
+
+TEMPLATE_TEST_CASE("matrix operations [100000000x100000000]", "[SparseSymMatProd]", float, double)
+{
+    std::srand(123);
+    constexpr Index n = 100000000; //adding another 0 will explode it, so perhaps it can be more efficient
+
+    Eigen::SparseMatrix<TestType> mat = generate_random_sparse_entries<TestType>(n, n, 500);
+    Eigen::SparseMatrix<TestType> mat1 = mat + Eigen::SparseMatrix<TestType>(mat.transpose());  // It needs to be symetric
+    Eigen::SparseMatrix<TestType> mat2 = generate_random_sparse_entries<TestType>(n, n, 1000);
+
+    SparseSymMatProd<TestType> sparse1(mat1);
+    INFO("The matrix-matrix product must not explode.")
+    Eigen::SparseMatrix<TestType> xs = sparse1 * mat2;
+    if (xs.coeff(50000,3042) == 0) {
+        xs.coeffRef(50000,3242) = 1; // Dont let eigen optimalize it away
+    } else
+    {
+        xs.coeffRef(50000,3242) = 2; // Dont let eigen optimalize it away
+    }
+    
+
+    INFO("The accessor operator must not explode")
+    REQUIRE(mat1.coeff(45, 22) == sparse1(45, 22));
+}