D&C的C实现比Naive矩阵乘法解更快？

Question

Divide&Conquer(D&C)解决方案和用于矩阵乘法的Naive解决方案都是使用C编程语言"就地"实现的.所以根本没有动态内存分配.

正如我们已经了解了两种解决方案,它们实际上具有相同的时间复杂度,即O(n ^ 3).现在它们共享相同的空间复杂性,因为它们都是就地实现的.那怎么可能比另一个快得多？

使用clock_gettime来获取时间.

凭借核心i7笔记本电脑的Windows 7上的Cygwin,D&C解决方案的运行速度比Naive解决方案(删除冗余日志)快得多:

编辑:

"algo0"表示Naive解决方案,而"algo1"表示D&C解决方案.

"len"表示矩阵的宽度和高度.矩阵是NxN矩阵.

"00:00:00:000:003:421"表示:"小时:分钟:秒:毫秒:微秒:纳赛克".

[alg0]time cost[0, len=00000002]: 00:00:00:000:003:421 (malloc_cnt=0)
[alg1]time cost[0, len=00000002]: 00:00:00:000:000:855 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[1, len=00000004]: 00:00:00:000:001:711 (malloc_cnt=0)
[alg1]time cost[1, len=00000004]: 00:00:00:000:001:711 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[2, len=00000008]: 00:00:00:000:009:408 (malloc_cnt=0)
[alg1]time cost[2, len=00000008]: 00:00:00:000:008:553 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[3, len=00000016]: 00:00:00:000:070:134 (malloc_cnt=0)
[alg1]time cost[3, len=00000016]: 00:00:00:000:065:858 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[4, len=00000032]: 00:00:00:000:564:066 (malloc_cnt=0)
[alg1]time cost[4, len=00000032]: 00:00:00:000:520:873 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[5, len=00000064]: 00:00:00:004:667:337 (malloc_cnt=0)
[alg1]time cost[5, len=00000064]: 00:00:00:004:340:188 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[6, len=00000128]: 00:00:00:009:662:680 (malloc_cnt=0)
[alg1]time cost[6, len=00000128]: 00:00:00:008:139:403 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[7, len=00000256]: 00:00:00:080:031:116 (malloc_cnt=0)
[alg1]time cost[7, len=00000256]: 00:00:00:065:395:329 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[8, len=00000512]: 00:00:00:836:392:576 (malloc_cnt=0)
[alg1]time cost[8, len=00000512]: 00:00:00:533:799:924 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[9, len=00001024]: 00:00:09:942:086:780 (malloc_cnt=0)
[alg1]time cost[9, len=00001024]: 00:00:04:307:021:362 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[10, len=00002048]: 00:02:53:413:046:992 (malloc_cnt=0)
[alg1]time cost[10, len=00002048]: 00:00:35:588:289:832 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[11, len=00004096]: 00:25:46:154:930:041 (malloc_cnt=0)
[alg1]time cost[11, len=00004096]: 00:04:38:196:205:661 (malloc_cnt=0)

即使在只有一个ARM内核的Raspberry Pi上,结果也是类似的(同样,删除了冗余数据):

[alg0]time cost[0, len=00000002]: 00:00:00:000:005:999 (malloc_cnt=0)
[alg1]time cost[0, len=00000002]: 00:00:00:000:051:997 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[1, len=00000004]: 00:00:00:000:004:999 (malloc_cnt=0)
[alg1]time cost[1, len=00000004]: 00:00:00:000:008:000 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[2, len=00000008]: 00:00:00:000:014:999 (malloc_cnt=0)
[alg1]time cost[2, len=00000008]: 00:00:00:000:023:999 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[3, len=00000016]: 00:00:00:000:077:996 (malloc_cnt=0)
[alg1]time cost[3, len=00000016]: 00:00:00:000:157:991 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[4, len=00000032]: 00:00:00:000:559:972 (malloc_cnt=0)
[alg1]time cost[4, len=00000032]: 00:00:00:001:248:936 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[5, len=00000064]: 00:00:00:005:862:700 (malloc_cnt=0)
[alg1]time cost[5, len=00000064]: 00:00:00:010:739:450 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[6, len=00000128]: 00:00:00:169:060:336 (malloc_cnt=0)
[alg1]time cost[6, len=00000128]: 00:00:00:090:290:373 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[7, len=00000256]: 00:00:03:207:909:599 (malloc_cnt=0)
[alg1]time cost[7, len=00000256]: 00:00:00:771:870:443 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[8, len=00000512]: 00:00:35:725:494:551 (malloc_cnt=0)
[alg1]time cost[8, len=00000512]: 00:00:08:139:712:988 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[9, len=00001024]: 00:06:29:762:101:314 (malloc_cnt=0)
[alg1]time cost[9, len=00001024]: 00:01:50:964:568:907 (malloc_cnt=0)
------------------------------------------------------
[alg0]time cost[10, len=00002048]: 00:52:03:950:717:474 (malloc_cnt=0)
[alg1]time cost[10, len=00002048]: 00:14:19:222:020:444 (malloc_cnt=0)

