[讨论]复数乘法SSE,SSE2优化

[shines77]

float版的SSE优化来源于AMD官方，double版的SSE2优化

源代码（包含初级版的FFT大数乘法，源代码已经测试exe）：
http://www.77studio.net/files/cmplx_test.rar

SSE2版本没有经过测试（我的CPU不支持），欢迎测试，讨论

下面是float版的复数乘法，一次同时计算4个：

__declspec(naked)
void cmplx_multiply_sse(float *x, float *y, int num_cmplx_elem, float *prod)
{
__asm
{
;
; TO ASSEMBLE INTO *.obj DO THE FOLLOWING:
; ml.exe -coff -c cmplx_multiply_sse.asm
;
;.586
;.K3D
;.XMM
;_TEXT SEGMENT
;PUBLIC _cmplx_multiply_sse
;_cmplx_multiply_sse PROC NEAR
;==============================================================================
; INSTRUCTIONS BELOW SAVE THE REGISTER STATE WITH WHICH THIS ROUTINE WAS ENTERED
; REGISTERS (EAX, ECX, EDX ARE CONSIDERED VOLATILE AND ASSUMED TO BE CHANGED)
; WHILE THE REGISTERS BELOW MUST BE PRESERVED IF THE USER IS CHANGING THEM
push ebp
mov ebp, esp
;==============================================================================
; parameters passed into routine:
; [ebp+8] = ->x
; [ebp+12] = ->y
; [ebp+16] = num_cmplx_elem
; [ebp+20] = ->prod
;==============================================================================
push ebx ; preserve contents in ebx,esi, and edi on stack
push esi ;
push edi ;
;===============================================================================
; THE CODE BELOW PUTS THE FLOATING POINT SIGN MASK
; [800000000000000800000000000000h]
; TO FLIP THE SIGN OF PACKED SINGLE PRECISION NUMBERS BY USING XORPS
;==============================================================================
mov eax, esp ; Copy stack pointer into EAX.
mov ebx, 16
sub esp, 32 ; Subtract 32 bytes from stack pointer.
and eax, 15 ; AND old stack pointer address with 15 to
; determine # of bytes the address is past a
; 16-byte-aligned address.
sub ebx, eax ; EBX = # of bytes above ESP to next
; 16-byte-aligned address
mov edi, 0h ; EDI = 00000000h
mov esi, 80000000h ; ESI = 80000000h
shr ebx, 2 ; EBX = # of DWORDs past 16-byte-aligned address
mov [esp+4*ebx+12], esi ; Move into address esp+4*ebx the single-precision
mov [esp+4*ebx+8], edi ; floating-point sign mask.
mov [esp+4*ebx+4], esi
mov [esp+4*ebx], edi
;==============================================================================
; THE 4 ASM LINES BELOW LOAD THE FUNCTION s ARGUMENTS INTO GENERAL-PURPOSE
; REGISTERS (GPRS)
; esi = address of array "x"
; edi = address of array "y"
; ecx = # of cmplx products to compute
; eax = address of product to which results are stored
;==============================================================================
mov esi, [ebp+8] ; esi = ->x
mov edi, [ebp+12] ; edi = ->y
mov ecx, [ebp+16] ; ecx = num_cmplx_elem
mov eax, [ebp+20] ; eax = ->prod
;==============================================================================
; THE 6 ASM LINES BELOW OFFSET THE ADDRESS TO THE ARRAYS x[] AND y[] SUCH
; THAT THEY CAN BE ACCESSED IN THE MOST EFFICIENT MANNER AS ILLUSTRATED
; BELOW IN THE LOOP mult4cmplxnum_loop WITH THE MINIMUM NUMBER OF
; ADDRESS INCREMENTS
;==============================================================================
mov edx, ecx ; edx = num_cmplx_elem
neg ecx ; ecx = -num_cmplx_elem
shl edx, 3 ; edx = 8 * num_cmplx_elem = # bytes in x[] and y[] to multiply
add esi, edx ; esi = -> to last element of x[] to multiply
add edi, edx ; edi = -> to last element of y[] to multiply
add eax, edx ; eax = -> end of prod[] to calculate
;==============================================================================
; THIS LOOP MULTIPLIES 4 COMPLEX #s FROM "x[]" UPON 4 COMPLEX #s FROM "y[]"
; AND RETURNS THE PRODUCT IN "prod[]".
;==============================================================================
ALIGN 16 ; Align address of loop to a 16-byte boundary.
eight_cmplx_prod_loop:
///*
movaps xmm0, [esi+ecx*8] ; xmm0=[x1i,x1r,x0i,x0r]
movaps xmm1, [esi+ecx*8+16] ; xmm1=[x3i,x3r,x2i,x2r]
movaps xmm4, [edi+ecx*8] ; xmm4=[y1i,y1r,y0i,y0r]
movaps xmm5, [edi+ecx*8+16] ; xmm5=[y3i,y3r,y2i,y2r]
movaps xmm2, xmm0 ; xmm2=[x1i,x1r,x0i,x0r]
movaps xmm3, xmm1 ; xmm3=[x3i,x3r,x2i,x2r]
movaps xmm6, xmm4 ; xmm6=[y1i,y1r,y0i,y0r]
movaps xmm7, xmm5 ; xmm7=[y3i,y3r,y2i,y2r]
//*/

/*
// 数据无法对齐时可以用这个
movups xmm0, [esi+ecx*8] ; xmm0=[x1i,x1r,x0i,x0r]
movups xmm1, [esi+ecx*8+16] ; xmm1=[x3i,x3r,x2i,x2r]
movups xmm4, [edi+ecx*8] ; xmm4=[y1i,y1r,y0i,y0r]
movups xmm5, [edi+ecx*8+16] ; xmm5=[y3i,y3r,y2i,y2r]
movups xmm2, xmm0 ; xmm2=[x1i,x1r,x0i,x0r]
movups xmm3, xmm1 ; xmm3=[x3i,x3r,x2i,x2r]
movups xmm6, xmm4 ; xmm6=[y1i,y1r,y0i,y0r]
movups xmm7, xmm5 ; xmm7=[y3i,y3r,y2i,y2r]
//*/

shufps xmm0, xmm0, 10100000b ; xmm0=[x1r,x1r,x0r,x0r]
shufps xmm1, xmm1, 10100000b ; xmm1=[x3r,x3r,x2r,x2r]
shufps xmm2, xmm2, 11110101b ; xmm2=[x1i,x1i,x0i,x0i]
shufps xmm3, xmm3, 11110101b ; xmm3=[x3i,x3i,x2i,x2i]
xorps xmm6, [esp+4*ebx] ; xmm6=[-y1i,y1r,-y0i,y0r]
xorps xmm7, [esp+4*ebx] ; xmm7=[-y3i,y3r,-y2i,y2r]
mulps xmm0, xmm4 ; xmm0=[x1r*y1i,x1r*y1r,x0r*y0i,x0r*y0r]
mulps xmm1, xmm5 ; xmm1=[x3r*y3i,x3r*y3r,x2r*y2i,x2r*y2r]
shufps xmm7, xmm7, 10110001b ; xmm7=[y3r,-y3i,y2r,-y2i]
mulps xmm2, xmm6 ; xmm2=[x1i*y1r,-x1i*y1i,x0i*y0r,-x0i*y0i]
mulps xmm3, xmm7 ; xmm3=[x3i*y3r,-x3i*y3i,x2i*y2r,-x2i*y2i]
addps xmm0, xmm2 ; xmm0=[x1r*y1i+x1i*y1r,x1r*y1r-x1i*y1i,
; x0r*y0i+x0i*y0r,x0r*y0r-x0i*y0i]
addps xmm1, xmm3 ; xmm1=[x3r*y3i+x3i*y3r,x3r*y3r-x3i*y3i,
; x2r*y2i+x2i*y2r,x2r*y2r-x2i*y2i]
movntps [eax+ecx*8], xmm0 ; Stream XMM0 and XMM1 to representative
movntps [eax+ecx*8+16], xmm1 ; memory address of prod[].
add ecx, 4 ; ECX = ECX + 4, 一次计算了4个(float)复数
jnz eight_cmplx_prod_loop
sfence ; Finish all memory writes.
;==============================================================================
; INSTRUCTIONS BELOW RESTORE THE REGISTER STATE WITH WHICH THIS ROUTINE WAS
; ENTERED
; REGISTERS EAX, ECX, AND EDX ARE CONSIDERED VOLATILE AND ASSUMED TO BE CHANGED
; WHILE THE REGISTERS BELOW MUST BE PRESERVED IF THE USER IS CHANGING THEM
add esp, 32
pop edi
pop esi
pop ebx
mov esp, ebp
pop ebp
;==============================================================================
ret
}
}

