How to write a third-party library wrapper class around expression templates

Question

We are trying to implement a new C++ code in my research group to perform large numerical simulations (finite elements, finite difference methods, topology optimization, etc.) The software will be used by people from academia and industry alike.

For the dense linear algebra piece of the software, we want to use either Eigen or Armadillo. We wish to build a wrapper around these packages for two reasons: 1. to expose our own API to the users rather than the third-party API; and 2. in case we need to switch libraries in the future. I understand reason 2 is a very expensive form of insurance, but we encountered this situation with our previous simulation software.

The information I have encountered regarding wrapping third-party libraries comes from these sources:

• Should third-party types be exposed in my C++ library's API

• https://softwareengineering.stackexchange.com/questions/107338/using-third-party-libraries-always-use-a-wrapper

My question relates as to the best way to build this wrapper class. Ideally, a thin layer wrapper would be the best, as:

template< typename T >
class my_vec {
private:
    arma::Col< T > _arma_vec;
};

or its equivalent with an Eigen vector.

Then, my class would call the third-party library class as:

my_vec::foo() { return _arma_vec.foo(); }

I think (and I would like confirmation on this) that the issue with this thin layer is that I lose the speed gained from the expression templates these libraries have implemented under the hood. For example, in Armadillo, the following operation:

// Assuming these vectors were already populated.
a =  b + c + d;

becomes something like this:

for ( std::size_t i = 0; i < a.size(); ++i ) {
    a[i] = b[i] + c[i] + d[i];
}

without creating any temporaries due to their implementation of expression templates. The same situation applies to Eigen.

As far as I undertand, the reason I lose the power of expression templates is that while Armadillo or Eigen do not create temporaries of their own, my class my_vec does. The only way to circumvent this would be to build a thin layer wrapper around their expression templates as well. However, at this point, this seems to be a violation of the YAGNI principle.

This related question here:

• How to integrate a library that uses expression templates?

suggests using something like:

my_vec a, b, c;
// ... populate vectors
a._arma_vec = b._arma_vec + c._arma_vec;

Is it possible to use something like this instead?

template< typename T >
arma::Col< T > &
my_vec< T >::data() { return _arma_vec; }

a.data() = b.data() + c.data();

Or use some operator overloading to hide data() from the user? What other alternatives are there if we do not wish to use the libraries directly? Using macros? Using aliases if we decide to use C++11?

Or what would be the most convenient way to build this wrapper class?

Answer 1

Just for future reference, this is how I decided to implement my solution: I overloaded the operator+ in the following way:

template< typename T1, typename T2 >
auto
operator+(
        const my_vec< T1 > & X,
        const my_vec< T2 > & Y ) ->decltype( X.data() + Y.data() )
{
    return X.data() + Y.data();
}

template< typename T1, typename T2 >
auto
operator+(
        const my_vec< T1 > & X,
        const T2 &           Y ) ->decltype( X.data() + Y )
{
    return X.data() + Y;
}

template< typename T1, typename T2 >
auto
operator+(
        const T1 &           X,
        const my_vec< T2 > & Y ) ->decltype( X + Y.data() )
{
    return X + Y.data();
}

Then, I overloaded my operator= in my_vec class with the following:

template< typename T >
template< typename A >
const my_vec< T > &
my_vec< T >::operator=(
        const A & X )
{
    _arma_vec = X;

    return *this;
}