[Bug]: Improve runtime performance with Eigen types

[chiphogg]

Bug Description

To be clear: this isn't about mp-units' design choice to materialize the rep on every operation. That was a conscious design decision. Instead, I'm talking about some "low hanging fruit": extra performance penalties that don't need to be there.

In the implementations for +, -, %, and comparison, we make copies instead of taking references for the cases where we don't need to make a conversion. This is fine for scalar types, but can result in a measurable performance regression with Eigen types.

Claude tells me that this patch file shows one way to fix it. I'd take it as an illustration rather than a complete solution.

mp-units-no-operand-copies.patch

Note that this patch might be old, and may or may not apply directly. 😅 I had filing this bug on my todo list for a little while.

Separately: the cross product implementation in test/runtime/linear_algebra_test.cc uses numerical_value_in, not numerical_value_ref_in, which always forces copies. This makes the cross product about 9x slower than raw Eigen. It's just an example rather than end user facing code, but still, it's low hanging fruit to improve.

Steps to Reproduce

I had Claude whip up this MRE:

#include <chrono>
#include <cstdio>
#include <vector>
#include <Eigen/Core>
#include <mp-units/integrations/eigen.h>
#include <mp-units/systems/isq/space_and_time.h>
#include <mp-units/systems/si.h>

using namespace mp_units;
using namespace mp_units::si::unit_symbols;
using Vec = Eigen::Matrix<double, 15, 1>;

template <class T, class F> double ns_per_op(const std::vector<T>& a, const std::vector<T>& b, F op) {
    double best = 1e30;
    for (int rep = 0; rep < 5; ++rep) {
        auto t0 = std::chrono::steady_clock::now();
        for (int it = 0; it < 20000; ++it)
            for (size_t i = 0; i < a.size(); ++i) { auto r = op(a[i], b[i]); asm volatile("" : : "g"(&r) : "memory"); }
        double ns = std::chrono::duration<double, std::nano>(std::chrono::steady_clock::now() - t0).count();
        best = std::min(best, ns / (20000.0 * a.size()));
    }
    return best;
}

int main() {
    std::vector<Vec> a(32), b(32);
    for (auto& v : a) v.setRandom();
    for (auto& v : b) v.setRandom();
    std::vector<Q> qa, qb;
    for (auto& v : a) qa.push_back(v * isq::displacement[m]);
    for (auto& v : b) qb.push_back(v * isq::displacement[m]);

    double raw = ns_per_op(a, b, [](const Vec& x, const Vec& y) { return Vec(x + y); });
    double mp  = ns_per_op(qa, qb, [](const Q& x, const Q& y) { return x + y; });
    double eq_raw = ns_per_op(a, a, [](const Vec& x, const Vec& y) { return x == y; });
    double eq_mp  = ns_per_op(qa, qa, [](const Q& x, const Q& y) { return x == y; });
    std::printf("a + b   : Eigen %.2f ns   mp-units %.2f ns   (%.2fx)\n", raw, mp, mp / raw);
    std::printf("a == b  : Eigen %.2f ns   mp-units %.2f ns   (%.2fx)\n", eq_raw, eq_mp, eq_mp / eq_raw);
}

Output from one run:

a + b   : Eigen 2.51 ns   mp-units 8.04 ns   (3.20x)
a == b  : Eigen 3.27 ns   mp-units 7.41 ns   (2.27x)

With the attached patch (operands read by reference when they already have the common type):

a + b   : Eigen 1.89 ns   mp-units 2.15 ns   (1.14x)
a == b  : Eigen 3.16 ns   mp-units 3.16 ns   (1.00x)

🚀 Compiler Explorer Link (Recommended)

No response

Compiler

GCC 12

C++ Standard

C++20

mp-units Version

e2d6b39

Build Configuration

• Using Conan
• Using CMake directly
• Header-only mode
• C++ modules enabled
• Freestanding mode
• Debug build
• Release build

Additional Context

No response

Checklist

• I have read the documentation and checked existing issues