求如何优化这段代码

时间:2023-01-03 21:23:04
这个函数的作用就是用mtx来影射rgba的值,把结果保存在rgba数组中
rgba[0]、rgba[1]、rgba[2]、rgba[3]的值分别代表red green blue alpha,范围在0~255
mtx的值本来是小数的,但我让它们乘上256转换为int,所以计算到最后有一个>>8的操作
最终的返回值必须在0~rgba[3]之间,rgba[3]是不变的
void map_rgba(int mtx[4][4], int rgba[4])
{
    int r, g, b;
    r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
    g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba[3] * mtx[3][1];
    b = rgba[0] * mtx[0][2] + rgba[1] * mtx[1][2] + rgba[2] * mtx[2][2] + rgba[3] * mtx[3][2];
    rgba[0] = min(max(0, r >> 8), rgba[3] );
    rgba[1] = min(max(0, g >> 8), rgba[3]);
    rgba[2] = min(max(0, b >> 8), rgba[3]);
}

void main()
{
    const int nTimes = 1000000;
    int mtx[4][4] = {
        0, 1, 2, 3,
        4, 5, 6, 7,
        8, 9, 10,11,
        12, 13, 14, 15
    };
    int rgba[4] = {255, 128, 64, 128};
    clock_t t1 = clock();
    for (int i = 0; i < 1000000; i++)
    {
     map_rgba(mtx, rgba);
    }
    clock_t t2 = clock();

    printf("excute %d times, time = %d\n", nTimes, t2 - t1);
    printf("rgb(%d, %d, %d)\n", rgba[0], rgba[1], rgba[2]);
}
这段代码在我机器(双核,4G内存)上运行耗时为18ms
如何优化map_rgba函数以提高效率呢?
可以使用MMX SSE 或者汇编,请大侠们把完整代码贴上来。(我对MMX SSE 和汇编都不熟悉的)

29 个解决方案

#1


