已知矩阵尺寸,最快迭代矩阵的方法是什么?

3
我想知道在c/c++中迭代矩阵的最快方法是什么。
到目前为止,我想到的最好方法是将矩阵映射到一维数组。
然后使用指针算术运算,是否有其他更快的方法?
维度在运行时已知但不在编译时已知,矩阵已完全填充。
#include <iostream>
#include <time.h>
#define XMAX 500
#define YMAX 400
#define ZMAX 300

int main()
{
    srand(0);   
    register double sum = 0;
    register int i;
    register int j;
    register int k;

    double *arr_ptr;
    arr_ptr = new double[XMAX*YMAX*ZMAX];

    for (i=0; i<XMAX*YMAX*ZMAX; ++i)
    {
        *(arr_ptr+i) = rand()/double(RAND_MAX);
    }

    clock_t start, finish;
    start = clock();

    for (i=0; i<XMAX; ++i)
    {
        for (j=0; j<YMAX; ++j)
        {
            for (k=0; k<ZMAX; ++k)
            {
                sum += *(arr_ptr+i*YMAX*ZMAX+j*ZMAX+k);
            }
        }
    }

    finish = clock();
    std::cout << "sum: " << sum << "\telapsed: " << finish - start << std::endl;
    std::cin.get();

    delete[] arr_ptr;
}

2
你说得对。将其表示为连续的一维数组,并使用线性索引进行迭代。 - Dima
2
我认为 *(array+i)array[i] 快不了多少... 编译器会以相同的方式翻译它。 - Simon
1
@Simon是正确的,array[i]被定义为与*(array + i)完全相同,因此这两种形式之间没有区别。 - caf
8个回答

7
为什么要费心去写三个嵌套的for循环?你可以直接这样做:
for (i=0; i<XMAX*YMAX*ZMAX; ++i)
{
    sum += *(arr_ptr+i);
}

或者只需 sum += arr_ptr[i]。请注意,对于非常大的数组,XMAX*YMAX*ZMAX 可能会导致类型溢出(在这种情况下,将其转换为 size_t 是适当的)。 - Stephen Canon
嵌套的循环是为了模拟预期的使用方式,应该更加清晰地说明这一点。 - darckeen
@user615174:你正在按线性顺序访问基本上是线性缓冲区;我认为单个for循环非常清晰。 - Stephen Canon

3

根据ideone.com编译器的测试结果,针对XMAX 500YMAX 400ZMAX 100运行100次,这比您的代码快了650毫秒。

double *p_current, *p_end;

p_current = arr_ptr;
p_end = (arr_ptr + XMAX*YMAX*ZMAX);
while(p_current != p_end) {
    sum += *p_current++;
}

请查看:旧版本新版本


2

实际上这并不重要,因为编译器会进行优化。所以arr[i][j][k]*(arr_ptr+i*YMAX*ZMAX+j*ZMAX+k)的速度是相同的。


一般而言并非如此。连续的数组将从缓存局部性中受益。动态分配的嵌套数组不会是连续的,因此您将遇到由于缓存未命中而产生的运行时开销。当然,如果在编译时已知维度并创建了静态数组,则该数组将是连续的,那么您是正确的。 - Konrad Rudolph
是的,我尝试了三维数组,但单维度在迭代时大约快了40%,如果考虑到释放内存,单维度的速度就像快了1000%。 - darckeen
编译时无法确定尺寸,只能在运行时确定。 - darckeen
@user615174:是什么让你觉得数组不会是连续的? - Jens Gustedt
只有在使用C99并且数组在堆栈上分配为连续数组时,此答案才是正确的...但那时它将是arr,而不是arr_ptr。 - Jim Balter

1

OpenCV使用指针算术:

double *ptr = arr_ptr;
for (i=0; i<XMAX*YMAX*ZMAX; ++i)
{
    sum += *ptr++;
}

我猜可能会更快一些。试试看,然后展示一下时间吧!


