c**********e
发帖数: 2007 | 1 主要问题是，两维数组如何在函数里实现。
一个矩阵当然可以用两维数组实现。比如说，我要做一个矩阵加法。
我可以定义三个两维数组
int** x=new int*[m];
for(int i=0;i
int** y=new int*[m];
for(int i=0;i
int** z=new int*[m];
for(int i=0;i
这个函数似乎可以定义为
int** matSum(int** x, int** y) {
for(int i=0;i
for(int j=0;j
z[i][j]=x[i][j]+y[i][j];
}
}
return z;
}
这样做显然是有问题的。大牛们怎么解决这个问题的呢？ | t****t
发帖数: 6806 | 2 那你说说, 这个函数显然有什么问题呢?
【在 c**********e 的大作中提到】
: 主要问题是，两维数组如何在函数里实现。
: 一个矩阵当然可以用两维数组实现。比如说，我要做一个矩阵加法。
: 我可以定义三个两维数组
: int** x=new int*[m];
: for(int i=0;i: int** y=new int*[m];
: for(int i=0;i: int** z=new int*[m];
: for(int i=0;i: 这个函数似乎可以定义为
| b********e
发帖数: 58 | 3 why don't define a matrix class? you will have row number, column number,
data type, data as it members. | C***m
发帖数: 120 | 4 你说的问题是m,n的size吗？
是不是可以把矩阵定义成class,然后把相应的information加进class
我自己练习写矩阵运算时候就这么弄的。
【在 c**********e 的大作中提到】
: 主要问题是，两维数组如何在函数里实现。
: 一个矩阵当然可以用两维数组实现。比如说，我要做一个矩阵加法。
: 我可以定义三个两维数组
: int** x=new int*[m];
: for(int i=0;i: int** y=new int*[m];
: for(int i=0;i: int** z=new int*[m];
: for(int i=0;i: 这个函数似乎可以定义为
| d****n
发帖数: 1637 | 5 okay,没代码没真相。
here is my shabby codes and I want to share with you all.
output and analysis will be post soon
//file name matrix_op.c
#include
#include
#include
#include
#include
#include
#define _MX_INIT(type_t, row, col) ({ int i;\
type_t ** ret;\
ret=(type_t ** )malloc( row *sizeof(type_t*));\
for(i=0;i
(ret);\
})
#define _MX_FREE(mx, row) ({int i; for(i=0;i
free(mx)
;})
#define _MX_GEN(type_t, row, col, randfunc )({\
int i, j;\
type_t **ret;\
ret=_MX_INIT(type_t, row,col);\
for(i=0;i
for(j=0;j
ret[i][j]=(randfunc) ;\
(ret);\
})
int **mx_add_naive(int row, int col, int *const *a, int *const *b)
{
int i, j;
int **m = _MX_INIT(int, row, col);
for (i = 0; i < row; ++i)
for (j = 0; j < col; ++j)
m[i][j] = a[i][j] + b[i][j];
return m;
}
int **mx_add_better(int row, int col, int *const *a, int *const *b)
{
int i, j;
int **m = _MX_INIT(int, row, col);
for (i = 0; i < row; ++i) {
int *ap = a[i], *bp = b[i], *mp = m[i]; //let compiler know where
it
can be buffered
for (j = 0; j < col; ++j)
mp[j] = ap[j] + bp[j];
}
return m;
}
int **mx_add_sse2(int row, int col, int *const *a, int *const *b)
{
int i, j, k;
int width = 4;
int **m = _MX_INIT(int, row, col);
int x[width];
for (i = 0; i < row; ++i) {
int *ap = a[i], *bp = b[i], *mp = m[i];
__m128i t;
for (j = 0; j < width * (col / width); j += width) {
//the function here is for 64 bits integet operation.
//if need only 32 bits, then change the function names
correspondingly.
t = _mm_add_epi64(_mm_loadu_si128((__m128i *) (ap + j)),
_mm_loadu_si128((__m128i *) (bp + j)));
_mm_store_si128((__m128i *) x, t);
mp[j] = x[0], mp[j + 1] = x[1]; //no bound check for simple
}
//finishing the leftovers
for (j = width * (col / width); j < col; j++) {
m[i][j] = a[i][j] + b[i][j];
}
}
return m;
}
int main(int argc, char *argv[])
{
clock_t t;
int col = 100, row = 100;
if (argc > 1)
row = col = atoi(argv[1]);
int **a, **b, **m;
srand(time(NULL));
a = _MX_GEN(int, row, col, rand());
b = _MX_GEN(int, row, col, rand());
t = clock();
m = mx_add_naive(row, col, a, b);
fprintf(stderr, "naive add: %lf second; m[%d][%d]=%d \n",
(double) (clock() - t) / CLOCKS_PER_SEC, row / 2, col / 2,
m[row / 2][col / 2]);
_MX_FREE(m, row);
t = clock();
m = mx_add_better(row, col, a, b);
fprintf(stderr, "better add: %lf second; m[%d][%d]=%d \n",
(double) (clock() - t) / CLOCKS_PER_SEC, row / 2, col / 2,
m[row / 2][col / 2]);
_MX_FREE(m, row);
t = clock();
m = mx_add_sse2(row, col, a, b);
fprintf(stderr, "sse2 add: %lf second; m[%d][%d]=%d \n",
(double) (clock() - t) / CLOCKS_PER_SEC, row / 2, col / 2,
m[row / 2][col / 2]);
_MX_FREE(m, row);
_MX_FREE(a, row);
_MX_FREE(b, row);
} | d****n
发帖数: 1637 | 6 //note : all output are in my system
1. compiled as
gcc matrix_op.c -lm -o matrix_op -msse2
[No -O3 and -msse2 yet]
commandline: ./matrix_op 10000
output:
naive add: 0.910000 second; m[5000][5000]=-1850619072
better add: 0.720000 second; m[5000][5000]=-1850619072
sse2 add: 0.670000 second; m[5000][5000]=-1850619072
2. compiled as :
gcc matrix_op.c -lm -o matrix_op_sse -O3 -msse2
commandline: ./matrix_op 10000
output:
naive add: 0.600000 second; m[5000][5000]=-1162569274
better add: 0.590000 second; m[5000][5000]=-1162569274
sse2 add: 0.610000 second; m[5000][5000]=-1162569274 | d****n
发帖数: 1637 | 7 My comments:
1. compiler level optimization
in function "mx_add_bette" implementation specific function
there is some restrictions but it works better than "naive" version
2. sse2
sse2 is better than non-sse for sure.
optimization is still buggy and not very generic typing.
3. gcc -O3 plus -msse2 can do better than my hand coding
save alot time even using a naive approach.
4. if matrix size is small(<1000 in my system), the results may change.
I think this is because the sse instruction overhead need to be
considered.
5. multiplications will make more sense using different approach.
welcome to give your suggestions and post your results. I would like to
study more as well.
thanks | d****n
发帖数: 1637 | 8 here are more output from my machine.
I realized that for a large matrix size it eats too much memory.
$ ./mx_int32 20000
naive add: 3.650000 second; m[10000][10000]=-1256136404
better add: 3.100000 second; m[10000][10000]=-1256136404
sse2 add: 2.970000 second; m[10000][10000]=-1256136404
$ ./mx_int32 30000
naive add: 7.990000 second; m[15000][15000]=1948902602
better add: 6.350000 second; m[15000][15000]=1948902602
sse2 add: 5.790000 second; m[15000][15000]=1948902602
$ ./mx_int32 40000
naive add: 14.300000 second; m[20000][20000]=-202077897
better add: 11.490000 second; m[20000][20000]=-202077897
sse2 add: 10.780000 second; m[20000][20000]=-202077897
$ ./mx_int32 50000
naive add: 23.310000 second; m[25000][25000]=2016613794
better add: 19.410000 second; m[25000][25000]=2016613794
sse2 add: 18.480000 second; m[25000][25000]=2016613794 | t****t
发帖数: 6806 | 9 你这是10年前的写法了, 现在没人这么写, 都靠编译器了. 你也看到了, 你做得并没有
比编译器好, 而且我可以告诉你, 你并没有给编译器完整的提示.
【在 d****n 的大作中提到】
: okay,没代码没真相。
: here is my shabby codes and I want to share with you all.
: output and analysis will be post soon
: //file name matrix_op.c
: #include
: #include
: #include
: #include
: #include
: #include
| d****n
发帖数: 1637 | 10 again, 没代码没真相。provide your code here
PS: I am not showing how good it is, but how this works.
【在 t****t 的大作中提到】
: 你这是10年前的写法了, 现在没人这么写, 都靠编译器了. 你也看到了, 你做得并没有
: 比编译器好, 而且我可以告诉你, 你并没有给编译器完整的提示.
| | | t****t
发帖数: 6806 | 11 你不是有code了么, 我就随便跑一跑. by specifying too much detail on
implementation, you actually INTERFERE with compiler optimization.
$ gcc -O3 -fprefetch-loop-arrays -march=native -funroll-loops -ffast-math 21
.c
$ a.out 20000
naive add: 0.920000 second; m[10000][10000]=2069822843
better add: 0.930000 second; m[10000][10000]=2069822843
sse2 add: 0.990000 second; m[10000][10000]=2069822843
【在 d****n 的大作中提到】
: again, 没代码没真相。provide your code here
: PS: I am not showing how good it is, but how this works.
| t****t
发帖数: 6806 | 12 you are actually showing how this DOES NOT work...
【在 d****n 的大作中提到】
: again, 没代码没真相。provide your code here
: PS: I am not showing how good it is, but how this works.
| d****n
发帖数: 1637 | 13 "how the optimization and coding style affect the speed" works.
Since earlier ppl talking about sse optimization things, I just made my
implementation here.
I am not good enough to beat gnu-gcc developers. but in the right way.
If you read carefully and understand my shabby English above,
as a non-psyche developer, who would like to do such way ?
BTW:
I am not interested arguing with you.
again, do a full illustration not just pop up argument with no support,
which I and I believe most ppl here dont like that.
【在 t****t 的大作中提到】
: you are actually showing how this DOES NOT work...
| t****t
发帖数: 6806 | 14 i've shown the illustration, didn't you see? by specifying correct compiler
switch, your most "NAIVE" implementation runs the FASTEST; your most "
advanced" implementation runs the slowest. the conclusion is, adding too
much detail interferes the optimization.
btw i am not gcc developer, merely a common gcc user.
i will show you what is a correct compiler hint in a while. i am not arguing
with you either, but i am sorry, that's not the right way.
【在 d****n 的大作中提到】
: "how the optimization and coding style affect the speed" works.
: Since earlier ppl talking about sse optimization things, I just made my
: implementation here.
: I am not good enough to beat gnu-gcc developers. but in the right way.
: If you read carefully and understand my shabby English above,
: as a non-psyche developer, who would like to do such way ?
: BTW:
: I am not interested arguing with you.
: again, do a full illustration not just pop up argument with no support,
: which I and I believe most ppl here dont like that.
| d****n
发帖数: 1637 | 15 Sure I know you are not gcc developer.
here is my result
$ gcc -O3 -fprefetch-loop-arrays -funroll-loops -ffast-math matrix_op.c
-lm
$ ./a.out 20000
naive add: 2.200000 second; m[10000][10000]=1836335890
better add: 2.200000 second; m[10000][10000]=1836335890
sse2 add: 2.280000 second; m[10000][10000]=1836335890
Be honest, I never play with these fancy flags using gcc.
and I think the flags step on the mine.
##############man page for gcc####################
-fprefetch-loop-arrays
If supported by the target machine, generate instructions to
prefetch memory to improve the performance of loops that access large arrays.
These options may generate better or worse code; results are
highly dependent on the structure of loops within the source code.
###################################################
gcc is such a good tool only optimize the naive code, isnt it?
how fun ! | t****t
发帖数: 6806 | 16 Add the following function into comparison. Note, this is almost the same as
"naive" version, but i added hint: "restrict", meaning a and b are not
alias of each other.
int **mx_add_restrict(int row, int col, int * restrict const * restrict a,
int * restrict const * restrict b)
{
int i, j;
int **m = _MX_INIT(int, row, col);
for (i = 0; i < row; ++i)
for (j = 0; j < col; ++j)
m[i][j] = a[i][j] + b[i][j];
return m;
}
compile with -std=c99 since restrict is a c99 keyword. Note: this is not
necessary faster than naive version. my point is, that's the correct
type of compiler HINT you should use in your code. you HINT the compiler
by supporting the knowledge it doesn't know, e.g. a and b does not overlap.
you do not, do the job FOR the compiler, because in 99% of cases, you
can not do better.
$ a.out 20000
naive add: 0.950000 second; m[10000][10000]=2005599955
restrict add: 0.940000 second; m[10000][10000]=2005599955
better add: 0.950000 second; m[10000][10000]=2005599955
sse2 add: 1.030000 second; m[10000][10000]=2005599955
【在 d****n 的大作中提到】
: okay,没代码没真相。
: here is my shabby codes and I want to share with you all.
: output and analysis will be post soon
: //file name matrix_op.c
: #include
: #include
: #include
: #include
: #include
: #include
| p**o
发帖数: 3409 | | d****n
发帖数: 1637 | 18 Okay, learned - > restrict.(Sorry thrust, I have to say this is lame.)
Conclusion:
-O3 for the win.
Got to go back to work.
same
as
not
restrict a,
【在 t****t 的大作中提到】
: Add the following function into comparison. Note, this is almost the same as
: "naive" version, but i added hint: "restrict", meaning a and b are not
: alias of each other.
: int **mx_add_restrict(int row, int col, int * restrict const * restrict a,
: int * restrict const * restrict b)
: {
: int i, j;
: int **m = _MX_INIT(int, row, col);
: for (i = 0; i < row; ++i)
: for (j = 0; j < col; ++j)
| p**o
发帖数: 3409 | 19 矩阵运算还是调BLAS/Lapack好了，提供了C的接口，调起来也方便。
考虑了对不同cpu的各种优化策略（尤其是对mem bus的优化），
光是矩阵相乘(level-3 BLAS)的实现就有十几万行fortran代码，
你自己裸写拿来练手可以，生产环境基本没法用的。
【在 d****n 的大作中提到】
: Sure I know you are not gcc developer.
: here is my result
: $ gcc -O3 -fprefetch-loop-arrays -funroll-loops -ffast-math matrix_op.c
: -lm
: $ ./a.out 20000
: naive add: 2.200000 second; m[10000][10000]=1836335890
: better add: 2.200000 second; m[10000][10000]=1836335890
: sse2 add: 2.280000 second; m[10000][10000]=1836335890
: Be honest, I never play with these fancy flags using gcc.
: and I think the flags step on the mine.
| d****n
发帖数: 1637 | 20 you are right.
I couldnt write such code here.
but the above is only from understanding point of view.
【在 p**o 的大作中提到】
: 矩阵运算还是调BLAS/Lapack好了，提供了C的接口，调起来也方便。
: 考虑了对不同cpu的各种优化策略（尤其是对mem bus的优化），
: 光是矩阵相乘(level-3 BLAS)的实现就有十几万行fortran代码，
: 你自己裸写拿来练手可以，生产环境基本没法用的。
| | | p**o
发帖数: 3409 | 21 btw, 楼主貌似问的是C++
可以考虑Boost uBLAS，提供了现成的BLAS的C++封装
www.boost.org/libs/numeric | t****t
发帖数: 6806 | 22 其实大家不过是借楼灌水, 按照我的理解, 楼主想问的是二维数组在C/C++里的common
practice calling convention而已, which no one really cares to answer anyway.
【在 p**o 的大作中提到】
: btw, 楼主貌似问的是C++
: 可以考虑Boost uBLAS，提供了现成的BLAS的C++封装
: www.boost.org/libs/numeric
| d**t
发帖数: 183 | 23 Try Eigen which is a nice C++ linear algebra library.
【在 c**********e 的大作中提到】
: 主要问题是，两维数组如何在函数里实现。
: 一个矩阵当然可以用两维数组实现。比如说，我要做一个矩阵加法。
: 我可以定义三个两维数组
: int** x=new int*[m];
: for(int i=0;i: int** y=new int*[m];
: for(int i=0;i: int** z=new int*[m];
: for(int i=0;i: 这个函数似乎可以定义为
| d****n
发帖数: 1637 | 24 I spent time on my replies and i dont think it was such fluid.
At least, it has something "from 10 years ago", as someone said.
who else? | d**t
发帖数: 183 | 25 uBlas is very slow compared to other libraries. If speed is not a problem,
it is a nice choice.
【在 p**o 的大作中提到】
: btw, 楼主貌似问的是C++
: 可以考虑Boost uBLAS，提供了现成的BLAS的C++封装
: www.boost.org/libs/numeric
| p**o
发帖数: 3409 | 26 开源很好，赞一个。如果你要继续优化你的代码的话，提两条建设性意见：
1. 你现在的乘法实现是非常紧的双重循环，矩阵元素是一个一个从RAM通过FSB拿到
L2 cache，再通过BSB拿到L1 cache，再进寄存器来运算的。可以预见的结果就是
cpu core繁忙，而FSB/BSB饥饿。一个聪明的实现，至少应该在RAM中根据L2 cache
的大小把矩阵分块，按块而不是单个元素来向L2 cache里传。类似的，L2 cache和
L1 cache之间、L1 cache和寄存器之间，也应该有相应的I/O优化策略，以及应对
各种cache miss的情形。据我所知，gcc目前的优化帮不上什么忙。
2. 稀疏矩阵怎么处理？是不是应该重新排序，然后分成若干块稠密矩阵来算？
其实这些在很多生产级别的开源库中都实现了，拿来用即可，又快又稳定。
不过如果你对进一步优化你的代码并深入学习感兴趣，不妨考虑一下上述建议。
X86架构下一个简单的实现，只用几千行C代码就可以基本搞定，对于规模稍大的
矩阵可以看到明显的提速。
【在 d****n 的大作中提到】
: I spent time on my replies and i dont think it was such fluid.
: At least, it has something "from 10 years ago", as someone said.
: who else?
| c**********e
发帖数: 2007 | 27 你们把俺的楼弄歪了，也不回答俺的问题。不过还是谢谢诸位大侠灌水，把这个歪楼建得高高的。
不过俺好像明白了。就是用双重指针，加上各维长度，m，n，等。
common
anyway.
【在 t****t 的大作中提到】
: 其实大家不过是借楼灌水, 按照我的理解, 楼主想问的是二维数组在C/C++里的common
: practice calling convention而已, which no one really cares to answer anyway.
| h****g
发帖数: 71 | 28 用库吧。
大矩阵乘法不是直接三重循环的，有更好的方法，低于O(n^3)。 | u****u
发帖数: 229 | 29 你的第一个建议是错误的，除了副作用外不会有任何的好处。你的建议是赤果果的蔑视
广大cpu/memory/cache设计者几个decades以来的努力。
【在 p**o 的大作中提到】
: 开源很好，赞一个。如果你要继续优化你的代码的话，提两条建设性意见：
: 1. 你现在的乘法实现是非常紧的双重循环，矩阵元素是一个一个从RAM通过FSB拿到
: L2 cache，再通过BSB拿到L1 cache，再进寄存器来运算的。可以预见的结果就是
: cpu core繁忙，而FSB/BSB饥饿。一个聪明的实现，至少应该在RAM中根据L2 cache
: 的大小把矩阵分块，按块而不是单个元素来向L2 cache里传。类似的，L2 cache和
: L1 cache之间、L1 cache和寄存器之间，也应该有相应的I/O优化策略，以及应对
: 各种cache miss的情形。据我所知，gcc目前的优化帮不上什么忙。
: 2. 稀疏矩阵怎么处理？是不是应该重新排序，然后分成若干块稠密矩阵来算？
: 其实这些在很多生产级别的开源库中都实现了，拿来用即可，又快又稳定。
: 不过如果你对进一步优化你的代码并深入学习感兴趣，不妨考虑一下上述建议。
| t****t
发帖数: 6806 | 30 喔, 有人有不同意见, 太好了, 搬个板凳前排学习中
【在 u****u 的大作中提到】
: 你的第一个建议是错误的，除了副作用外不会有任何的好处。你的建议是赤果果的蔑视
: 广大cpu/memory/cache设计者几个decades以来的努力。
| | | N***m
发帖数: 4460 | 31 不带你这么翅果果的调戏人家的
【在 t****t 的大作中提到】
: 喔, 有人有不同意见, 太好了, 搬个板凳前排学习中
| t****t
发帖数: 6806 | 32 你们这些人哪, 就是TSSN!
这方面我确实不熟, 我是真心来学习的.
【在 N***m 的大作中提到】
: 不带你这么翅果果的调戏人家的
| N***m
发帖数: 4460 | 33 what is TSSN?
【在 t****t 的大作中提到】
: 你们这些人哪, 就是TSSN!
: 这方面我确实不熟, 我是真心来学习的.
| t****t
发帖数: 6806 | 34 难道你们没有听过江core的名言么...
【在 N***m 的大作中提到】
: what is TSSN?
| N***m
发帖数: 4460 | 35 ...
高啊！
【在 t****t 的大作中提到】
: 难道你们没有听过江core的名言么...
| p**o
发帖数: 3409 | 36
表达观点要有论据或论证。我不是做体系结构这块的，说错话、甚至
讲外行话的可能性也有。我在这里把我的依据以及理解稍微展开一下，
希望你反驳的时候言之有据，无论对错大家都能受益。
“广大cpu/memory/cache设计者几个decades以来的努力”，的确为
上层提供了一些必要的intrinsics，但各级I/O策略是具体依赖上层
软件来控制实现的，硬件没有聪明到能自动优化这些。
不要说硬件，就连编译器都很难识别上层的语义做优化。比如GCC，
就连对矩阵indexing的处理都非常笨拙；Intel的ICC/IFC可以向量化
内循环，从而更方便用SSE指令集做优化，但是没有专门对矩阵相关
的各级I/O做优化——这些优化是Intel放在MKL数值库里做的。
我的根据之一是这篇科普文章（文章虽长，但值得每个程序员一读）：
Ulrich Drepper, "What Every Programmer Should Know About Memory"
http://people.freebsd.org/~lstewart/articles/cpumemory.pdf
文章前部讲RAM/cache电路实现的可以先快速扫过，第6章可详细读。
首先，以p.49那个naive的矩阵乘法的C代码为例：
for (i = 0; i < N; ++i)
for (j = 0; j < N; ++j)
for (k = 0; k < N; ++k)
res[i][j] += mul1[i][k] * mul2[k][j];
经过gcc编译后，mul1和mul2都是行优先的内存布局，这导致对
mul2的内存访问不是sequential的，造成大幅性能下降。
一个解决办法是把mul2转置，变成p.49右下的代码，这个要靠
程序员，通用编译器和底层硬件还没有聪明到能帮忙做这些。
关于L1 cache的优化，p.50提到怎么设gcc参数来hardcode
cache line size (CLS)。C代码要围绕这个CLS来写，程序变成
6重循环。这些都是写矩阵库的程序员自己要关心的事情。
（p.97有这个6重循环的比较完整的实例代码）
至于L2/L3 cache的优化，主要思想就是根据cache大小给矩阵
分块计算，以及在此基础上设计合适的I/O策略。摘录p.59几句话：
"To avoid the high costs of cache misses, the working
set size should be matched to the cache size."
"A program has to perform its job even if the data set
is too large. It is the programmer’s job to do the work
in a way which minimizes cache misses. For last-level
caches this is possible–just as for L1 caches–by working
on the job in smaller pieces."
"What this means is that the code must dynamically adjust
itself to the cache line size. This is an optimization
specific to the program."
这些事情都要写矩阵库的程序员自己来忙活。
第二个依据是GotoBLAS库的作者自述其设计原则的文章，通过阅读
我大致了解其各级cache levels下的分块策略，以及自己徒手实现
的“不可能性”（实在对矩阵库实现有兴趣的同学不妨也读一读）：
K. Goto and R. Geijn,
"Anatomy of high-performance matrix multiplication," ACM TOMS'08.
http://dl.acm.org/citation.cfm?id=1356053
【在 u****u 的大作中提到】
: 你的第一个建议是错误的，除了副作用外不会有任何的好处。你的建议是赤果果的蔑视
: 广大cpu/memory/cache设计者几个decades以来的努力。
| t****t
发帖数: 6806 | 37 听上去很有道理(note: 我没说很对或很不对), 但是你提出观点的时候针对的是dryden
先前的sample code, which is a simple element-by-element add. 这跟乘法的
access pattern有本质区别. 所以你的论据和论点关系不大.
【在 p**o 的大作中提到】
:
: 表达观点要有论据或论证。我不是做体系结构这块的，说错话、甚至
: 讲外行话的可能性也有。我在这里把我的依据以及理解稍微展开一下，
: 希望你反驳的时候言之有据，无论对错大家都能受益。
: “广大cpu/memory/cache设计者几个decades以来的努力”，的确为
: 上层提供了一些必要的intrinsics，但各级I/O策略是具体依赖上层
: 软件来控制实现的，硬件没有聪明到能自动优化这些。
: 不要说硬件，就连编译器都很难识别上层的语义做优化。比如GCC，
: 就连对矩阵indexing的处理都非常笨拙；Intel的ICC/IFC可以向量化
: 内循环，从而更方便用SSE指令集做优化，但是没有专门对矩阵相关
| x*******1
发帖数: 28835 | 38 他第一条有什么错？ blas不就是block cache + SSE2么？
【在 u****u 的大作中提到】
: 你的第一个建议是错误的，除了副作用外不会有任何的好处。你的建议是赤果果的蔑视
: 广大cpu/memory/cache设计者几个decades以来的努力。
| x*******1
发帖数: 28835 | 39 你讲的ijk, ikj... 不同的顺序造成的访问memory locality不一样在GCC上编译出的
code是对的。 如果用icc, ifort, 编出来的都是一样的性能。 icc认出这块code做什
么， 直接translate最优的code了。
【在 p**o 的大作中提到】
:
: 表达观点要有论据或论证。我不是做体系结构这块的，说错话、甚至
: 讲外行话的可能性也有。我在这里把我的依据以及理解稍微展开一下，
: 希望你反驳的时候言之有据，无论对错大家都能受益。
: “广大cpu/memory/cache设计者几个decades以来的努力”，的确为
: 上层提供了一些必要的intrinsics，但各级I/O策略是具体依赖上层
: 软件来控制实现的，硬件没有聪明到能自动优化这些。
: 不要说硬件，就连编译器都很难识别上层的语义做优化。比如GCC，
: 就连对矩阵indexing的处理都非常笨拙；Intel的ICC/IFC可以向量化
: 内循环，从而更方便用SSE指令集做优化，但是没有专门对矩阵相关
| p**o
发帖数: 3409 | 40
dryden
这个确实是我看走眼了... 没注意到dryden的代码只有加法。
对于矩阵加法，两个矩阵的access pattern完全相同，
后一个矩阵不必转置，所以我提的“转置”那一段不适用；
而且对于这种pattern，编译器很容易把循环向量化
（gcc -O3会打开-ftree-vectorize）。
加法对矩阵分块也似乎没有什么要求，cache优化或许可以
完全依赖编译器。回头我查文档确认一下。
【在 t****t 的大作中提到】
: 听上去很有道理(note: 我没说很对或很不对), 但是你提出观点的时候针对的是dryden
: 先前的sample code, which is a simple element-by-element add. 这跟乘法的
: access pattern有本质区别. 所以你的论据和论点关系不大.
| | | u****u
发帖数: 229 | 41 我只是对于你如何得出1建议的前提表示很疑惑.什么叫做"矩阵元素是一个一个从RAM通过FSB拿到L2 cache，再通过BSB拿到L1 cache"?对于一段连续的矩阵元素来说,在稍微现代一点的处理器上是不可能出现这种情况的.再说了,越是性能好的程序就越是应该"cpu core繁忙，而FSB/BSB饥饿",我不懂你为什么说要避免这种情况.让core忙死而FSB闲死可以说是cache设计者的final fantasy.
cache的工作都是硬件执行的,软件一般不应该试图操作cache"按块而不是单个元素来向L2 cache里传".软件的设计应该是试图让存储方式和cache合拍,而不是去控制.现代的cache控制是非常复杂的,在这么一个高度复杂的系统外另加一套控制的做法我一般认为是只会适得其反的.
比如说,你说可以用指令让cache多读点数据进来.问题是数据读的多了,别的东西就cache的少了,你可能就把重要的系统指令段的数据给踢出去了,然后系统狂miss,可能一miss就freeze了系统,然后你就发现所有的程序都停下来了.这种东西太复杂了,谁也说不清的.
【在 p**o 的大作中提到】
:
: dryden
: 这个确实是我看走眼了... 没注意到dryden的代码只有加法。
: 对于矩阵加法，两个矩阵的access pattern完全相同，
: 后一个矩阵不必转置，所以我提的“转置”那一段不适用；
: 而且对于这种pattern，编译器很容易把循环向量化
: （gcc -O3会打开-ftree-vectorize）。
: 加法对矩阵分块也似乎没有什么要求，cache优化或许可以
: 完全依赖编译器。回头我查文档确认一下。
| x*******1
发帖数: 28835 | 42 大哥， 指令cache和数据cache都是分开的。 | t****t
发帖数: 6806 | 43 on x86, L1 are separate while L2 are unified...
【在 x*******1 的大作中提到】
: 大哥， 指令cache和数据cache都是分开的。
| r*********r
发帖数: 3195 | | E*******1
发帖数: 3464 | 45 why not use the overload operator to define m(i,j)? Do not think in a way of
C, you are using C++
class matrix{
const size_t _size;
double * _b;
public:
double operator ()(size_i i, size_t j){return _b{j*_size+i};}
//build some constructor thing below
............
};
main()
{matrix m;
//m(i,j) is the matrix element in ith row and jth column
} | n******t
发帖数: 4406 | 46 大的矩阵转置这种事情，不配合cache做分块，编译器怎么自动优化啊？
通过FSB拿到L2 cache，再通过BSB拿到L1 cache"?对于一段连续的矩阵元素来说,在稍
微现代一点的处理器上是不可能出现这种情况的.再说了,越是性能好的程序就越是应该
"cpu core繁忙，而FSB/BSBbr />
cache"按块而不是单个元素来向L2 cache里传".软件的设计应该是试图让存储方式和
cache合拍,而不是去控制.现代的cache控制是非常复杂的,在这么一个高度复杂的系统
外另加一套控制的做法我一般认为是只会适得其反的.
cache的少了,你可能就把重要的系统指令段的数据给踢出去了,然后系统狂miss,可能一
miss就freeze了系统,然后你就发现所有的程序都停下来了.这种东西太复杂了,谁也说
不清的.
【在 u****u 的大作中提到】
: 我只是对于你如何得出1建议的前提表示很疑惑.什么叫做"矩阵元素是一个一个从RAM通过FSB拿到L2 cache，再通过BSB拿到L1 cache"?对于一段连续的矩阵元素来说,在稍微现代一点的处理器上是不可能出现这种情况的.再说了,越是性能好的程序就越是应该"cpu core繁忙，而FSB/BSB饥饿",我不懂你为什么说要避免这种情况.让core忙死而FSB闲死可以说是cache设计者的final fantasy.
: cache的工作都是硬件执行的,软件一般不应该试图操作cache"按块而不是单个元素来向L2 cache里传".软件的设计应该是试图让存储方式和cache合拍,而不是去控制.现代的cache控制是非常复杂的,在这么一个高度复杂的系统外另加一套控制的做法我一般认为是只会适得其反的.
: 比如说,你说可以用指令让cache多读点数据进来.问题是数据读的多了,别的东西就cache的少了,你可能就把重要的系统指令段的数据给踢出去了,然后系统狂miss,可能一miss就freeze了系统,然后你就发现所有的程序都停下来了.这种东西太复杂了,谁也说不清的.
| E*******1
发帖数: 3464 | 47 sorry, I messed up in typing brackets
in the class the overload operator () should be obviously written as
double operator ()(size_i i, size_t j){return _b[j*_size+i];}
Then you can access to the matrix element as in Fortran
of
【在 E*******1 的大作中提到】
: why not use the overload operator to define m(i,j)? Do not think in a way of
: C, you are using C++
: class matrix{
: const size_t _size;
: double * _b;
: public:
: double operator ()(size_i i, size_t j){return _b{j*_size+i};}
: //build some constructor thing below
: ............
: };
|
|
