Merge pull request #16 from SpM-lab/bugfix20241106

Work around an issue in Eigen3

include/sparseir/_linalg.hpp

@@ -199,23 +199,6 @@ xprec::DDouble _copysign(xprec::DDouble x, xprec::DDouble y) {
     return xprec::copysign(x, y);
 }
 
-double _sqrt(double x) {
-    return std::sqrt(x);
-}
-
-xprec::DDouble _sqrt(xprec::DDouble x) {
-    return xprec::sqrt(x);
-}
-
-double _abs(double x) {
-    return std::abs(x);
-}
-
-xprec::DDouble _abs(xprec::DDouble x) {
-    return xprec::abs(x);
-}
-
-
 /*
 This implementation is based on Julia's LinearAlgebra.refrector! function.
 
@@ -285,6 +268,9 @@ void reflectorApply(Eigen::VectorBlock<Eigen::Block<Eigen::MatrixX<T>, -1, 1, tr
 
 template <typename T>
 std::pair<QRPivoted<T>, int> rrqr(MatrixX<T>& A, T rtol = std::numeric_limits<T>::epsilon()) {
+    using std::abs;
+    using std::sqrt;
+
     int m = A.rows();
     int n = A.cols();
     int k = std::min(m, n);
@@ -294,7 +280,7 @@ std::pair<QRPivoted<T>, int> rrqr(MatrixX<T>& A, T rtol = std::numeric_limits<T>
 
     Vector<T, Dynamic> xnorms = A.colwise().norm();
     Vector<T, Dynamic> pnorms = xnorms;
-    T sqrteps = _sqrt(std::numeric_limits<T>::epsilon());
+    T sqrteps = sqrt(std::numeric_limits<T>::epsilon());
 
     for (int i = 0; i < k; ++i) {
 
@@ -315,24 +301,21 @@ std::pair<QRPivoted<T>, int> rrqr(MatrixX<T>& A, T rtol = std::numeric_limits<T>
         reflectorApply<T>(Ainp, tau_i, block);
 
         for (int j = i + 1; j < n; ++j) {
-            T temp = _abs((A(i, j))) / pnorms[j];
+            T temp = abs((A(i, j))) / pnorms[j];
             temp = std::max<T>(T(0), (T(1) + temp) * (T(1) - temp));
             // 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;
             } else {
-                pnorms(j) = pnorms(j) * _sqrt(temp);
+                pnorms(j) = pnorms(j) * sqrt(temp);
             }
         }
 
-        // std::cout << "rtol" << rtol << std::endl;
-        // std::cout << _abs(A(i,i)) << " v.s " << rtol * _abs((A(0,0))) << std::endl;
-        // std::cout << (_abs(A(i,i)) < rtol * _abs((A(0,0)))) << std::endl;
-
-        if (_abs(A(i,i)) < rtol * _abs((A(0,0)))) {
+        if (abs(A(i,i)) < rtol * abs((A(0,0)))) {
             A.bottomRightCorner(m - i, n - i).setZero();
             taus.tail(k - i).setZero();
             rk = i;

test/_linalg.cxx

@@ -236,10 +236,13 @@ TEST_CASE("RRQR Trunc", "[linalg]") {
         }
         // DDouble
         {
+                using std::pow;
                 VectorX<DDouble> x = VectorX<DDouble>::LinSpaced(101, -1, 1);
                 MatrixX<DDouble> Aorig(101, 21);
-                for (int i = 0; i < 21; i++) {
-                        Aorig.col(i) = x.array().pow(i);
+                for (int i = 0; i < Aorig.cols(); i++) {
+                        for (int j = 0; j < Aorig.rows(); j++) {
+                                Aorig(j, i) = pow(x(j), i); // xprec::pow is called.
+                        }
                 }
 
                 MatrixX<DDouble> A = Aorig;
@@ -248,17 +251,14 @@ TEST_CASE("RRQR Trunc", "[linalg]") {
                 QRPivoted<DDouble> A_qr;
                 int k;
                 std::tie(A_qr, k) = rrqr<DDouble>(A, DDouble(0.00001));
-                /*
                 REQUIRE(k < std::min(m, n));
-                // This fails
                 REQUIRE(k == 17);
                 auto QR = truncate_qr_result<DDouble>(A_qr, k);
                 auto Q = QR.first;
                 auto R = QR.second;
 
                 MatrixX<DDouble> A_rec = Q * R * getPropertyP(A_qr).transpose();
                 REQUIRE(A_rec.isApprox(Aorig, 1e-5 * A.norm()));
-                */
         }
 }