如何围绕表达式模板编写第三方库包装器 class

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

我们正在尝试在我的研究小组中实现一个新的 C++ 代码来执行大型数值模拟(有限元、有限差分法、拓扑优化等)。该软件将被学术界和工业界的人们使用。

对于软件的密集线性代数部分,我们想使用 Eigen 或 Armadillo。我们希望围绕这些包构建一个包装器,原因有二:1. 向用户公开我们自己的 API 而不是第三方 API; 2. 以防我们将来需要切换库。我理解原因2是一种非常昂贵的保险形式,但是我们以前的模拟软件遇到过这种情况。

我遇到的关于包装第三方库的信息来自这些来源:

我的问题涉及构建此包装器的最佳方法 class。理想情况下,薄层包装纸是最好的,因为:

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

或其与特征向量的等价物。

然后,我的class调用第三方库class为:

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

我认为(并且我想确认这一点)这个薄层的问题是我失去了从这些库在引擎盖下实现的表达式模板中获得的速度。例如在犰狳中,如下操作:

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

变成这样:

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

由于 expression templates 的实现而没有创建任何临时对象。同样的情况也适用于Eigen。

据我所知,我失去表达模板能力的原因是 Armadillo 或 Eigen 不会创建它们自己的临时对象,而我的 class my_vec 会。避免这种情况的唯一方法是在它们的表达式模板周围构建一个薄层包装器。然而,在这一点上,这似乎违反了YAGNI原则。

这里有这个相关问题:

建议使用类似的东西:

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

是否可以改用这样的东西?

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

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

或者使用一些运算符重载来对用户隐藏 data()? 如果我们不想直接使用这些库,还有哪些其他选择?使用宏?如果我们决定使用 C++11,是否使用别名?

或者构建此包装器最方便的方法是什么class?

仅供将来参考,这就是我决定实施我的解决方案的方式:我按以下方式重载了 operator+:

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

然后,我在 my_vec class 中重载了我的 operator= 如下:

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

    return *this;
}