我的第一个猜测是它必须由GCC完成一些优化.但具体如何？

以下是Naive解决方案和D&C解决方案的代码.天真的解决方案:

void ClassicalMulti(int const * const mat1,
                    int const * const mat2,
                    int * const matrix,
                    const int n) {
    if (!mat1 || !mat2 || n<=0) {
        printf("ClassicalMulti: Invalid Input\n");
        return;
    }

    int cnt, row, col;

    for (row=0;row<n;++row) {
        for (col=0;col<n;++col) {
            for (cnt=0;cnt<n;++cnt) {
                matrix[row*n+col] += mat1[row*n+cnt] * mat2[cnt*n+col];
            }
        }
    }
}

分而治之的解决方案:

void DCMulti(int const * const mat1,
             int const * const mat2,
             int * const matrix,
             const int p1,
             const int p2,
             const int pn,
             const int n) {
    if (!mat1 || !mat2 || !matrix || n<2 || p1<0 || p2 <0 || pn<2) {
        printf("DCMulti: Invalid Input\n");
        return;
    }

    if (pn == 2) {
        int pos = (p1/n)*n + p2%n;
        matrix[pos]     += mat1[p1]*mat2[p2] + mat1[p1+1]*mat2[p2+n];
        matrix[pos+1]   += mat1[p1]*mat2[p2+1] + mat1[p1+1]*mat2[p2+1+n];
        matrix[pos+n]   += mat1[p1+n]*mat2[p2] + mat1[p1+1+n]*mat2[p2+n];
        matrix[pos+1+n] += mat1[p1+n]*mat2[p2+1] + mat1[p1+1+n]*mat2[p2+1+n];
    } else {
        int a = p1;
        int b = p1 + pn/2;
        int c = p1 + pn*n/2;
        int d = p1 + pn*(n+1)/2;
        int e = p2;
        int f = p2 + pn/2;
        int g = p2 + pn*n/2;
        int h = p2 + pn*(n+1)/2;
        DCMulti(mat1, mat2, matrix, a, e, pn/2, n);   // a*e
        DCMulti(mat1, mat2, matrix, b, g, pn/2, n);   // b*g
        DCMulti(mat1, mat2, matrix, a, f, pn/2, n);   // a*f
        DCMulti(mat1, mat2, matrix, b, h, pn/2, n);   // b*h 
        DCMulti(mat1, mat2, matrix, c, e, pn/2, n);   // c*e 
        DCMulti(mat1, mat2, matrix, d, g, pn/2, n);   // d*g 
        DCMulti(mat1, mat2, matrix, c, f, pn/2, n);   // c*f 
        DCMulti(mat1, mat2, matrix, d, h, pn/2, n);   // d*h 
    }
}

Answer 1

这两种方法的不同之处仅在于存储器访问模式.即缓存局部性; 对于大型矩阵,尤其是行竞争相同的高速缓存行并且导致对高速缓存未命中的越来越大的惩罚.最后,D&C -strategy得到了回报,尽管全局更好的方法是将问题分成8x8块 - 一种称为循环平铺的技术.(毫不奇怪,矩阵乘法在维基百科文章中作为拱形示例呈现......)