evol128[Blog]

I am the bone of my code

为什么转置一个512x512的矩阵,会比513x513的矩阵慢很多?

evol128 posted @ 2012年9月05日 12:38 in arichitecture with tags architecture cache , 10788 阅读

谨以此文,纪念刚退休的Professor Sibert以及Professor Goel。你们尽管已年过70,却还仍然坚持在教导学生,实在令人钦佩。我今天所拥有的编程知识,经验,技巧,很大一部分是从你们那儿学来的。谢谢你们。

 

问题的出处:http://stackoverflow.com/questions/11413855/why-is-transposing-a-matrix-of-512x512-much-slower-than-transposing-a-matrix-of

事情的起因是这样的,先看下面这段代码:

 

 

 

define SAMPLES 1000
#define MATSIZE 512

#include <time.h>
#include <iostream>
int mat[MATSIZE][MATSIZE];

void transpose()
{
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
   {
       int aux = mat[i][j];
       mat[i][j] = mat[j][i];
       mat[j][i] = aux;
   }
}

int main()
{
   //initialize matrix
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
       mat[i][j] = i+j;

   int t = clock();
   for ( int i = 0 ; i < SAMPLES ; i++ )
       transpose();
   int elapsed = clock() - t;

   std::cout << "Average for a matrix of " << MATSIZE << ": " << elapsed / SAMPLES;

 

很普通的一个求矩阵转置的程序。但是,当MATSIZE取512和513的时候,出现了非常有意思的结果:

512 平均 2.19ms

513 平均 0.57ms

很让人惊讶吧,513竟然比512快。更进一步的研究发现,size=512的时候,运算速度会比同数量级的其它数字慢很多很多。这是怎么一回事呢?

stackoverflow上大牛给的解答非常正确,但是这次,我不想做翻译了。我从Professor Sibert那里,从Professor Goel那里,学到的知识,足够帮我解决这个问题了,我不是一个人。下面是我的解答:

很容易就联想到,造成这个问题的原因是CPU cache,我们有很多种方式来存储cache,具体可以参考这里

原作者没有给出他的CPU型号,但是如今的pc几乎都是采用的set associative的cache结构,下面我用2-way set associate来做例子,讲解一下cache的工作原理。

(图片取自Professor Sibert的讲义,这可是纯ascii画的哦= =)


一个内存地址,可以划分为block,tag,word,byte 4个部分。10bits的block,对应了1024个cache set,内存地址的block固定了,就必须存储在相应的set里面,这样可以把查询cache的事件从O(n)缩短为O(1)。

举个例子,block是1023(1111111111),你的数据就放在第1023个set里面。可能有人会觉得奇怪,为什么block不是取的最前面的10bits,这当然是有道理的,通常在内存里数据都是连续存放的,就是说,同一段程序用的数据,他们前10位几乎都是一样的,如果用前10位来定位block,那么collision的发生率非常高,cache效率非常低下,所以才选了后面的10位来定位block。

当然,每个set里面有多条记录,2-way是2条,你得遍历这两条记录,比较前面50位的tag,如果tag一样,并且Valid bit(V)=1,那么恭喜你,你的数据在cache里面,接着就可以通过word和byte来取数据了。

如果遍历完这两条记录,还是没有找到tag的话,那么很遗憾,你的数据不在cache里,得从内存里读。从内存里获取相应的数据,然后把它存到对应的cache set里,如果set里有空位的话最好,如果没有的话,用LRU来替换。因为一个set里只有2条数据,所以实现LRU仅仅需要一个额外bit就可以了,非常高效。

好了,背景知识介绍的差不多了,让我们回到这个问题上来。为什么512大小的矩阵,会比其它数字慢那么多?

让我们来计算一下,512x512的int矩阵,在内存里是连续存放的。每个cache line是16bytes,对应4个int,所以一个n阶矩阵的row可以填充n/4个cache set。假设第一个数据a[0][0]正好对应cache set 0,那么其中每一个数据a[i][j]对应的cache set是(512i+j)/4%1024=(128i+j/4)%1024。可以看到,前面的系数正好可以整除。很不巧的是,在进行矩阵转置的运算时,在第2个for循环中,我们需要依次访问每一个row中对应i的值。这样会造成下面的结果:假设i=0,set(a[0][0])=0, set(a[1][0])=128, set(a[2][0])=256...set(a[7][0])=896,set(a[0][0])=0,后面开始重复了,到a[15][0]的时候刚好填完整个cache的所有128整数倍的set,当读取a[16][0]的时候,将会发生replace,把a[0][0]从cache里移除。这样,当源程序的i=1时,将完全重复i=0的计算过程,每次取数据都需要先从memory读到cache中来,cache的作用完全没有体现。

而当size=513的时候,事情就不一样了,mat[i][j]对应的cache set是(513i+j)/4%1024,前面的系数除不尽了,每递增4次结果会比size=512时偏差1。例如:set(a[0][0])=0, set(a[1][0])=128, set(a[2][0])=256,set(a[3][0])=384, set(a[4][0])=513...这样就很微妙的把cache set给错开了。a[16][0]不在第0行而是第4行,不会覆盖之前的数据。即使将全部的a[0-15][i]都读入cache,也不会发生碰撞。之后,由于一个cache有4个word,a[0-15][i+1],a[0-15][i+2],a[0-15][i+3]也同时被读进cache里了,所以计算i+1,i+2,i+3时,仅仅需要读对应行的数据就可以了,同一行的数据都是连续的,所以碰撞率很低。这个计算过程很好的利用了cache,如果不考虑其他因素(实际上,这个已经是影响运行时间的最大因素了),理论上我们可以节省75%的运行时间,可以看到,这个理论预测是和提问者给的数据相符合的。

总之,当你的data size是128的整数倍的时候,得特别小心,搞不好cache collision就把你的程序给拖慢了呢

 

Update 1: 原代码有逻辑错误,这点大家都不要吐槽了,代码不是我写的= =

Update 2:帅哥问我,为什么可以加速这么多。这个循环包括4次读cache的操作,2次写cache的操作,以及0-2次replace操作。每次replace操作会有一次memory read,有可能会有memory write(假设它是write back)。前面的读写cache时间和读写内存相比,几乎可以忽略,对效率产生显著影响的是后面的内存读写。如果cache的hit率高了,那么内存读写的次数就少了,程序运行时间是会产生很大影响的

Update 3:当然,具体效果还视乎CPU架构而定,我自己试验的只有节省25%左右时间

Update 4: 有人提出了用划分矩阵(把大矩阵分成若干个小矩阵分别计算)的方法来求转置。划分矩阵可以解决类似的问题(譬如说求两个矩阵乘积),但是对解决这个问题没有任何帮助。因为求转置的时候,每个数据只用到了一次,没有重复访问;即便划分成更小的矩阵,在cache里面的位置也没有发生改变。

Update 5: 据说,Professor Goel只是因病休息几个学期,没有退休。。。(原来你还要回来教课!!!)

 

 

Avatar_small
nonoob 说:
2012年9月08日 12:18

前几天看到一个由于指令预取导致效率的差异,这次又看到cache命中带来的差异 ==

Avatar_small
Arch 说:
2012年9月20日 22:25

好害怕啊,c编程这么步步惊心啊

Avatar_small
zmx 说:
2012年9月24日 22:11

我机器上, 513的比512的慢....

Avatar_small
evol128 说:
2012年9月25日 05:32

@zmx: 这个比较奇怪,能把操作系统,CPU,编译器信息发给我看看么?

Avatar_small
zmx 说:
2012年9月25日 10:55

@evol128: Ubuntu11.10 gcc 4.6.1 Intel(R) Core(TM) i3 CPU 530 @ 2.93GHz
512个元素的矩阵, 时间是2.47 , 513的时间是2.58

Avatar_small
evol128 说:
2012年9月26日 01:33

@zmx: 我刚找了redhat + i5 + gcc试一下,并没有出现你这种结果。可以把你用的程序以及gcc -S的结果发给我看一下吗

Avatar_small
zmx 说:
2012年9月26日 14:03

@evol128: 恩,我找了别的几台i3的机器,也没有这样的结果...

.file "test_cache.cpp"
.local _ZStL8__ioinit
.comm _ZStL8__ioinit,1,1
.globl mat
.bss
.align 32
.type mat, @object
.size mat, 1048576
mat:
.zero 1048576
.text
.globl _Z9transposev
.type _Z9transposev, @function
_Z9transposev:
.LFB966:
.cfi_startproc
pushl %ebp
.cfi_def_cfa_offset 8
.cfi_offset 5, -8
movl %esp, %ebp
.cfi_def_cfa_register 5
subl $16, %esp
movl $0, -12(%ebp)
jmp .L2
.L5:
movl $0, -8(%ebp)
jmp .L3
.L4:
movl -12(%ebp), %eax
sall $9, %eax
addl -8(%ebp), %eax
movl mat(,%eax,4), %eax
movl %eax, -4(%ebp)
movl -8(%ebp), %eax
sall $9, %eax
addl -12(%ebp), %eax
movl mat(,%eax,4), %edx
movl -12(%ebp), %eax
sall $9, %eax
addl -8(%ebp), %eax
movl %edx, mat(,%eax,4)
movl -8(%ebp), %eax
sall $9, %eax
addl -12(%ebp), %eax
movl -4(%ebp), %edx
movl %edx, mat(,%eax,4)
addl $1, -8(%ebp)
.L3:
cmpl $511, -8(%ebp)
setle %al
testb %al, %al
jne .L4
addl $1, -12(%ebp)
.L2:
cmpl $511, -12(%ebp)
setle %al
testb %al, %al
jne .L5
leave
.cfi_restore 5
.cfi_def_cfa 4, 4
ret
.cfi_endproc
.LFE966:
.size _Z9transposev, .-_Z9transposev
.section .rodata
.LC1:
.string "Average for a matrix of "
.LC2:
.string ": "
.text
.globl main
.type main, @function
main:
.LFB967:
.cfi_startproc
pushl %ebp
.cfi_def_cfa_offset 8
.cfi_offset 5, -8
movl %esp, %ebp
.cfi_def_cfa_register 5
andl $-16, %esp
subl $64, %esp
movl $0, 44(%esp)
jmp .L7
.L10:
movl $0, 48(%esp)
jmp .L8
.L9:
movl 48(%esp), %eax
movl 44(%esp), %edx
addl %eax, %edx
movl 44(%esp), %eax
sall $9, %eax
addl 48(%esp), %eax
movl %edx, mat(,%eax,4)
addl $1, 48(%esp)
.L8:
cmpl $511, 48(%esp)
setle %al
testb %al, %al
jne .L9
addl $1, 44(%esp)
.L7:
cmpl $511, 44(%esp)
setle %al
testb %al, %al
jne .L10
call clock
movl %eax, 56(%esp)
movl $0, 52(%esp)
jmp .L11
.L12:
call _Z9transposev
addl $1, 52(%esp)
.L11:
cmpl $999, 52(%esp)
setle %al
testb %al, %al
jne .L12
call clock
subl 56(%esp), %eax
movl %eax, 60(%esp)
fildl 60(%esp)
flds .LC0
fdivrp %st, %st(1)
fstps 28(%esp)
movl $.LC1, 4(%esp)
movl $_ZSt4cout, (%esp)
call _ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_PKc
movl $512, 4(%esp)
movl %eax, (%esp)
call _ZNSolsEi
movl $.LC2, 4(%esp)
movl %eax, (%esp)
call _ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_PKc
flds 28(%esp)
fstps 4(%esp)
movl %eax, (%esp)
call _ZNSolsEf
movl $_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_, 4(%esp)
movl %eax, (%esp)
call _ZNSolsEPFRSoS_E
movl $0, %eax
leave
.cfi_restore 5
.cfi_def_cfa 4, 4
ret
.cfi_endproc
.LFE967:
.size main, .-main
.type _Z41__static_initialization_and_destruction_0ii, @function
_Z41__static_initialization_and_destruction_0ii:
.LFB977:
.cfi_startproc
pushl %ebp
.cfi_def_cfa_offset 8
.cfi_offset 5, -8
movl %esp, %ebp
.cfi_def_cfa_register 5
subl $24, %esp
cmpl $1, 8(%ebp)
jne .L13
cmpl $65535, 12(%ebp)
jne .L13
movl $_ZStL8__ioinit, (%esp)
call _ZNSt8ios_base4InitC1Ev
movl $_ZNSt8ios_base4InitD1Ev, %eax
movl $__dso_handle, 8(%esp)
movl $_ZStL8__ioinit, 4(%esp)
movl %eax, (%esp)
call __cxa_atexit
.L13:
leave
.cfi_restore 5
.cfi_def_cfa 4, 4
ret
.cfi_endproc
.LFE977:
.size _Z41__static_initialization_and_destruction_0ii, .-_Z41__static_initialization_and_destruction_0ii
.type _GLOBAL__sub_I_mat, @function
_GLOBAL__sub_I_mat:
.LFB978:
.cfi_startproc
pushl %ebp
.cfi_def_cfa_offset 8
.cfi_offset 5, -8
movl %esp, %ebp
.cfi_def_cfa_register 5
subl $24, %esp
movl $65535, 4(%esp)
movl $1, (%esp)
call _Z41__static_initialization_and_destruction_0ii
leave
.cfi_restore 5
.cfi_def_cfa 4, 4
ret
.cfi_endproc
.LFE978:
.size _GLOBAL__sub_I_mat, .-_GLOBAL__sub_I_mat
.section .ctors,"aw",@progbits
.align 4
.long _GLOBAL__sub_I_mat
.weakref _ZL20__gthrw_pthread_oncePiPFvvE,pthread_once
.weakref _ZL27__gthrw_pthread_getspecificj,pthread_getspecific
.weakref _ZL27__gthrw_pthread_setspecificjPKv,pthread_setspecific
.weakref _ZL22__gthrw_pthread_createPmPK14pthread_attr_tPFPvS3_ES3_,pthread_create
.weakref _ZL20__gthrw_pthread_joinmPPv,pthread_join
.weakref _ZL21__gthrw_pthread_equalmm,pthread_equal
.weakref _ZL20__gthrw_pthread_selfv,pthread_self
.weakref _ZL22__gthrw_pthread_detachm,pthread_detach
.weakref _ZL22__gthrw_pthread_cancelm,pthread_cancel
.weakref _ZL19__gthrw_sched_yieldv,sched_yield
.weakref _ZL26__gthrw_pthread_mutex_lockP15pthread_mutex_t,pthread_mutex_lock
.weakref _ZL29__gthrw_pthread_mutex_trylockP15pthread_mutex_t,pthread_mutex_trylock
.weakref _ZL31__gthrw_pthread_mutex_timedlockP15pthread_mutex_tPK8timespec,pthread_mutex_timedlock
.weakref _ZL28__gthrw_pthread_mutex_unlockP15pthread_mutex_t,pthread_mutex_unlock
.weakref _ZL26__gthrw_pthread_mutex_initP15pthread_mutex_tPK19pthread_mutexattr_t,pthread_mutex_init
.weakref _ZL29__gthrw_pthread_mutex_destroyP15pthread_mutex_t,pthread_mutex_destroy
.weakref _ZL30__gthrw_pthread_cond_broadcastP14pthread_cond_t,pthread_cond_broadcast
.weakref _ZL27__gthrw_pthread_cond_signalP14pthread_cond_t,pthread_cond_signal
.weakref _ZL25__gthrw_pthread_cond_waitP14pthread_cond_tP15pthread_mutex_t,pthread_cond_wait
.weakref _ZL30__gthrw_pthread_cond_timedwaitP14pthread_cond_tP15pthread_mutex_tPK8timespec,pthread_cond_timedwait
.weakref _ZL28__gthrw_pthread_cond_destroyP14pthread_cond_t,pthread_cond_destroy
.weakref _ZL26__gthrw_pthread_key_createPjPFvPvE,pthread_key_create
.weakref _ZL26__gthrw_pthread_key_deletej,pthread_key_delete
.weakref _ZL30__gthrw_pthread_mutexattr_initP19pthread_mutexattr_t,pthread_mutexattr_init
.weakref _ZL33__gthrw_pthread_mutexattr_settypeP19pthread_mutexattr_ti,pthread_mutexattr_settype
.weakref _ZL33__gthrw_pthread_mutexattr_destroyP19pthread_mutexattr_t,pthread_mutexattr_destroy
.section .rodata
.align 4
.LC0:
.long 1315859240
.ident "GCC: (Ubuntu/Linaro 4.6.1-9ubuntu3) 4.6.1"
.section .note.GNU-stack,"",@progbits

Avatar_small
evol128 说:
2012年9月27日 08:10

@zmx: 我编译出来的assembly和这个totally different...我需要时间仔细check一下

Avatar_small
nwpumabian 说:
2012年10月18日 11:11

512太小了吧。换成1024试试。
原因可能是cache line size不同吧。

Avatar_small
foool 说:
2013年9月22日 19:14

你这说的有问题吧,

Avatar_small
coolgod 说:
2013年11月11日 20:38

1、

"set(a[1][0])=128, set(a[2][0])=256...set(a[7][0])=896,set(a[0][0])=0,后面开始重复了"

---这里是笔误吗,应该是从set(a[8][0])开始就重复了吧?应该是从set(a[8][0])就填完了所有128倍数的set了吧?

 

2、

"当源程序的i=1时,将完全重复i=0的计算过程"

——这里也没看懂,i=1?笔误吗?源程序j表示列,i表示行的呀。这里的i是指列吧?

 

Avatar_small
后 说:
2017年8月03日 17:07

这两教授是哪个学校的

Avatar_small
gaodq 说:
2017年12月19日 20:47

跟cache line有关,512 x 512的矩阵是1MB数据,在我这里L2 cache有1MB,L3 cache 4MB,按照博主描述的那种现象,如果矩阵小点的话,小于1MB,性能差别就不是很明显的


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter