c - 计算矩阵每个元素的指数的最有效方法

Question

我正在从 Matlab 迁移到 C + GSL，我想知道计算矩阵 B 的最有效方法是什么:

B[i][j] = exp(A[i][j])

其中 i 在 [0, Ny] 中，j 在 [0, Nx] 中。

注意这与矩阵指数不同:

B = exp(A)

这可以通过 GSL (linalg.h) 中的一些不稳定/不受支持的代码来完成。

我刚刚找到了蛮力解决方案(几个“for”循环)，但是有没有更聪明的方法来做到这一点？

编辑

来自 Drew Hall 的解决方案帖子的结果

所有结果都来自一个 1024x1024 for(for) 循环，在该循环中，在每次迭代中都分配了两个 double 值(一个复数)。 时间是100次执行的平均时间。

• 考虑到 {Row,Column}-Major 模式来存储矩阵时的结果:
  • 226.56 毫秒，当以行优先模式(情况 1)循环遍历内部循环中的行时。
  • 223.22 毫秒，当以行优先模式(情况 2)在内部循环中循环列时。
  • 使用 GSL 提供的 gsl_matrix_complex_set 函数时为 224.60 毫秒(案例 3)。

案例1的源码:

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix[2*(i*s_tda + j)] = GSL_REAL(c_value);
        matrix[2*(i*s_tda + j)+1] = GSL_IMAG(c_value);
    }
}

案例2的源码:

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix->data[2*(j*s_tda + i)] = GSL_REAL(c_value);
        matrix->data[2*(j*s_tda + i)+1] = GSL_IMAG(c_value);
    }
}

案例3的源码:

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        gsl_matrix_complex_set(matrix, i, j, c_value);
    }
}

Answer 1

无法避免遍历所有元素并对每个元素调用 exp() 或等效方法。但是有更快和更慢的迭代方式。

特别是，您的目标应该是尽量减少缓存未命中。查看您的数据是以行优先还是列优先的顺序存储的，并确保安排您的循环，使内部循环迭代内存中连续存储的元素，并且外部循环大步前进到下一行(如果行主要)或列(如果列主要)。虽然这看起来微不足道，但它可以在性能上产生巨大的差异(取决于矩阵的大小)。

处理完缓存后，您的下一个目标是消除循环开销。第一步(如果您的矩阵 API 支持)是从嵌套循环(M 和 N 边界)到迭代基础数据的单个循环(MN 边界)。您需要获取指向底层内存块的原始指针(即 double 而不是 double **)才能执行此操作。

最后，加入一些循环展开(即，为循环的每次迭代执行 8 或 16 个元素)以进一步减少循环开销，这可能是您能做到的最快速度。您可能需要一个带有 fall-through 的最终 switch 语句来清理剩余元素(当您的数组大小 % block 大小!= 0 时)。