[shines77]

SSE2版的double复数乘法，是我自己改的，没有经过测试，欢迎指正

要完整代码，请到以下地址去下载
http://www.77studio.net/files/cmplx_src.rar （上面地址是错的！）

__declspec(naked)
void cmplx_multiply_sse(double *x, double *y, int num_cmplx_elem, double *prod)
{
__asm
{
;==============================================================================
; INSTRUCTIONS BELOW SAVE THE REGISTER STATE WITH WHICH THIS ROUTINE WAS ENTERED
; REGISTERS (EAX, ECX, EDX ARE CONSIDERED VOLATILE AND ASSUMED TO BE CHANGED)
; WHILE THE REGISTERS BELOW MUST BE PRESERVED IF THE USER IS CHANGING THEM

push ebp
mov ebp, esp

;==============================================================================
; parameters passed into routine:
; [ebp+8] = ->x
; [ebp+12] = ->y
; [ebp+16] = num_cmplx_elem
; [ebp+20] = ->prod
;==============================================================================

push ebx
push esi
push edi

;===============================================================================
; THE CODE BELOW PUTS THE FLOATING POINT SIGN MASK
; [800000000000000000000000000000h]
; TO FLIP THE SIGN OF PACKED DOUBLE PRECISION NUMBERS BY USING XORPS OR XORPD
;===============================================================================

mov eax, esp ; Copy stack pointer into EAX.
mov ebx, 16
sub esp, 32 ; Subtract 32 bytes from stack pointer.
and eax, 15 ; AND old stack pointer address with 15 to
; determine # of bytes the address is past a
; 16-byte-aligned address.
sub ebx, eax ; EBX = # of bytes above ESP to next
; 16-byte-aligned address
mov edi, 0h ; EDI = 00000000h
mov esi, 80000000h ; ESI = 80000000h
shr ebx, 2 ; EBX = # of DWORDs past 16-byte-aligned address
mov [esp+4*ebx+12], esi ; Move into address esp+4*ebx the double-precision
mov [esp+4*ebx+8], edi ; floating-point sign mask.
mov [esp+4*ebx+4], edi
mov [esp+4*ebx], edi

;==============================================================================
; THE 4 ASM LINES BELOW LOAD THE FUNCTION s ARGUMENTS INTO GENERAL-PURPOSE
; REGISTERS (GPRS)
; esi = address of array "x"
; edi = address of array "y"
; ecx = # of cmplx products to compute
; eax = address of product to which results are stored
;==============================================================================

mov esi, [ebp+8] ; esi = ->x
mov edi, [ebp+12] ; edi = ->y
mov ecx, [ebp+16] ; ecx = num_cmplx_elem
mov eax, [ebp+20] ; eax = ->prod

;==============================================================================
; THE 6 ASM LINES BELOW OFFSET THE ADDRESS TO THE ARRAYS x[] AND y[] SUCH
; THAT THEY CAN BE ACCESSED IN THE MOST EFFICIENT MANNER AS ILLUSTRATED
; BELOW IN THE LOOP mult4cmplxnum_loop WITH THE MINIMUM NUMBER OF
; ADDRESS INCREMENTS
;==============================================================================

mov edx, ecx ; edx = num_cmplx_elem
neg ecx ; ecx = -num_cmplx_elem
shl edx, 4 ; edx = 16 * num_cmplx_elem = # bytes in x[] and y[] to multiply
shl ecx, 2 ; 注：因为不支持 [esi+ecx*16]，只能用 [esi+ecx*8]，所以把ecx放大一倍
add esi, edx ; esi = -> to last element of x[] to multiply
add edi, edx ; edi = -> to last element of y[] to multiply
add eax, edx ; eax = -> end of prod[] to calculate

;==============================================================================
; THIS LOOP MULTIPLIES 4 COMPLEX #s FROM "x[]" UPON 4 COMPLEX #s FROM "y[]"
; AND RETURNS THE PRODUCT IN "prod[]".
;==============================================================================
ALIGN 16 ; Align address of loop to a 16-byte boundary.

two_cmplx_prod_loop:

// 数据无法16字节对齐的时候用这个 ^_^
//movupd xmm0, [esi+ecx*8] ; xmm0=[x0i,x0r]
//movupd xmm1, [esi+ecx*8+16] ; xmm0=[x1i,x1r]
//movupd xmm4, [edi+ecx*8] ; xmm0=[y0i,y0r]
//movupd xmm5, [edi+ecx*8+16] ; xmm0=[y1i,y1r]

movapd xmm0, [esi+ecx*8] ; xmm0=[x0i,x0r]
movapd xmm1, [esi+ecx*8+16] ; xmm1=[x1i,x1r]
movapd xmm4, [edi+ecx*8] ; xmm4=[y0i,y0r]
movapd xmm5, [edi+ecx*8+16] ; xmm5=[y1i,y1r]

movapd xmm2, xmm0 ; xmm2=[x0i,x0r]
movapd xmm3, xmm1 ; xmm3=[x1i,x1r]
movapd xmm6, xmm4 ; xmm6=[y0i,y0r]
movapd xmm7, xmm5 ; xmm7=[y1i,y1r]

/*
// SSE版本的变换，没有测试环境，不知道哪个快
shufps xmm0, xmm0, 01000100b ; xmm0=[x0r,x0r]
shufps xmm1, xmm1, 01000100b ; xmm1=[x1r,x1r]
shufps xmm2, xmm2, 11101110b ; xmm2=[x0i,x0i]
shufps xmm3, xmm3, 11101110b ; xmm3=[x1i,x1i]
//*/

// SSE2版本的变换
shufpd xmm0, xmm0, 00000000b ; xmm0=[x0r,x0r]
shufpd xmm1, xmm1, 00000000b ; xmm1=[x1r,x1r]
shufpd xmm2, xmm2, 00000011b ; xmm2=[x0i,x0i]
shufpd xmm3, xmm3, 00000011b ; xmm3=[x1i,x1i]

/*
// SSE版本，在这里的实现效果应该和xorpd是一样的
// [esp+4*ebx] = [800000000000000000000000000000h]
xorps xmm6, [esp+4*ebx] ; xmm6=[-y0i,y0r]
xorps xmm7, [esp+4*ebx] ; xmm7=[-y1i,y1r]
//*/

// SSE2版本的变号
// [esp+4*ebx] = [800000000000000000000000000000h]
xorpd xmm6, [esp+4*ebx] ; xmm6=[-y0i,y0r]
xorpd xmm7, [esp+4*ebx] ; xmm7=[-y1i,y1r]

mulpd xmm0, xmm4 ; xmm0=[x0r*y0i,x0r*y0r]
mulpd xmm1, xmm5 ; xmm1=[x1r*y1i,x1r*y1r]

/*
// SSE版本的变换
shufps xmm6, xmm6, 01001110b ; xmm6=[y0r,-y0i]
shufps xmm7, xmm7, 01001110b ; xmm7=[y1r,-y1i]
//*/

// SSE2版本的变换
shufpd xmm6, xmm6, 00000001b ; xmm6=[y0r,-y0i]
shufpd xmm7, xmm7, 00000001b ; xmm7=[y1r,-y1i]

mulpd xmm2, xmm6 ; xmm2=[x0i*y0r,-x0i*y0i]
mulpd xmm3, xmm7 ; xmm3=[x1i*y1r,-x1i*y1i]

addpd xmm0, xmm2 ; xmm0=[x0r*y0i+x0i*y0r,x0r*y0r-x0i*y0i]
addpd xmm1, xmm3 ; xmm1=[x1r*y1i+x1i*y1r,x1r*y1r-x1i*y1i]

/*
// 同样的，这个是SSE版本，和下面应该是等效的
movntps [eax+ecx*8], xmm0 ; Stream XMM0 and XMM1 to representative
movntps [eax+ecx*8+16], xmm1 ; memory address of prod[].
//*/

movntpd [eax+ecx*8], xmm0 ; Stream XMM0 and XMM1 to representative
movntpd [eax+ecx*8+16], xmm1 ; memory address of prod[].

add ecx, 4 ; ECX = ECX + 4, 一次计算了2个(double)复数
; ? 为什么是加4，而不是加2，因为前面 ecx 放大了一倍
jnz two_cmplx_prod_loop
sfence ; Finish all memory writes.

;==============================================================================
; INSTRUCTIONS BELOW RESTORE THE REGISTER STATE WITH WHICH THIS ROUTINE WAS
; ENTERED
; REGISTERS EAX, ECX, AND EDX ARE CONSIDERED VOLATILE AND ASSUMED TO BE CHANGED
; WHILE THE REGISTERS BELOW MUST BE PRESERVED IF THE USER IS CHANGING THEM

add esp, 32
pop edi
pop esi
pop ebx
mov esp, ebp
pop ebp

;==============================================================================

ret
}
}

