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?

Ellingston

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;
}

Ellingston

Just stumbled onto this question - and I am also trying to accomplish a similar objective. A few questions came to mind looking at your code: 1. Does your wrapper class only have a vector private member? For a matrix member, did you write another wrapper class? 2. How did you accomplish matrix vector multiplication? 3. Does your implementation preserve the efficiency of expression templates? – Polonaise

1. Yes. Or a Row, or a Matrix private member. You can either write a wrapper for each class (Col, Row, Mat), or just use a Matrix for everything and write a single wrapper. 2. You can write multiple wrappers, or a single wrapper, as I said in 1. Then you would need to overload operators. 3. It's efficient, but never as efficient as the raw implementation. I would say we lose anywhere between 1 and 10% efficiency in computation time. What we gain is the ability to switch between libraries. Right now, we can easily switch between Eigen and Armadillo, without changing the public API. – Ellingston

I have been trying to write a single Matrix wrapper class for both matrices and vectors just so I can put the overloaded operators in one wrapper file. However, Eigen does not make it easy to do so. I would like to have just one private member - an uninitialized Matrix<T, Dynamic, Dynamic> - that can be resized to be both a MxN matrix or a Mx1 vector, but Eigen vector operations (like dot prod) don't work because a vector needs to be initialized as a Matrix<T, Dynamic, 1>. – Polonaise

You can still have a single Matrix wrapper. However, for some functions, you will need to convert the Eigen matrix into an Eigen vector before calling the Eigen API. This is never exposed to the user, it is done inside the function. – Ellingston