diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp
index 0655767ea..9faaf9553 100644
--- a/examples/correlation_matrices/sampler.cpp
+++ b/examples/correlation_matrices/sampler.cpp
@@ -90,28 +90,6 @@ void write_to_file(std::string filename, std::vector<PointType> const& randPoint
     std::cout.rdbuf(coutbuf);
 }
 
-bool is_correlation_matrix(const MT& matrix, const double tol = 1e-8){
-    //check if all the diagonal elements are ones
-    for(int i=0 ; i<matrix.rows() ; i++)
-    {
-        if(std::abs(matrix(i, i)-1.0) > tol)
-	{
-	    return false;
-	}
-    }
-
-    //check if the matrix is positive semidefinite
-    using NT = double;
-    using MatrixType = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>;
-    EigenvaluesProblems<NT, MatrixType, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;
-
-    if(solver.isPositiveSemidefinite(matrix))
-    {
-        return true;
-    }
-    return false;
-}
-
 template<typename WalkType>
 void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){
 
@@ -138,12 +116,12 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned
     std::cout << walkname << " samples uniformly "<< num_points << " correlation matrices of size " << n << " with matrix PointType" << std::endl;
     std::chrono::steady_clock::time_point start, end;
     double time;
-    std::vector<PointMT> randPoints;
+    std::list<MT> randCorMatrices;
     unsigned int walkL = 1;
 
     start = std::chrono::steady_clock::now();
 
-    uniform_correlation_sampling_MT<WalkType, PointMT, RNGType>(n, randPoints, walkL, num_points, 0);
+    uniform_correlation_sampling_MT<WalkType, PointMT, RNGType>(n, randCorMatrices, walkL, num_points, 0);
 
     end = std::chrono::steady_clock::now();
     time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
@@ -151,13 +129,21 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned
 
     int valid_points = 0;
     EigenvaluesProblems<NT, MT, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;
-    for(const auto& points : randPoints){
-        if(solver.is_correlation_matrix(points.mat)){
+    for(const auto& matrix : randCorMatrices){
+    	if(solver.is_correlation_matrix(matrix)){
             valid_points++;
-        }
+    	}
     }
+
     std::cout << "Number of valid points = " << valid_points << std::endl;
 
+    std::vector<PointMT> randPoints;
+    for(const auto &mat : randCorMatrices){
+        PointMT p;
+    	p.mat = mat;
+    	randPoints.push_back(p);
+    }
+	
     write_to_file<PointMT>(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints);
 }
 
diff --git a/include/sampling/sample_correlation_matrices.hpp b/include/sampling/sample_correlation_matrices.hpp
index 9a1de2319..50b7a45da 100644
--- a/include/sampling/sample_correlation_matrices.hpp
+++ b/include/sampling/sample_correlation_matrices.hpp
@@ -41,19 +41,27 @@ template
     typename WalkTypePolicy,
     typename PointType,
     typename RNGType,
-    typename PointList
+    typename MT
 >
 void uniform_correlation_sampling_MT(   const unsigned int &n,
-                                        PointList &randPoints,
+                                        std::list<MT> &randCorMatrices,
                                         const unsigned int &walkL,
                                         const unsigned int &num_points,
                                         unsigned int const& nburns){
+    using PointList = std::list<PointType>;
+    PointList randPoints;
+    
     CorrelationSpectrahedron_MT<PointType> P(n);
     const unsigned int d = P.dimension();
     PointType startingPoint(n);
     RNGType rng(d);
 
     uniform_sampling<WalkTypePolicy>(randPoints, P, rng, walkL, num_points, startingPoint, nburns);
+
+    for(const auto&p : randPoints){
+        MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n);
+    	randCorMatrices.push_back(final_cor_mat);
+    }
 }
 
 template
@@ -83,21 +91,29 @@ template
     typename WalkTypePolicy,
     typename PointType,
     typename RNGType,
-    typename PointList,
+    typename MT,
     typename NT
 >
 void gaussian_correlation_sampling_MT(  const unsigned int &n,
-                                        PointList &randPoints,
+                                        std::list<MT> &randCorMatrices,
                                         const unsigned int &walkL,
                                         const unsigned int &num_points,
                                         const NT &a,
                                         unsigned int const& nburns = 0){
+    using PointList = std::list<PointType>;
+    PointList randPoints;
+    
     CorrelationSpectrahedron_MT<PointType> P(n);
     const unsigned int d = P.dimension();
     PointType startingPoint(n);
     RNGType rng(d);
 
     gaussian_sampling<WalkTypePolicy>(randPoints, P, rng, walkL, num_points, a, startingPoint, nburns);
+
+    for(const auto&p : randPoints){
+        MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n);
+        randCorMatrices.push_back(final_cor_mat);
+    }
 }
 
 template
@@ -130,17 +146,20 @@ template
         typename WalkTypePolicy,
         typename PointType,
         typename RNGType,
