How to execute mathematical operations between two boost::multi_arrays (C++)?

Question

How can I execute a mathematical operation between two boost::multi_arrays?

Example of adding two arrays with value type double:

auto array1 = boost::multi_array<double, 2>(boost::extents[10][10]);
auto array2 = boost::multi_array<double, 2>(boost::extents[10][10]);

auto array3 = array1  + array2; //Does not compile

One possibility I know is the Intel IPP library. Adding two matrices can be done with e.g. ippiAdd_. But Intel IPP does unfortunately not support double values for adding.

So does someone know another library than Intel IPP respectively a solution to overcome the shortcomings of restricted value types in Intel IPP?

Answer 1

You have to overload the + operator for those your types of objects boost::multi_array<double, 2> with your desired implementation.

EDIT I just tried real-quick, apparently it was not so hard, but maybe needs more testing and review ;)

Here you go:

boost::multi_array<double, 2> operator+(boost::multi_array<double, 2> a, boost::multi_array<double, 2> b) {
    boost::multi_array<double, 2> result = a;
    for (size_t i=0; i<a.size(); ++i) {
        for (size_t j=0; j<a[i].size(); ++j) {
            result[i][j] = a[i][j] + b[i][j];
        }
    }
    return result;
}

Answer 2

You could "just write it":

namespace ArrayOperators {
    template <typename L, typename R>
    static inline auto operator+(L const& l, R const& r) {
        return ArrayOp {std::plus<>{}} (l, r); }

    template <typename L, typename R>
    static inline auto operator-(L const& l, R const& r) {
        return ArrayOp {std::minus<>{}} (l, r); }

    template <typename L, typename R>
    static inline auto operator/(L const& l, R const& r) {
        return ArrayOp {std::divides<>{}} (l, r); }

    template <typename L, typename R>
    static inline auto operator*(L const& l, R const& r) {
        return ArrayOp {std::multiplies<>{}} (l, r); }
}

Of course, this requires us to actually implement the ArrayOp calleable. I took the liberty to

• implement it for heterogeneous arrays (so when left and right hand have different element type)
• implement it for the case where the right-hand side is not an array, in which case the scalar operand will be applied to every element of the left-hand-side
• I didn't support
  • in-place operations
  • array ref/array (const) view
  • arrays of differing shapes or dimensionality

Here goes:

template <typename Op> struct ArrayOp {
    Op op;
    explicit ArrayOp(Op op) : op(op) {}

    template <typename T, typename Scalar, size_t Dim> auto operator()(
        boost::multi_array<T, Dim> const& l,
        Scalar const& v) const
    {
        std::array<int, Dim> shape;
        std::copy_n(l.shape(), Dim, shape.data());

        using R = boost::multi_array<decltype(op(T{}, v)), Dim>;
        R result(shape);
        
        std::transform(
           l.data(), l.data()+l.num_elements(),
           result.data(),
           [&op=op,v](auto const& el) { return op(el, v); });

        return result;
    }

    template <typename T, typename U, size_t Dim> auto operator()(
        boost::multi_array<T, Dim> const& l,
        boost::multi_array<U, Dim> const& r) const
    {
        std::array<int, Dim> shape;
        std::copy_n(l.shape(), Dim, shape.data());
        assert(std::equal(shape.begin(), shape.end(), r.shape()));

        using R = boost::multi_array<decltype(op(T{}, U{})), Dim>;
        R result(shape);
        
        std::transform(
           l.data(), l.data()+l.num_elements(),
           r.data(), result.data(),
           [&op=op](auto const& v1, auto const& v2) { return op(v1, v2); });

        return result;
    }
};

Basically it comes down to

• deduce resulting array element type and shape
• do a unary or binary transform (depending on scalar/array rhs)

Now we can write the program:

Live On Compiler Explorer

int main() {
    using MA = boost::multi_array<int, 2>;

    auto shape = boost::extents[3][3];
    MA array1(shape), array2(shape);

    std::generate_n(array1.data(), array1.num_elements(),
            [n = 0]() mutable { return n+=100; });
    std::generate_n(array2.data(), array2.num_elements(),
            [n = 0]() mutable { return n+=1; });

    fmt::print("array1:\n\t{}\n", fmt::join(array1,"\n\t"));
    fmt::print("array2:\n\t{}\n", fmt::join(array2,"\n\t"));

    using namespace ArrayOperators;
    auto array3 = (array1 + array2)/100.0;
    fmt::print("array3:\n\t{}\n", fmt::join(array3,"\n\t"));
}

And it prints

array1:
    {100, 200, 300}
    {400, 500, 600}
    {700, 800, 900}
array2:
    {1, 2, 3}
    {4, 5, 6}
    {7, 8, 9}
array3:
    {1.01, 2.02, 3.03}
    {4.04, 5.05, 6.06}
    {7.07, 8.08, 9.09}

BUT WAIT, WHAT ARE YOU SOLVING

• If you want matrix (not "array") operations use Boost uBlas, Eigen, Armadillo
• If you want utmost perf, using SIMD/AVX2/GPU instructions, you can use Boost Compute