[shines77]

复数乘法的计算公式是这样的：

product.real = z1.real * z2.real - z1.imag * z2.imag
product.imag = z1.real * z2.imag + z1.imag * z2.real

例如：
z1 * z2 = (4 + 3i)(5 + 2i) = [4 * 5 - 3 * 2] + [3 * 5 + 4 * 2]i = 14 + 23i

[yaos]

感觉有点复杂

是不是简化一下，先别用一次4个复数乘法
有时候，多不是快阿

[shines77]

现在我在朋友的P4电脑试了一下我的SSE2程序,好象输入参数就自动关闭了:(

郁闷,无法调试

[yaos]

是啊，我执行了你的程序，出内存泄漏错误
仔细找吧 ：）

[shines77]

我没法调啊, 我回去看看

to yaos:
你输入0, 正常退出是不会有内存泄露的

[yaos]

嘿嘿
你自己的源代码，我都看着头大 ：）

测试帮着做，修改 ：） 还是我修改自己书上的C代码吧

[yaos]

根本不能运算的 :)

[shines77]

SSE的float版本是一次计算4个复数，一个float是4字节，一个float复数是8字节，一个XMM寄存器是128bit=16字节，一共有8个XMM寄存器，配合u,v流水线，一次计算4个复数是非常合理的，那是官方的版本，不容置疑，AMD的文章最后说，可能对于你的SSE清零的函数也有些帮助