-        typename PointList,
+        typename MT,
         typename NT,
         typename VT
 >
 void exponential_correlation_sampling_MT(   const unsigned int &n,
-                                            PointList &randPoints,
+                                            std::list<MT> &randCorMatrices,
                                             const unsigned int &walkL,
                                             const unsigned int &num_points,
                                             const VT &c,
                                             const NT &T,
                                             unsigned int const& nburns = 0){
+    using PointList = std::list<PointType>;
+    PointList randPoints;
+    
     CorrelationSpectrahedron_MT<PointType> P(n);
     const unsigned int d = P.dimension();
     PointType startingPoint(n);
@@ -148,6 +167,11 @@ void exponential_correlation_sampling_MT(   const unsigned int &n,
     PointType _c(c);
 
     exponential_sampling<WalkTypePolicy>(randPoints, P, rng, walkL, num_points, _c, T, startingPoint, nburns);
+
+    for(const auto&p : randPoints){
+        MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n);
+        randCorMatrices.push_back(final_cor_mat);
+    }
 }
 
 #endif //VOLESTI_SAMPLING_SAMPLE_CORRELATION_MATRICES_HPP
diff --git a/test/sampling_correlation_matrices_test.cpp b/test/sampling_correlation_matrices_test.cpp
index 39bb0e605..a17c36ce6 100644
--- a/test/sampling_correlation_matrices_test.cpp
+++ b/test/sampling_correlation_matrices_test.cpp
@@ -36,6 +36,21 @@ MT rebuildMatrix(const VT &xvector, const unsigned int n){
     return mat;
 }
 
+template<typename NT, typename MT>
+Eigen::Matrix<NT, Eigen::Dynamic, 1> getCoefficientsFromMatrix(const MT& mat) {
+	int n = mat.rows();
+	int d = n * (n - 1) / 2;
+	Eigen::Matrix<NT, Eigen::Dynamic, 1> coeffs(d);
+	int k = 0;
+	for (int i = 0; i < n; ++i) {
+    	for (int j = 0; j < i; ++j) {
+        	coeffs(k) = mat(i, j);
+        	++k;
+    	}
+	}
+	return coeffs;
+}
+
 template<typename NT, typename VT, typename MT, typename PointList>
 void check_output(PointList &randPoints, int num_points, int n){
     int d = n*(n-1)/2, count = 0;
@@ -64,6 +79,33 @@ void check_output(PointList &randPoints, int num_points, int n){
     CHECK(score.maxCoeff() < 1.1);
 }
 
+template<typename NT, typename VT, typename MT>
+void check_output_MT(std::list<MT> &randCorMatrices, int num_points, int n){
+    int d = n*(n-1)/2, count = 0;
+    MT A;
+    Eigen::LDLT<MT> mat_ldlt;
+    for(auto& mat : randCorMatrices){
+        mat_ldlt = Eigen::LDLT<MT>(mat);
+    	if(mat_ldlt.info() == Eigen::NumericalIssue || !mat_ldlt.isPositive()){
+            ++count;
+    	}
+    }
+    std::cout << "Fails " << count << " / " << num_points << " samples\n";
+    CHECK(count == 0);
+
+    MT samples(d, num_points);
+    unsigned int jj = 0;
+    for(const auto& mat : randCorMatrices){
+        samples.col(jj) = getCoefficientsFromMatrix<NT, MT>(mat);
+    	jj++;
+    }
+
+    VT score = univariate_psrf<NT, VT>(samples);
+    std::cout << "psrf = " << score.maxCoeff() << std::endl;
+
+    CHECK(score.maxCoeff() < 1.1);
+}
+
 template <typename NT>
 void test_corre_spectra_classes(unsigned int const n){
     typedef Cartesian<NT>                                       Kernel;
@@ -136,18 +178,18 @@ void test_new_uniform_MT(const unsigned int n, const unsigned int num_points = 1
     std::cout << "Test new sampling 2 : "<< num_points << " uniform correlation matrices of size " << n << std::endl;
     std::chrono::steady_clock::time_point start, end;
     double time;
-    std::vector<Point> randPoints;
+    std::list<MT> randCorMatrices;
     unsigned int walkL = 1;
 
     start = std::chrono::steady_clock::now();
 
-    uniform_correlation_sampling_MT<WalkType, Point, RNGType>(n, randPoints, walkL, num_points, 0);
+    uniform_correlation_sampling_MT<WalkType, Point, RNGType>(n, randCorMatrices, walkL, num_points, 0);
 
     end = std::chrono::steady_clock::now();
     time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
     std::cout << "Elapsed time : " << time << " (ms)" << std::endl;
 
-    check_output<NT, VT, MT>(randPoints, num_points, n);
+    check_output_MT<NT, VT, MT>(randCorMatrices, num_points, n);
 }
 
 int n = 3;