diff --git a/include/sparseir/_linalg.hpp b/include/sparseir/_linalg.hpp
index 138e52c..53d1f07 100644
--- a/include/sparseir/_linalg.hpp
+++ b/include/sparseir/_linalg.hpp
@@ -306,7 +306,6 @@ std::pair<QRPivoted<T>, int> rrqr(MatrixX<T>& A, T rtol = std::numeric_limits<T>
             // abs2
             T temp2 = temp * (pnorms(j) / xnorms(j)) * (pnorms(j) / xnorms(j));
             if (temp2 < sqrteps) {
-                std::cout << "Called" << std::endl;
                 auto recomputed = A.col(j).tail(m - i - 1).norm();
                 pnorms(j) = recomputed;
                 xnorms(j) = recomputed;
@@ -506,7 +505,13 @@ tsvd(Eigen::MatrixX<T>& A, T rtol = std::numeric_limits<T>::epsilon()) {
 
     // Step 3: Compute SVD of R_trunc
 
-    Eigen::JacobiSVD<Eigen::MatrixX<T>> svd(R_trunc.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV);
+    //Eigen::JacobiSVD<Eigen::MatrixX<T>> svd(R_trunc.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV);
+    // There seems to be a bug in the latest version of Eigen3.
+    // Please first construct a Jacobi SVD and then compare the results.
+    // Do not use the svd_jacobi function directly.
+    // Better to write a wrrapper function for the SVD.
+    Eigen::JacobiSVD<decltype(R_trunc)> svd;
+    svd.compute(R_trunc.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV);
 
     Eigen::PermutationMatrix<Dynamic, Dynamic> perm(p.size());
     perm.indices() = invperm(p);
diff --git a/test/_linalg.cxx b/test/_linalg.cxx
index 60b27e4..481605b 100644
--- a/test/_linalg.cxx
+++ b/test/_linalg.cxx
@@ -263,13 +263,17 @@ TEST_CASE("RRQR Trunc", "[linalg]") {
 }
 
 TEST_CASE("TSVD", "[linalg]") {
+        using std::pow;
         // double
         {
             for (auto tol : {1e-14, 1e-13}) {
                 VectorX<double> x = VectorX<double>::LinSpaced(201, -1, 1);
                 MatrixX<double> Aorig(201, 51);
-                for (int i = 0; i < 51; i++) {
-                    Aorig.col(i) = x.array().pow(i);
+                for (int i = 0; i < Aorig.cols(); i++) {
+                    //Aorig.col(i) = x.array().pow(i);
+                    for (int j = 0; j < Aorig.rows(); j++) {
+                        Aorig(j, i) = pow(x(j), i);
+                    }
                 }
 
                 MatrixX<double> A = Aorig;
@@ -285,8 +289,9 @@ TEST_CASE("TSVD", "[linalg]") {
                 // std::cout << "S: " << S_diag.rows() << "," << S_diag.cols() << std::endl;
                 // std::cout << "V: " << V.rows() << "," << V.cols() << std::endl;
                 auto B = U * S_diag * V.transpose() - Aorig;
-                //std::cout << B << std::endl; // Oh...?
-                REQUIRE(1==1);
+                std::cout << Aorig.norm() << std::endl; // Oh...?
+                std::cout << B.norm() << std::endl; // Oh...?
+                REQUIRE(0==1);
                 //REQUIRE((U * S_diag * V).isApprox(Aorig, tol * A.norm()));
                 /*
                 REQUIRE((U.transpose() * U).isIdentity());