The last optimization in both implementations is the use of the MOVNTQ and MOVNTPS instructions, nontemporal writes to memory that stream data to main memory. These instructions increase throughput to memory and make more efficient use of the bandwidth provided by the processor and memory controller. Nontemporal writes, such as MOVNTQ, MOVNTPS, and MOVNTDQ, should only be used on data that is not going to be accessed again in the near future.

SSE2的double版本是仿照float的写的，虽没有测试过，但应该是对的，有空麻烦帮我测试一下速度和正确性，源码里有详细说明，顺便也提供了一个最基本最原始的FFT乘法的代码

再补充一点，的确，不一定是多就是快，但是根据充分利用流水线和XMM寄存器128bit的特性，只能这么写，就像写16bit, 24bit, 32bit的AlphaBlending函数一样，都是并行处理几个单元的，MMX版的16bit AlphaBlending一次处理4个象素，这是由指令处理的容量和CPU的处理流程决定的，当然不排除有更好的方案，但是常规的思想和做法都是这么走的，最终的目的都是最少的CPU时钟周期

而且理顺流程以后其实并不是太复杂：） 而且，AMD的代码有非常详细的注释

[shines77]

忘了说，关键是看SSE2的double版本啦，它一次只能计算2个double复数，因为一个XMM寄存器的个数有限，而一个XMM寄存器刚好能放下1个double复数，同时处理2个，刚好做到指令配对，你可以看到SSE的float版本基本上指令也是配对的

我原来给你的测试代码，我有个地方写成了

shufpd xmm7, xmm7, 00000001b ; xmm6=[y0r,-y0i] // 2句都是xmm7了，这句应该是xmm6
shufpd xmm7, xmm7, 00000001b ; xmm7=[y1r,-y1i]

所以慢可能是指令不配对，还有计算结果不正常造成的，不过现在我估计也快不了多少

再拜托一句，有谁有空代为测试一下，下载（有完整源码）：
http://www.77studio.net/files/cmplx_src.rar

欢迎你来修改