3

我正在从 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 解决方案帖子的结果

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

  • 考虑到 {Row,Column}-Major 模式来存储矩阵时的结果：
    • 在 Row-Major 模式下循环内循环中的行时为 226.56 毫秒（案例 1）。
    • 以 Row-Major 模式在内部循环中循环列时为 223.22 毫秒（案例 2）。
    • 使用 GSL 提供的功能时为 224.60 ms gsl_matrix_complex_set（案例 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);
    }
}
4

4 回答 4

5

没有办法避免迭代所有元素并exp()在每个元素上调用或等效。但是有更快和更慢的迭代方式。

特别是，您的目标应该是最小化缓存未命中。找出您的数据是以行优先还是列优先的顺序存储的，并确保安排循环，以便内部循环遍历内存中连续存储的元素，而外部循环则大步前进到下一行（如果行专业）或列（如果列专业）。尽管这看起来微不足道，但它可以在性能上产生巨大的差异（取决于矩阵的大小）。

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

最后，加入一些循环展开（即，为循环的每次迭代执行 8 或 16 个元素）以进一步减少循环开销，这可能与您能做到的一样快。您可能需要一个带有贯穿的最终 switch 语句来清理剩余元素（当您的数组大小 % 块大小！= 0 时）。

于 2010-07-24T02:48:27.980 回答
3

不，除非有一些我没有听说过的奇怪的数学怪癖，否则你几乎只需要用两个 for 循环遍历元素。

于 2010-07-23T21:38:58.460 回答
2

如果你只想应用exp到一个数字数组，真的没有捷径可走。你必须调用它 (Nx * Ny) 次。如果一些矩阵元素很简单，比如 0，或者有重复的元素，一些记忆可能会有所帮助。

但是，如果您真正想要的是矩阵指数（这非常有用），我们依赖的算法是DGPADM。它在 Fortran 中，但您可以使用f2c将其转换为 C。这是关于它的论文。

于 2010-07-24T02:04:45.763 回答
0

由于没有显示循环的内容，计算c_value的位我们不知道代码的性能是受内存带宽限制还是受CPU限制。唯一确定的方法是使用分析器，而且是一个复杂的分析器。它需要能够测量内存延迟，即 CPU 空闲等待数据从 RAM 到达的时间量。

如果您受到内存带宽的限制，那么一旦您按顺序访问内存，您就无能为力了。当按顺序获取数据时，CPU 和内存工作得最好。随机访问会影响吞吐量，因为数据更有可能必须从 RAM 中提取到缓存中。您总是可以尝试获得更快的 RAM。

如果您受到 CPU 的限制，那么您可以使用更多选项。使用 SIMD 是一种选择，手动编码浮点代码也是如此（由于许多原因，C/C++ 编译器在 FPU 代码方面表现不佳）。如果这是我，并且内部循环中的代码允许这样做，我将有两个指向数组的指针，一个在开始处，另一个在 4/5 处。每次迭代，将使用第一个指针执行 SIMD 操作，使用第二个指针执行标量 FPU 操作，以便循环的每次迭代执行五个值。然后，我将 SIMD 指令与 FPU 指令交错以降低延迟成本。这应该不会影响您的缓存，因为（至少在 Pentium 上）MMU 可以同时传输多达四个数据流（即

于 2010-07-25T23:12:54.723 回答