Simplify realization of the wrappers with Eigen array methods.

This should allow Eigen to do its internal optimizations based on storage order. We move the previous manual-loop implementation over to the tests so that we can confirm that we used the array methods correctly. Also fixed some of the wrapper tests from previous Workspace changes.
LTLA · Jul 3, 2024 · 78e6a6c · 78e6a6c
1 parent ab12fda
commit 78e6a6c
Show file tree

Hide file tree

Showing 2 changed files with 25 additions and 41 deletions.
diff --git a/include/irlba/wrappers.hpp b/include/irlba/wrappers.hpp
@@ -280,12 +280,7 @@ class Centered {
     template<class EigenMatrix_>
     Eigen::MatrixXd realize() const {
         auto output = wrapped_realize<EigenMatrix_>(my_matrix);
-        auto nc = output.cols(), nr = output.rows();
-        for (Eigen::Index c = 0; c < nc; ++c) {
-            for (Eigen::Index r = 0; r < nr; ++r) {
-                output(r, c) -= my_center[c];
-            }
-        }
+        output.array().rowwise() -= my_center.adjoint().array();
         return output;
     }
     /**
@@ -407,35 +402,18 @@ class Scaled {
     template<class EigenMatrix_>
     EigenMatrix_ realize() const {
         auto output = wrapped_realize<EigenMatrix_>(my_matrix);
-        auto nc = output.cols(), nr = output.rows();
 
         if constexpr(column_) {
             if (my_divide) {
-                for (Eigen::Index c = 0; c < nc; ++c) {
-                    for (Eigen::Index r = 0; r < nr; ++r) {
-                        output(r, c) /= my_scale[c];
-                    }
-                }
+                output.array().rowwise() /= my_scale.adjoint().array();
             } else {
-                for (Eigen::Index c = 0; c < nc; ++c) {
-                    for (Eigen::Index r = 0; r < nr; ++r) {
-                        output(r, c) *= my_scale[c];
-                    }
-                }
+                output.array().rowwise() *= my_scale.adjoint().array();
             }
         } else {
             if (my_divide) {
-                for (Eigen::Index c = 0; c < nc; ++c) {
-                    for (Eigen::Index r = 0; r < nr; ++r) {
-                        output(r, c) /= my_scale[r];
-                    }
-                }
+                output.array().colwise() /= my_scale.array();
             } else {
-                for (Eigen::Index c = 0; c < nc; ++c) {
-                    for (Eigen::Index r = 0; r < nr; ++r) {
-                        output(r, c) *= my_scale[r];
-                    }
-                }
+                output.array().colwise() *= my_scale.array();
             }
         }
 

diff --git a/tests/src/wrappers.cpp b/tests/src/wrappers.cpp
@@ -22,8 +22,10 @@ TEST_F(WrapperCenteringTest, Basic) {
 
     auto realized = centered.template realize<Eigen::MatrixXd>();
     Eigen::MatrixXd expected = A;
-    for (Eigen::Index r = 0; r < A.rows(); ++r) {
-        expected.row(r).array() -= B.array();
+    for (Eigen::Index c = 0; c < A.cols(); ++c) {
+        for (Eigen::Index r = 0; r < A.rows(); ++r) {
+            expected(r, c) -= B[c];
+        }
     }
     expect_equal_matrix(expected, realized);
 }
@@ -37,7 +39,7 @@ TEST_F(WrapperCenteringTest, Multiply) {
 
     auto wrk = centered.workspace();
     bool is_placeholder = std::is_same<decltype(wrk), bool>::value;
-    EXPECT_TRUE(is_placeholder); // just inherits the child.
+    EXPECT_FALSE(is_placeholder);
 
     Eigen::VectorXd output(20);
     centered.multiply(C, wrk, output);
@@ -58,7 +60,7 @@ TEST_F(WrapperCenteringTest, AdjointMultiply) {
 
     auto wrk = centered.adjoint_workspace();
     bool is_placeholder = std::is_same<decltype(wrk), bool>::value;
-    EXPECT_TRUE(is_placeholder); // just inherits the child.
+    EXPECT_FALSE(is_placeholder);
 
     Eigen::VectorXd output(10);
     centered.adjoint_multiply(C, wrk, output);
@@ -94,11 +96,13 @@ TEST_P(WrapperScalingTest, Basic) {
 
         auto realized = scaled.template realize<Eigen::MatrixXd>();
         Eigen::MatrixXd expected = A;
-        for (Eigen::Index r = 0; r < expected.rows(); ++r) {
-            if (divide) {
-                expected.row(r).array() /= B1.array();
-            } else {
-                expected.row(r).array() *= B1.array();
+        for (Eigen::Index c = 0; c < A.cols(); ++c) {
+            for (Eigen::Index r = 0; r < A.rows(); ++r) {
+                if (divide) {
+                    expected(r, c) /= B1[c];
+                } else {
+                    expected(r, c) *= B1[c];
+                }
             }
         }
         expect_equal_matrix(expected, realized);
@@ -111,11 +115,13 @@ TEST_P(WrapperScalingTest, Basic) {
 
         auto realized = scaled.template realize<Eigen::MatrixXd>();
         Eigen::MatrixXd expected = A;
-        for (Eigen::Index c = 0; c < expected.cols(); ++c) {
-            if (divide) {
-                expected.col(c).array() /= B2.array();
-            } else {
-                expected.col(c).array() *= B2.array();
+        for (Eigen::Index c = 0; c < A.cols(); ++c) {
+            for (Eigen::Index r = 0; r < A.rows(); ++r) {
+                if (divide) {
+                    expected(r, c) /= B2[r];
+                } else {
+                    expected(r, c) *= B2[r];
+                }
             }
         }
         expect_equal_matrix(expected, realized);