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

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} - 存储矩阵的主要模式时的结果:
  • 在Row-Major模式下循环内部循环中的行时为226.56 ms(情况1).
  • 在Row-Major模式下循环内循环中的列时为223.22 ms(情况2).
  • 使用gsl_matrix_complex_setGSL提供的功能时224.60 ms (案例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()在每个元素上调用或等效.但是有更快更慢的迭代方式.

特别是,您的目标应该是最大限度地减少缓存未命中.找出你的数据是以行主要顺序还是按列主顺序存储,并确保循环使得内部循环迭代在内存中连续存储的元素,并且外部循环将大步向下一行( if row major)或column(如果是major major).虽然这看起来微不足道,但它可以在性能上产生巨大的差异(取决于矩阵的大小).

处理完缓存后,您的下一个目标是消除循环开销.第一步(如果您的矩阵API支持它)是从嵌套循环(M&N边界)到迭代基础数据(M N界限)的单个循环.你需要获得一个指向底层内存块的原始指针(即双倍而不是双倍**)才能做到这一点.

最后,抛出一些循环展开(也就是说,为循环的每次迭代做8或16个元素)以进一步减少循环开销,这可能就像你可以做到的那样快.你可能需要一个带有fall-through的最终switch语句来清理其余的元素(当你的数组大小为%block size!= 0时).