4

我正在尝试将我以前用 python 编写的一些代码转移到 C++ 中,我目前正在测试 xtensor 以查看它是否可以比 numpy 更快地完成我需要的操作。

我的一个函数采用一个方阵 d 和一个标量 alpha,并执行元素运算alpha/(alpha+d)。背景:此函数用于测试哪个值alpha是“最佳”,因此它处于一个循环中,d始终相同,但alpha变化。

以下所有时间尺度都是运行该函数的 100 个实例的平均值。

在 numpy 中,执行此操作大约需要 0.27 秒,代码如下:

def kfun(d,alpha):
    k = alpha /(d+alpha)
    return k

但是 xtensor 大约需要 0.36 秒,代码如下所示:

xt::xtensor<double,2> xk(xt::xtensor<double,2> d, double alpha){
    return alpha/(alpha+d);
}

我也尝试过使用以下版本,std::vector但这是我不想长期使用的东西,即使它只花了 0.22 秒。

std::vector<std::vector<double>> kloops(std::vector<std::vector<double>> d, double alpha, int d_size){
    for (int i = 0; i<d_size; i++){
        for (int j = 0; j<d_size; j++){
            d[i][j] = alpha/(alpha + d[i][j]);
        }
    }
    return d;
}

我注意到operator/in xtensor 使用“延迟广播”,有没有办法让它立即生效?

编辑:

在 Python 中,函数调用如下,并使用“time”包进行计时

t0 = time.time()
for i in range(100):
    kk = k(dsquared,alpha_squared)
print(time.time()-t0)

在 C++ 中,我调用函数如下,并使用 chronos 计时:

//d is saved as a 1D npy file, an artefact from old code
auto sd2 = xt::load_npy<double>("/path/to/d.npy");

shape = {7084, 7084};
    xt::xtensor<double, 2> xd2(shape);
    for (int i = 0; i<7084;i++){
        for (int j=0; j<7084;j++){
            xd2(i,j) = (sd2(i*7084+j));
        }
    }

auto start = std::chrono::steady_clock::now();
for (int i = 0;i<10;i++){
    matrix<double> kk = kfun(xd2,4000*4000,7084);
}
auto end = std::chrono::steady_clock::now();
std::chrono::duration<double> elapsed_seconds = end-start;
std::cout << "k takes: " << elapsed_seconds.count() << "\n";

如果您希望运行此代码,我建议使用xd2对称的 7084x7084 随机矩阵,对角线为零。

函数的输出,一个名为 的矩阵k,然后继续用于其他函数,但我仍然需要d保持不变,因为它稍后会被重用。

结束编辑

要运行我的 C++ 代码,我在终端中使用以下行:

cd "/path/to/src/" && g++ -mavx2 -ffast-math -DXTENSOR_USE_XSIMD -O3 ccode.cpp -o ccode -I/path/to/xtensorinclude && "/path/to/src/"ccode

提前致谢!

4

1 回答 1

3

C++ 实现的一个问题可能是它创建了一个甚至两个可以避免的临时副本。第一个副本来自不通过引用(或完美转发)传递参数。如果不查看其余代码,很难判断这是否对性能有影响。如果保证在方法之后不使用,编译器可能会d进入方法xk(),但更有可能将数据复制到d.

要通过引用传递,可以将方法更改为

xt::xtensor<double,2> xk(const xt::xtensor<double,2>& d, double alpha){
    return alpha/(alpha+d);
}

要使用完美转发(并且还启用其他 xtensor 容器,例如xt::xarrayor xt::xtensor_fixed),可以将方法更改为

template<typename T>
xt::xtensor<double,2> xk(T&& d, double alpha){
    return alpha/(alpha+d);
}

此外,您可以避免为返回值保留内存。同样,如果不查看其余代码,很难判断。但是,如果该方法在循环内使用,并且返回值始终具有相同的形状,那么在循环外创建返回值并通过引用返回可能是有益的。为此,可以将方法更改为:

template<typename T, typename U>
void xk(T& r, U&& d, double alpha){
    r = alpha/(alpha+d);
}

如果保证d并且r不指向相同的内存,则可以进一步换rxt::noalias()以避免在分配结果之前进行临时复制。如果您不通过引用返回,则函数的返回值也是如此。

祝你好运,编码愉快!

于 2021-03-19T18:21:51.627 回答