这里没有理由指针算术应该更快。编译器应该为索引访问和指针代码生成完全相同的代码。另一方面,在循环外提升乘法可能是一个好主意。 - Konrad Rudolph
1
如果你在谈论 XMAX*YMAX*ZMAX 的乘法,那么这个肯定会被预处理器简化! - Simon
您所提到的乘法是指XMAX * YMAX * ZMAX吗?我一直很好奇,大多数编译器是否会通过预先计算来优化掉这个乘法,因为XMAX YMAX和ZMAX扩展为整数常量。如果它们是实际变量但在循环中没有更改,我想知道大多数编译器是否会预先计算结果? - Ken Wayne VanderLinde
它更快,并且提供与我提出的解决方案相同的性能增益 =) - Trinidad

1
double *ptr = arr_ptr;
for (int i=XMAX*YMAX*ZMAX; i>0; --i)
{
    sum += *ptr++;
} 

将循环变量与零进行比较,而不是与某个常数进行比较,可以为每次迭代节省一两个时钟周期(例如,在 Intel CPU 上使用 JNZ 指令)


我按照我的答案尝试了同样的方法,与零比较没有明显的收益,但当然比原来快。http://www.ideone.com/OhsS8 - Trinidad
不要这样做。在现代系统上,int通常是32位,而size_t是64位。初始值很容易溢出。 - Jens Gustedt

1
在您的示例中,边界是常量,因此正常的三维数组可以使用,无论是C还是C++。
然后,就动态分配具有可变边界的数组而言,C和C++是真正不同的语言,请不要混淆它们。对于C ++,请使用向量类等东西。它们是为此而制作的,应该很有效。
在C中,自C99以来就有VLA,即可变长度数组。与城市传说相反,如果您不将它们分配到堆栈上,则它们可能非常高效。像在C中为任何大块内存一样使用malloc。
double (*arr_ptr)[XMAX][YMAX][ZMAX]
  = malloc(sizeof(*arr_ptr));

for (register size_t i=0; i<XMAX; ++i)
  for (register size_t j=0; j<YMAX; ++j)
    for (register size_t k=0; k<ZMAX; ++k)
       (*arr_ptr)[i][j][k] = rand()/double(RAND_MAX);

.

free(arr_ptr);

现代处理器具有相当复杂的寻址方案,因此可能没有必要有效地进行完整的索引计算。你的编译器通常比你更懂。为了高效,更重要的是如何声明和处理循环变量。使用正确的类型进行索引,`size_t` 是正确的无符号类型。当计算三维扁平化索引并在此处使用有符号类型时,`int` 可能会很容易溢出并且没有太多意义。将这些变量尽可能声明为本地变量,可以使事情更清晰。`register` 只是与编译器签订协议,您永远不会获取这种索引的地址。通常这并不能太大程度地改进效率。但是,当您修改代码时,这可能会防止您做出低效的事情。最后但并非最不重要的是,如果您真的担心效率,请检查您的编译器生成了什么。例如,`gcc` 有选项`-S`来生成中间汇编程序。读取它而不是猜测效率。

0
需要说的第一件事是,堆栈分配的多维数组在内存中(在C和C++中)以行主序存储。也就是说,matrix[2][2] = {{1,2},{3,4}}将被存储在内存中,就像您实际声明了一个array[4] = {1,2,3,4},而matrix[][]语法只是*(matrix + i * 2 + j)的语法糖。
因此,遍历矩阵的最快方法取决于您如何遍历它:按行主序或列主序,并且矩阵的大小如何:
- 如果整个矩阵可以适合CPU缓存,那么遍历顺序并不重要; - 如果矩阵比CPU缓存大,则进行按行主顺序遍历会导致更少的CPU缓存未命中。
了解矩阵操作是否存在性能问题以及原因的最佳方法是对代码进行性能分析。

0
对于非常大的数据块,请考虑并行操作。在这种情况下,可以使用 gather 操作计算总和 -- 其形式将取决于您选择的并行框架。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接