它在VS2010编译出来的代码如何:
void map_rgba(int mtx[4][4], int rgba[4])
{
011F1000  push        ecx  
011F1001  mov         ecx,dword ptr [esp+8] 
    int r, g, b;
    r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
011F1005  mov         edx,dword ptr [eax+20h] 
011F1008  push        ebx  
011F1009  mov         ebx,dword ptr [ecx+4] 
011F100C  push        ebp  
011F100D  mov         ebp,dword ptr [ecx] 
011F100F  push        esi  
011F1010  mov         esi,dword ptr [ecx+8] 
011F1013  imul        edx,esi 
011F1016  push        edi  
011F1017  mov         edi,dword ptr [ecx+0Ch] 
011F101A  mov         ecx,dword ptr [eax+10h] 
011F101D  imul        ecx,ebx 
011F1020  add         ecx,edx 
011F1022  mov         edx,dword ptr [eax+30h] 
011F1025  imul        edx,edi 
011F1028  add         ecx,edx 
011F102A  mov         edx,dword ptr [eax] 
011F102C  imul        edx,ebp 
011F102F  add         ecx,edx 
    g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba[3] * mtx[3][1];
011F1031  mov         edx,dword ptr [eax+4] 
011F1034  imul        edx,ebp 
011F1037  mov         dword ptr [esp+10h],ebp 
011F103B  mov         ebp,dword ptr [eax+14h] 
011F103E  imul        ebp,ebx 
011F1041  add         edx,ebp 
011F1043  mov         ebp,dword ptr [eax+24h] 
011F1046  imul        ebp,esi 
011F1049  mov         esi,dword ptr [eax+34h] 
011F104C  imul        esi,edi 
011F104F  add         edx,ebp 
    b = rgba[0] * mtx[0][2] + rgba[1] * mtx[1][2] + rgba[2] * mtx[2][2] + rgba[3] * mtx[3][2];
011F1051  mov         ebp,dword ptr [eax+18h] 
011F1054  add         edx,esi 
011F1056  imul        ebp,ebx 
011F1059  mov         esi,dword ptr [eax+8] 
011F105C  imul        esi,dword ptr [esp+10h] 
011F1061  mov         ebx,dword ptr [eax+28h] 
011F1064  mov         eax,dword ptr [eax+38h] 
011F1067  add         esi,ebp 
011F1069  imul        eax,edi 
011F106C  mov         ebp,dword ptr [esp+18h] 
011F1070  imul        ebx,dword ptr [ebp+8] 
011F1074  add         esi,ebx 
011F1076  add         esi,eax 
    rgba[0] = min(max(0, r >> 8), rgba[3] );
011F1078  sar         ecx,8 
011F107B  xor         eax,eax 
011F107D  test        ecx,ecx 
011F107F  setl        al   
011F1082  dec         eax  
011F1083  and         eax,ecx 
011F1085  cmp         eax,edi 
011F1087  jge         map_rgba+95h (11F1095h) 
011F1089  xor         eax,eax 
011F108B  test        ecx,ecx 
011F108D  setl        al   
011F1090  dec         eax  
011F1091  and         ecx,eax 
011F1093  jmp         map_rgba+97h (11F1097h) 
011F1095  mov         ecx,edi 
    rgba[1] = min(max(0, g >> 8), rgba[3]);
011F1097  sar         edx,8 
011F109A  mov         eax,ebp 
011F109C  mov         dword ptr [eax],ecx 
011F109E  xor         ecx,ecx 
011F10A0  test        edx,edx 
011F10A2  setl        cl   
011F10A5  dec         ecx  
011F10A6  and         ecx,edx 
011F10A8  cmp         ecx,edi 
011F10AA  jge         map_rgba+0B8h (11F10B8h) 
011F10AC  xor         ecx,ecx 
011F10AE  test        edx,edx 
011F10B0  setl        cl   
011F10B3  dec         ecx  
011F10B4  and         ecx,edx 
011F10B6  jmp         map_rgba+0BAh (11F10BAh) 
011F10B8  mov         ecx,edi 
011F10BA  mov         dword ptr [eax+4],ecx 
    rgba[2] = min(max(0, b >> 8), rgba[3]);
011F10BD  sar         esi,8 
011F10C0  xor         ecx,ecx 
011F10C2  test        esi,esi 
011F10C4  setl        cl   
011F10C7  dec         ecx  
011F10C8  and         ecx,esi 
011F10CA  cmp         ecx,edi 
011F10CC  jge         map_rgba+0E1h (11F10E1h) 
011F10CE  xor         ecx,ecx 
011F10D0  test        esi,esi 
011F10D2  setl        cl   
011F10D5  pop         edi  
011F10D6  dec         ecx  
011F10D7  and         ecx,esi 
011F10D9  pop         esi  
011F10DA  pop         ebp  
011F10DB  mov         dword ptr [eax+8],ecx 
011F10DE  pop         ebx  

#2


hao

#3


可以用PMADDWD(SSE2)、PHADDD(SSSE3)指令。

#4


引用 3 楼 delphiguy 的回复:
可以用PMADDWD(SSE2)、PHADDD(SSSE3)指令。

能否具体点?

#5


建议先弄明白这些指令的含意,然后才知道往哪使劲。

#6


已经很底层了,指令级上估计没什么可优化了,算法应该还可以在优化一点。
比方说求 r g b中的rgba[x]数组,没有必要每次都这样在引用一下把,那么此时必然多一个mov指令,你可以将此值放到一个reg中,下次直接使用。
还有那个矩阵也一样,维数已定,你自己通过偏移add 应该比编译器生成的要快把。
当然了这些代码就必须由你来实现了。关键的是,就算你这样做了,估计也没有什么质的改变。除非还有更好的方法。
一些个人看法。。。

#7


PMADDWD可以一条指令做8个乘法,4个加法。

  r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
  g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba[3] * mtx[3][1];
  b = rgba[0] * mtx[0][2] + rgba[1] * mtx[1][2] + rgba[2] * mtx[2][2] + rgba[3] * mtx[3][2];

这种运算,可以用4条PMADDWD和3条PADDD指令实现(不包括装入数据的指令)。

#8


引用 7 楼 delphiguy 的回复:
PMADDWD可以一条指令做8个乘法,4个加法。

  r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
  g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba……


已经使用这两条指令进行优化了,但是效率并没有明显提高。请问还可以如何进一步优化呢?
#include "emmintrin.h"
void map_rgba_sse(const int mtx[4][4], int rgba[4])
{
    const int *p = mtx[0];
    __m128i mmtx1 = _mm_set_epi16(p[0], p[4], p[1], p[5], p[2], p[6], p[3], p[7]);
    __m128i mmtx2 = _mm_set_epi16(p[8], p[12], p[9], p[13], p[10], p[14], p[11], p[15]);

    __m128i rg = _mm_set1_epi32((rgba[0] << 16) | rgba[1]);
    __m128i ba = _mm_set1_epi32((rgba[2] << 16) | rgba[3]);

    __m128i mrg = _mm_madd_epi16(rg, mmtx1);
    __m128i mba = _mm_madd_epi16(ba, mmtx2);
    __m128i mrgba = _mm_add_epi32(mrg, mba);
    mrgba = _mm_srli_epi32(mrgba, 8);

    rgba[0] = min(max(0, mrgba.m128i_i32[3]), rgba[3] );
    rgba[1] = min(max(0, mrgba.m128i_i32[2]), rgba[3]);
    rgba[2] = min(max(0, mrgba.m128i_i32[1]), rgba[3]);
}
min(max(0, mrgba.m128i_i32[3]), rgba[3] );这些代码又如何优化呢?

#9


优化的非常不错,不看算法我觉得全用协处理器代替处理器运算在现在已经可以处理绝多数情况了。
只是最后一片的跳转会对性能大打折扣,还没有那种超惊人的处理器论文能设计一套摆脱预跳转失误导致的性能损耗呢,特别是预处理了跳转目标后要等待慢腾腾的协处理器来返回自己有没有猜错的处理器更加不喜欢跳转指令了。

暂时想到的二个思路,首先可以继续用指令MINPD

然后呢,可能已经不实用了,记得有一代CPU的介绍里据说已经把预跳转出错处理到20%性能损耗内。
但思路是用赋值代替if,用符号位(是否为负数的标志bit)作为数组的index.有点累,代码明天写试试。

#10


在我应用的过程中,矩阵mtx是固定不变的,而rgba是变化的,所以我改进了程序,效率是最开始那个版本的2倍!
//在外面初始化这些固定的数据,因为map_rgba_sse会被调用非常非常多次的!
const int *p = mtx[0];
__m128i mmtx1 = _mm_set_epi16(p[0], p[4], p[1], p[5], p[2], p[6], p[3], p[7]);
__m128i mmtx2 = _mm_set_epi16(p[8], p[12], p[9], p[13], p[10], p[14], p[11], p[15]);
__m128i lowerbound = _mm_setzero_si128();

void map_rgba_sse(int rgba[4])
{
    __m128i rg = _mm_set1_epi32((rgba[0] << 16) | rgba[1]);
    __m128i ba = _mm_set1_epi32((rgba[2] << 16) | rgba[3]);

    __m128i mrg = _mm_madd_epi16(rg, mmtx1);
    __m128i mba = _mm_madd_epi16(ba, mmtx2);
    __m128i mrgba = _mm_add_epi32(mrg, mba);
    mrgba = _mm_srli_epi32(mrgba, 8);

    //__m128i lowerbound = _mm_setzero_si128();
    __m128i upperbound = _mm_set1_epi32(rgba[3]);
    mrgba = _mm_max_epi32(lowerbound, mrgba);
    mrgba = _mm_min_epi32(upperbound, mrgba);
    rgba[0] = mrgba.m128i_i32[3];
    rgba[1] = mrgba.m128i_i32[2];
    rgba[2] = mrgba.m128i_i32[1];
}
请问还可以如何进一步优化呢?

#11


既然map_rgba_sse会被调用非常非常多次,何不将其整合到一个大循环中,以免去不必要的函数调用。

#12


引用 11 楼 g_spider 的回复:
既然map_rgba_sse会被调用非常非常多次,何不将其整合到一个大循环中,以免去不必要的函数调用。

map_rgba_sse会被这样使用的
for (int i = 0; i < nTimes; i++)
{
    set_rgba(rgba);  
    map_rgba_sse(rgba);
}
如何整合?“以免去不必要的函数调用”,是指哪些函数调用?

#13


我还想问几个问题:
1、在编译的时候如何判断编译器是否支持SSE2 SSE4?
2、在运行的时候如何动态判断用户的CPU是否支持SSE2 SSE4?
3、类似于_mm_set1_epi32、_mm_max_epi32这样的函数能否在Linux下面使用?如何判断?

#14


打个比方:
for (int i = 0; i < nTimes; i++)
{
  set_rgba(rgba[***]);   
}
然后:
 map_rgba_sse(rgba[],n); 

//以上只是伪码。

#15


引用 10 楼 oyjian 的回复:
在我应用的过程中,矩阵mtx是固定不变的,而rgba是变化的,所以我改进了程序,效率是最开始那个版本的2倍!
//在外面初始化这些固定的数据,因为map_rgba_sse会被调用非常非常多次的!
const int *p = mtx[0];
__m128i mmtx1 = _mm_set_epi16(p[0], p[4], p[1], p[5], p[2], p[6], p[3], p[7……


  __m128i mrg = _mm_madd_epi16(rg, mmtx1);
  __m128i mba = _mm_madd_epi16(ba, mmtx2);

这可能运算溢出,结果不能保证正确。

#16


啊,如果保证mtx和rgba的元素都在0..255范围内也是可以的。

#17


引用 14 楼 g_spider 的回复:
打个比方:
for (int i = 0; i < nTimes; i++)
{
  set_rgba(rgba[***]);   
}
然后:
 map_rgba_sse(rgba[],n); 

//以上只是伪码。

这样改了之后有减少函数调用吗?
在map_rgba_sse(rgba[],n); 里面不也要进行一个for循环吗?
为什么这样改能够提高效率?

#18


感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。

#19


引用 18 楼 delphiguy 的回复:
感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。

rgba是uint,展开之后以int[4]传递进来,处理完毕之后再用<<与|压缩为uint。
如果直接传递uint进来的话,如何优化?

#20


int先强制转换成BYTE [],然后再操作。
首先把这个inline了罢。我昨天写着写着瞌睡就打岔了,但还有一个思路是,用不带跳转指令的几个指令代替掉最后的二个跳转指令,主要就是不需要跳转的求绝对值、判断减去一个数之后的正负,不过已经用_mm_min_epi32这种指令代替了,就不提这个了。
但外层是一个for,不管最终优化结构是把这个for的循环展开几次,也首先得修改成inline罢。
最后连续的最大、最小判断,我总觉得不对劲,也许是0的缘故.

#21


引用 19 楼 oyjian 的回复:
引用 18 楼 delphiguy 的回复:

感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。
rgba是uint,展开之后以int[4]传递进来,处理完毕之后再用<<与|压缩为uint。
如果直接传递uint进来的话,如何优化?


可以把连续的4个像素作为int128装入,处理完再一次写回。

#22


在大家指导下,我把处理一个像素的代码修改如下:
效率比当初好很多了。
谢谢大家的指导!
void map_rgba_sse(uint &argb)
{
    __m128i rg = _mm_set1_epi32((pb[2] << 16) | pb[1]);
    __m128i ba = _mm_set1_epi32((pb[0] << 16) | pb[3]);

    __m128i mrg = _mm_madd_epi16(rg, mmtx1);
    __m128i mba = _mm_madd_epi16(ba, mmtx2);
    __m128i mrgba = _mm_add_epi32(mrg, mba);
    mrgba = _mm_srli_epi32(mrgba, 8);
    
    __m128i upperbound = _mm_set1_epi32(pb[3]);
    mrgba = _mm_max_epi32(lowerbound, mrgba);
    mrgba = _mm_min_epi32(upperbound, mrgba);
    pb[2] = mrgba.m128i_i32[3];
    pb[1] = mrgba.m128i_i32[2];
    pb[0] = mrgba.m128i_i32[1];
}

#23


引用 21 楼 delphiguy 的回复:
引用 19 楼 oyjian 的回复:

引用 18 楼 delphiguy 的回复:

感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。
rgba是uint,展开之后以int[4]传递进来,处理完……

这也是一个很好的建议啊。Qt也是用这种方法来处理图片像素的。
不过我这里采用这种方法好像很难做到啊。

#24


不复杂,4个像素装入到一个xmm寄存器,然后每个像素用pshuflw、pshufhw、pshufd(或者用punpcklbw、pshufd也可以)组合成你需要的00grgrgr、00ababab形式,运算是一样的,只是重复做4次。
4个结果通过pslldq、por组合在一起,再做pshuf排好序就可以写回了。

#25


4个结果通过pslldq、por组合在一起,再做pshuf排好序就可以写回了。
================================================================
4个结果通过pslldq、por组合在一起,再做 pshufb排好序就可以写回了。

#26


引用 24 楼 delphiguy 的回复:
不复杂,4个像素装入到一个xmm寄存器,然后每个像素用pshuflw、pshufhw、pshufd(或者用punpcklbw、pshufd也可以)组合成你需要的00grgrgr、00ababab形式,运算是一样的,只是重复做4次。
4个结果通过pslldq、por组合在一起,再做pshuf排好序就可以写回了。

修改之后效率大幅度提升,代码如下:
void map_rgba_sse4_4(const int mtx[4][4], UINT *argb, int length)
{
    const int *p = mtx[0];
    __m128i *pSrc = reinterpret_cast<__m128i*>(&argb);
    __m128i *pEnd = pSrc + length / 4;
    const __m128i alphaMask = _mm_set1_epi32(0xff000000);
    const __m128i colorMask = _mm_set1_epi32(0x00ff00ff);
    __m128i mmtxred1 = _mm_set1_epi32((p[0] << 16) | p[8]);
    __m128i mmtxred2 = _mm_set1_epi32((p[12] << 16) | p[4]);
    __m128i mmtxg1 = _mm_set1_epi32((p[1] << 16) | p[9]);
    __m128i mmtxg2 = _mm_set1_epi32((p[13] << 16) | p[5]);
    __m128i mmtxb1 = _mm_set1_epi32((p[2] << 16) | p[10]);
    __m128i mmtxb2 = _mm_set1_epi32((p[14] << 16) | p[6]);
    for (; pSrc < pEnd; pSrc++)
    {
         const __m128i msrc = _mm_loadu_si128(pSrc);//有时候在这里崩溃
        const __m128i mrb = _mm_and_si128(msrc, colorMask);
        const __m128i mag = _mm_srli_epi16(msrc, 8);
        const __m128i malpha = _mm_and_si128(msrc, alphaMask);
        const __m128i upperbound = _mm_srli_epi32(malpha, 24);

        const __m128i mred1 = _mm_madd_epi16(mrb, mmtxred1);
        const __m128i mred2 = _mm_madd_epi16(mag, mmtxred2);
        __m128i mreds = _mm_add_epi32(mred1, mred2);
        mreds = _mm_srli_epi32(mreds, 8);
        mreds = _mm_max_epi32(lowerbound, mreds);
        mreds = _mm_min_epi32(upperbound, mreds);
        mreds = _mm_slli_epi32(mreds, 16);

        const __m128i mg1 = _mm_madd_epi16(mrb, mmtxg1);
        const __m128i mg2 = _mm_madd_epi16(mag, mmtxg2);
        __m128i mgs = _mm_add_epi32(mred1, mred2);
        mgs = _mm_srli_epi32(mgs, 8);
        mgs = _mm_max_epi32(lowerbound, mgs);
        mgs = _mm_min_epi32(upperbound, mgs);
        mgs = _mm_slli_epi32(mgs, 8);

        const __m128i mb1 = _mm_madd_epi16(mrb, mmtxb1);
        const __m128i mb2 = _mm_madd_epi16(mag, mmtxb2);
        __m128i mbs = _mm_add_epi32(mred1, mred2);
        mbs = _mm_srli_epi32(mbs, 8);
        mbs = _mm_max_epi32(lowerbound, mbs);
        mbs = _mm_min_epi32(upperbound, mbs);

        __m128i result = _mm_or_si128(malpha, mreds);
        result = _mm_or_si128(result, _mm_or_si128(mgs, mbs));

        _mm_storeu_si128(pSrc, result);
    }

    switch (length % 4)
    {
    case 3: map_rgba_sse4_unpack(mtx, argb[--length]);
    case 2: map_rgba_sse4_unpack(mtx, argb[--length]);
    case 1: map_rgba_sse4_unpack(mtx, argb[--length]);
    }
}
但是,有时候在const __m128i msrc = _mm_loadu_si128(pSrc);崩溃了。
这行代码的汇编为:
        const __m128i msrc = _mm_loadu_si128(pSrc);
00B12924  mov         eax,dword ptr [ebp-18h] 
00B12927  movdqu      xmm0,xmmword ptr [eax] 
00B1292B  movdqa      xmmword ptr [ebp-740h],xmm0 
00B12933  movdqa      xmm0,xmmword ptr [ebp-740h] 
00B1293B  movdqa      xmmword ptr [msrc],xmm0 

崩溃时,pSrc的值为:0x0024fff4(没有对齐16字节)
错误码为:First-chance exception at 0x00b12927 in sse.exe: 0xC0000005: Access violation reading location 0x00250000.
Unhandled exception at 0x77c36344 (ntdll.dll) in sse.exe: 0xC0000005: Access violation.
这样的错误,我前两天也遇到过,是因为字节没有对齐的原因。详细见: http://topic.csdn.net/u/20110527/20/8c9a7317-f968-49b3-a55a-ddb4f3a3a588.html?seed=1879088504&r=73573469#r_73573469
但现在这里的话,用的是_mm_loadu_si128函数,应该没有要求对齐16字节吧?但程序为什么会崩溃呢?

#27


它用了movdqa呀,movdqa要求地址16字节对齐。也许是_mm_loadu_si128的一个BUG。

另外,
00B12927 movdqu xmm0,xmmword ptr [eax] 
00B1292B movdqa xmmword ptr [ebp-740h],xmm0  // !!!
00B12933 movdqa xmm0,xmmword ptr [ebp-740h]  // !!!
00B1293B movdqa xmmword ptr [msrc],xmm0  

中间转储了一个也很令人费解,
直接:
movdqu xmm0,xmmword ptr [eax] 
movdq u xmmword ptr [msrc],xmm0  
不是很好吗。

#28


引用 27 楼 delphiguy 的回复:
它用了movdqa呀,movdqa要求地址16字节对齐。也许是_mm_loadu_si128的一个BUG。

另外,
00B12927 movdqu xmm0,xmmword ptr [eax] 
00B1292B movdqa xmmword ptr [ebp-740h],xmm0  // !!!
00B12933 movdqa xmm0,xmmword ptr [ebp-740h]……

崩溃的原因找到了!
__m128i *pSrc = reinterpret_cast<__m128i*>( &argb);
多了一个&,argb本身就是数组指针的。
这样看来,_mm_loadu_si128并没有bug。
因为movdqu xmm0,xmmword ptr [eax]已经用了u版本,我猜想,执行完这条指令之后xmm0是16字节对齐的,所以后面就可以使用movdqa了!

#29


即便地址[ebp-740h]是16字节对齐的,多做了两次movdqa也不可能快呀。

#1


它在VS2010编译出来的代码如何:
void map_rgba(int mtx[4][4], int rgba[4])
{
011F1000  push        ecx  
011F1001  mov         ecx,dword ptr [esp+8] 
    int r, g, b;
    r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
011F1005  mov         edx,dword ptr [eax+20h] 
011F1008  push        ebx  
011F1009  mov         ebx,dword ptr [ecx+4] 
011F100C  push        ebp  
011F100D  mov         ebp,dword ptr [ecx] 
011F100F  push        esi  
011F1010  mov         esi,dword ptr [ecx+8] 
011F1013  imul        edx,esi 
011F1016  push        edi  
011F1017  mov         edi,dword ptr [ecx+0Ch] 
011F101A  mov         ecx,dword ptr [eax+10h] 
011F101D  imul        ecx,ebx 
011F1020  add         ecx,edx 
011F1022  mov         edx,dword ptr [eax+30h] 
011F1025  imul        edx,edi 
011F1028  add         ecx,edx 
011F102A  mov         edx,dword ptr [eax] 
011F102C  imul        edx,ebp 
011F102F  add         ecx,edx 
    g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba[3] * mtx[3][1];
011F1031  mov         edx,dword ptr [eax+4] 
011F1034  imul        edx,ebp 
011F1037  mov         dword ptr [esp+10h],ebp 
011F103B  mov         ebp,dword ptr [eax+14h] 
011F103E  imul        ebp,ebx 
011F1041  add         edx,ebp 
011F1043  mov         ebp,dword ptr [eax+24h] 
011F1046  imul        ebp,esi 
011F1049  mov         esi,dword ptr [eax+34h] 
011F104C  imul        esi,edi 
011F104F  add         edx,ebp 
    b = rgba[0] * mtx[0][2] + rgba[1] * mtx[1][2] + rgba[2] * mtx[2][2] + rgba[3] * mtx[3][2];
011F1051  mov         ebp,dword ptr [eax+18h] 
011F1054  add         edx,esi 
011F1056  imul        ebp,ebx 
011F1059  mov         esi,dword ptr [eax+8] 
011F105C  imul        esi,dword ptr [esp+10h] 
011F1061  mov         ebx,dword ptr [eax+28h] 
011F1064  mov         eax,dword ptr [eax+38h] 
011F1067  add         esi,ebp 
011F1069  imul        eax,edi 
011F106C  mov         ebp,dword ptr [esp+18h] 
011F1070  imul        ebx,dword ptr [ebp+8] 
011F1074  add         esi,ebx 
011F1076  add         esi,eax 
    rgba[0] = min(max(0, r >> 8), rgba[3] );
011F1078  sar         ecx,8 
011F107B  xor         eax,eax 
011F107D  test        ecx,ecx 
011F107F  setl        al   
011F1082  dec         eax  
011F1083  and         eax,ecx 
011F1085  cmp         eax,edi 
011F1087  jge         map_rgba+95h (11F1095h) 
011F1089  xor         eax,eax 
011F108B  test        ecx,ecx 
011F108D  setl        al   
011F1090  dec         eax  
011F1091  and         ecx,eax 
011F1093  jmp         map_rgba+97h (11F1097h) 
011F1095  mov         ecx,edi 
    rgba[1] = min(max(0, g >> 8), rgba[3]);
011F1097  sar         edx,8 
011F109A  mov         eax,ebp 
011F109C  mov         dword ptr [eax],ecx 
011F109E  xor         ecx,ecx 
011F10A0  test        edx,edx 
011F10A2  setl        cl   
011F10A5  dec         ecx  
011F10A6  and         ecx,edx 
011F10A8  cmp         ecx,edi 
011F10AA  jge         map_rgba+0B8h (11F10B8h) 
011F10AC  xor         ecx,ecx 
011F10AE  test        edx,edx 
011F10B0  setl        cl   
011F10B3  dec         ecx  
011F10B4  and         ecx,edx 
011F10B6  jmp         map_rgba+0BAh (11F10BAh) 
011F10B8  mov         ecx,edi 
011F10BA  mov         dword ptr [eax+4],ecx 
    rgba[2] = min(max(0, b >> 8), rgba[3]);
011F10BD  sar         esi,8 
011F10C0  xor         ecx,ecx 
011F10C2  test        esi,esi 
011F10C4  setl        cl   
011F10C7  dec         ecx  
011F10C8  and         ecx,esi 
011F10CA  cmp         ecx,edi 
011F10CC  jge         map_rgba+0E1h (11F10E1h) 
011F10CE  xor         ecx,ecx 
011F10D0  test        esi,esi 
011F10D2  setl        cl   
011F10D5  pop         edi  
011F10D6  dec         ecx  
011F10D7  and         ecx,esi 
011F10D9  pop         esi  
011F10DA  pop         ebp  
011F10DB  mov         dword ptr [eax+8],ecx 
011F10DE  pop         ebx  

#2


hao

#3


可以用PMADDWD(SSE2)、PHADDD(SSSE3)指令。

#4


引用 3 楼 delphiguy 的回复:
可以用PMADDWD(SSE2)、PHADDD(SSSE3)指令。

能否具体点?

#5


建议先弄明白这些指令的含意,然后才知道往哪使劲。

#6


已经很底层了,指令级上估计没什么可优化了,算法应该还可以在优化一点。
比方说求 r g b中的rgba[x]数组,没有必要每次都这样在引用一下把,那么此时必然多一个mov指令,你可以将此值放到一个reg中,下次直接使用。
还有那个矩阵也一样,维数已定,你自己通过偏移add 应该比编译器生成的要快把。
当然了这些代码就必须由你来实现了。关键的是,就算你这样做了,估计也没有什么质的改变。除非还有更好的方法。
一些个人看法。。。

#7


PMADDWD可以一条指令做8个乘法,4个加法。

  r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
  g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba[3] * mtx[3][1];
  b = rgba[0] * mtx[0][2] + rgba[1] * mtx[1][2] + rgba[2] * mtx[2][2] + rgba[3] * mtx[3][2];

这种运算,可以用4条PMADDWD和3条PADDD指令实现(不包括装入数据的指令)。

#8


引用 7 楼 delphiguy 的回复:
PMADDWD可以一条指令做8个乘法,4个加法。

  r = rgba[0] * mtx[0][0] + rgba[1] * mtx[1][0] + rgba[2] * mtx[2][0] + rgba[3] * mtx[3][0];
  g = rgba[0] * mtx[0][1] + rgba[1] * mtx[1][1] + rgba[2] * mtx[2][1] + rgba……


已经使用这两条指令进行优化了,但是效率并没有明显提高。请问还可以如何进一步优化呢?
#include "emmintrin.h"
void map_rgba_sse(const int mtx[4][4], int rgba[4])
{
    const int *p = mtx[0];
    __m128i mmtx1 = _mm_set_epi16(p[0], p[4], p[1], p[5], p[2], p[6], p[3], p[7]);
    __m128i mmtx2 = _mm_set_epi16(p[8], p[12], p[9], p[13], p[10], p[14], p[11], p[15]);

    __m128i rg = _mm_set1_epi32((rgba[0] << 16) | rgba[1]);
    __m128i ba = _mm_set1_epi32((rgba[2] << 16) | rgba[3]);

    __m128i mrg = _mm_madd_epi16(rg, mmtx1);
    __m128i mba = _mm_madd_epi16(ba, mmtx2);
    __m128i mrgba = _mm_add_epi32(mrg, mba);
    mrgba = _mm_srli_epi32(mrgba, 8);

    rgba[0] = min(max(0, mrgba.m128i_i32[3]), rgba[3] );
    rgba[1] = min(max(0, mrgba.m128i_i32[2]), rgba[3]);
    rgba[2] = min(max(0, mrgba.m128i_i32[1]), rgba[3]);
}
min(max(0, mrgba.m128i_i32[3]), rgba[3] );这些代码又如何优化呢?

#9


优化的非常不错,不看算法我觉得全用协处理器代替处理器运算在现在已经可以处理绝多数情况了。
只是最后一片的跳转会对性能大打折扣,还没有那种超惊人的处理器论文能设计一套摆脱预跳转失误导致的性能损耗呢,特别是预处理了跳转目标后要等待慢腾腾的协处理器来返回自己有没有猜错的处理器更加不喜欢跳转指令了。

暂时想到的二个思路,首先可以继续用指令MINPD

然后呢,可能已经不实用了,记得有一代CPU的介绍里据说已经把预跳转出错处理到20%性能损耗内。
但思路是用赋值代替if,用符号位(是否为负数的标志bit)作为数组的index.有点累,代码明天写试试。

#10


在我应用的过程中,矩阵mtx是固定不变的,而rgba是变化的,所以我改进了程序,效率是最开始那个版本的2倍!
//在外面初始化这些固定的数据,因为map_rgba_sse会被调用非常非常多次的!
const int *p = mtx[0];
__m128i mmtx1 = _mm_set_epi16(p[0], p[4], p[1], p[5], p[2], p[6], p[3], p[7]);
__m128i mmtx2 = _mm_set_epi16(p[8], p[12], p[9], p[13], p[10], p[14], p[11], p[15]);
__m128i lowerbound = _mm_setzero_si128();

void map_rgba_sse(int rgba[4])
{
    __m128i rg = _mm_set1_epi32((rgba[0] << 16) | rgba[1]);
    __m128i ba = _mm_set1_epi32((rgba[2] << 16) | rgba[3]);

    __m128i mrg = _mm_madd_epi16(rg, mmtx1);
    __m128i mba = _mm_madd_epi16(ba, mmtx2);
    __m128i mrgba = _mm_add_epi32(mrg, mba);
    mrgba = _mm_srli_epi32(mrgba, 8);

    //__m128i lowerbound = _mm_setzero_si128();
    __m128i upperbound = _mm_set1_epi32(rgba[3]);
    mrgba = _mm_max_epi32(lowerbound, mrgba);
    mrgba = _mm_min_epi32(upperbound, mrgba);
    rgba[0] = mrgba.m128i_i32[3];
    rgba[1] = mrgba.m128i_i32[2];
    rgba[2] = mrgba.m128i_i32[1];
}
请问还可以如何进一步优化呢?

#11


既然map_rgba_sse会被调用非常非常多次,何不将其整合到一个大循环中,以免去不必要的函数调用。

#12


引用 11 楼 g_spider 的回复:
既然map_rgba_sse会被调用非常非常多次,何不将其整合到一个大循环中,以免去不必要的函数调用。

map_rgba_sse会被这样使用的
for (int i = 0; i < nTimes; i++)
{
    set_rgba(rgba);  
    map_rgba_sse(rgba);
}
如何整合?“以免去不必要的函数调用”,是指哪些函数调用?

#13


我还想问几个问题:
1、在编译的时候如何判断编译器是否支持SSE2 SSE4?
2、在运行的时候如何动态判断用户的CPU是否支持SSE2 SSE4?
3、类似于_mm_set1_epi32、_mm_max_epi32这样的函数能否在Linux下面使用?如何判断?

#14


打个比方:
for (int i = 0; i < nTimes; i++)
{
  set_rgba(rgba[***]);   
}
然后:
 map_rgba_sse(rgba[],n); 

//以上只是伪码。

#15


引用 10 楼 oyjian 的回复:
在我应用的过程中,矩阵mtx是固定不变的,而rgba是变化的,所以我改进了程序,效率是最开始那个版本的2倍!
//在外面初始化这些固定的数据,因为map_rgba_sse会被调用非常非常多次的!
const int *p = mtx[0];
__m128i mmtx1 = _mm_set_epi16(p[0], p[4], p[1], p[5], p[2], p[6], p[3], p[7……


  __m128i mrg = _mm_madd_epi16(rg, mmtx1);
  __m128i mba = _mm_madd_epi16(ba, mmtx2);

这可能运算溢出,结果不能保证正确。

#16


啊,如果保证mtx和rgba的元素都在0..255范围内也是可以的。

#17


引用 14 楼 g_spider 的回复:
打个比方:
for (int i = 0; i < nTimes; i++)
{
  set_rgba(rgba[***]);   
}
然后:
 map_rgba_sse(rgba[],n); 

//以上只是伪码。

这样改了之后有减少函数调用吗?
在map_rgba_sse(rgba[],n); 里面不也要进行一个for循环吗?
为什么这样改能够提高效率?

#18


感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。

#19


引用 18 楼 delphiguy 的回复:
感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。

rgba是uint,展开之后以int[4]传递进来,处理完毕之后再用<<与|压缩为uint。
如果直接传递uint进来的话,如何优化?

#20


int先强制转换成BYTE [],然后再操作。
首先把这个inline了罢。我昨天写着写着瞌睡就打岔了,但还有一个思路是,用不带跳转指令的几个指令代替掉最后的二个跳转指令,主要就是不需要跳转的求绝对值、判断减去一个数之后的正负,不过已经用_mm_min_epi32这种指令代替了,就不提这个了。
但外层是一个for,不管最终优化结构是把这个for的循环展开几次,也首先得修改成inline罢。
最后连续的最大、最小判断,我总觉得不对劲,也许是0的缘故.

#21


引用 19 楼 oyjian 的回复:
引用 18 楼 delphiguy 的回复:

感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。
rgba是uint,展开之后以int[4]传递进来,处理完毕之后再用<<与|压缩为uint。
如果直接传递uint进来的话,如何优化?


可以把连续的4个像素作为int128装入,处理完再一次写回。

#22


在大家指导下,我把处理一个像素的代码修改如下:
效率比当初好很多了。
谢谢大家的指导!
void map_rgba_sse(uint &argb)
{
    __m128i rg = _mm_set1_epi32((pb[2] << 16) | pb[1]);
    __m128i ba = _mm_set1_epi32((pb[0] << 16) | pb[3]);

    __m128i mrg = _mm_madd_epi16(rg, mmtx1);
    __m128i mba = _mm_madd_epi16(ba, mmtx2);
    __m128i mrgba = _mm_add_epi32(mrg, mba);
    mrgba = _mm_srli_epi32(mrgba, 8);
    
    __m128i upperbound = _mm_set1_epi32(pb[3]);
    mrgba = _mm_max_epi32(lowerbound, mrgba);
    mrgba = _mm_min_epi32(upperbound, mrgba);
    pb[2] = mrgba.m128i_i32[3];
    pb[1] = mrgba.m128i_i32[2];
    pb[0] = mrgba.m128i_i32[1];
}

#23


引用 21 楼 delphiguy 的回复:
引用 19 楼 oyjian 的回复:

引用 18 楼 delphiguy 的回复:

感觉max(0, r >> 8)这种运算是不需要的,rgb都是4对8位数乘积之和,作为int型数据永远不会成为负数。

另外,rgba应该本来是int型吧,你作为int[4]处理,先要展开一次,处理完还要压缩一次,这也是很浪费时间的。
rgba是uint,展开之后以int[4]传递进来,处理完……

这也是一个很好的建议啊。Qt也是用这种方法来处理图片像素的。
不过我这里采用这种方法好像很难做到啊。

#24


不复杂,4个像素装入到一个xmm寄存器,然后每个像素用pshuflw、pshufhw、pshufd(或者用punpcklbw、pshufd也可以)组合成你需要的00grgrgr、00ababab形式,运算是一样的,只是重复做4次。
4个结果通过pslldq、por组合在一起,再做pshuf排好序就可以写回了。

#25


4个结果通过pslldq、por组合在一起,再做pshuf排好序就可以写回了。
================================================================
4个结果通过pslldq、por组合在一起,再做 pshufb排好序就可以写回了。

#26


引用 24 楼 delphiguy 的回复:
不复杂,4个像素装入到一个xmm寄存器,然后每个像素用pshuflw、pshufhw、pshufd(或者用punpcklbw、pshufd也可以)组合成你需要的00grgrgr、00ababab形式,运算是一样的,只是重复做4次。
4个结果通过pslldq、por组合在一起,再做pshuf排好序就可以写回了。

修改之后效率大幅度提升,代码如下:
void map_rgba_sse4_4(const int mtx[4][4], UINT *argb, int length)
{
    const int *p = mtx[0];
    __m128i *pSrc = reinterpret_cast<__m128i*>(&argb);
    __m128i *pEnd = pSrc + length / 4;
    const __m128i alphaMask = _mm_set1_epi32(0xff000000);
    const __m128i colorMask = _mm_set1_epi32(0x00ff00ff);
    __m128i mmtxred1 = _mm_set1_epi32((p[0] << 16) | p[8]);
    __m128i mmtxred2 = _mm_set1_epi32((p[12] << 16) | p[4]);
    __m128i mmtxg1 = _mm_set1_epi32((p[1] << 16) | p[9]);
    __m128i mmtxg2 = _mm_set1_epi32((p[13] << 16) | p[5]);
    __m128i mmtxb1 = _mm_set1_epi32((p[2] << 16) | p[10]);
    __m128i mmtxb2 = _mm_set1_epi32((p[14] << 16) | p[6]);
    for (; pSrc < pEnd; pSrc++)
    {
         const __m128i msrc = _mm_loadu_si128(pSrc);//有时候在这里崩溃
        const __m128i mrb = _mm_and_si128(msrc, colorMask);
        const __m128i mag = _mm_srli_epi16(msrc, 8);
        const __m128i malpha = _mm_and_si128(msrc, alphaMask);
        const __m128i upperbound = _mm_srli_epi32(malpha, 24);

        const __m128i mred1 = _mm_madd_epi16(mrb, mmtxred1);
        const __m128i mred2 = _mm_madd_epi16(mag, mmtxred2);
        __m128i mreds = _mm_add_epi32(mred1, mred2);
        mreds = _mm_srli_epi32(mreds, 8);
        mreds = _mm_max_epi32(lowerbound, mreds);
        mreds = _mm_min_epi32(upperbound, mreds);
        mreds = _mm_slli_epi32(mreds, 16);

        const __m128i mg1 = _mm_madd_epi16(mrb, mmtxg1);
        const __m128i mg2 = _mm_madd_epi16(mag, mmtxg2);
        __m128i mgs = _mm_add_epi32(mred1, mred2);
        mgs = _mm_srli_epi32(mgs, 8);
        mgs = _mm_max_epi32(lowerbound, mgs);
        mgs = _mm_min_epi32(upperbound, mgs);
        mgs = _mm_slli_epi32(mgs, 8);

        const __m128i mb1 = _mm_madd_epi16(mrb, mmtxb1);
        const __m128i mb2 = _mm_madd_epi16(mag, mmtxb2);
        __m128i mbs = _mm_add_epi32(mred1, mred2);
        mbs = _mm_srli_epi32(mbs, 8);
        mbs = _mm_max_epi32(lowerbound, mbs);
        mbs = _mm_min_epi32(upperbound, mbs);

        __m128i result = _mm_or_si128(malpha, mreds);
        result = _mm_or_si128(result, _mm_or_si128(mgs, mbs));

        _mm_storeu_si128(pSrc, result);
    }

    switch (length % 4)
    {
    case 3: map_rgba_sse4_unpack(mtx, argb[--length]);
    case 2: map_rgba_sse4_unpack(mtx, argb[--length]);
    case 1: map_rgba_sse4_unpack(mtx, argb[--length]);
    }
}
但是,有时候在const __m128i msrc = _mm_loadu_si128(pSrc);崩溃了。
这行代码的汇编为:
        const __m128i msrc = _mm_loadu_si128(pSrc);
00B12924  mov         eax,dword ptr [ebp-18h] 
00B12927  movdqu      xmm0,xmmword ptr [eax] 
00B1292B  movdqa      xmmword ptr [ebp-740h],xmm0 
00B12933  movdqa      xmm0,xmmword ptr [ebp-740h] 
00B1293B  movdqa      xmmword ptr [msrc],xmm0 

崩溃时,pSrc的值为:0x0024fff4(没有对齐16字节)
错误码为:First-chance exception at 0x00b12927 in sse.exe: 0xC0000005: Access violation reading location 0x00250000.
Unhandled exception at 0x77c36344 (ntdll.dll) in sse.exe: 0xC0000005: Access violation.
这样的错误,我前两天也遇到过,是因为字节没有对齐的原因。详细见: http://topic.csdn.net/u/20110527/20/8c9a7317-f968-49b3-a55a-ddb4f3a3a588.html?seed=1879088504&r=73573469#r_73573469
但现在这里的话,用的是_mm_loadu_si128函数,应该没有要求对齐16字节吧?但程序为什么会崩溃呢?

#27


它用了movdqa呀,movdqa要求地址16字节对齐。也许是_mm_loadu_si128的一个BUG。

另外,
00B12927 movdqu xmm0,xmmword ptr [eax] 
00B1292B movdqa xmmword ptr [ebp-740h],xmm0  // !!!
00B12933 movdqa xmm0,xmmword ptr [ebp-740h]  // !!!
00B1293B movdqa xmmword ptr [msrc],xmm0  

中间转储了一个也很令人费解,
直接:
movdqu xmm0,xmmword ptr [eax] 
movdq u xmmword ptr [msrc],xmm0  
不是很好吗。

#28


引用 27 楼 delphiguy 的回复:
它用了movdqa呀,movdqa要求地址16字节对齐。也许是_mm_loadu_si128的一个BUG。

另外,
00B12927 movdqu xmm0,xmmword ptr [eax] 
00B1292B movdqa xmmword ptr [ebp-740h],xmm0  // !!!
00B12933 movdqa xmm0,xmmword ptr [ebp-740h]……

崩溃的原因找到了!
__m128i *pSrc = reinterpret_cast<__m128i*>( &argb);
多了一个&,argb本身就是数组指针的。
这样看来,_mm_loadu_si128并没有bug。
因为movdqu xmm0,xmmword ptr [eax]已经用了u版本,我猜想,执行完这条指令之后xmm0是16字节对齐的,所以后面就可以使用movdqa了!

#29


即便地址[ebp-740h]是16字节对齐的,多做了两次movdqa也不可能快呀。