1 // Copyright ©2016 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
5 //+build !noasm,!appengine
10 #define MOVDDUP_X2_X3 BYTE $0xF2; BYTE $0x0F; BYTE $0x12; BYTE $0xDA
12 #define MOVDDUP_X4_X5 BYTE $0xF2; BYTE $0x0F; BYTE $0x12; BYTE $0xEC
14 #define MOVDDUP_X6_X7 BYTE $0xF2; BYTE $0x0F; BYTE $0x12; BYTE $0xFE
16 #define MOVDDUP_X8_X9 BYTE $0xF2; BYTE $0x45; BYTE $0x0F; BYTE $0x12; BYTE $0xC8
19 #define ADDSUBPD_X2_X3 BYTE $0x66; BYTE $0x0F; BYTE $0xD0; BYTE $0xDA
21 #define ADDSUBPD_X4_X5 BYTE $0x66; BYTE $0x0F; BYTE $0xD0; BYTE $0xEC
23 #define ADDSUBPD_X6_X7 BYTE $0x66; BYTE $0x0F; BYTE $0xD0; BYTE $0xFE
25 #define ADDSUBPD_X8_X9 BYTE $0x66; BYTE $0x45; BYTE $0x0F; BYTE $0xD0; BYTE $0xC8
27 // func AxpyUnitary(alpha complex128, x, y []complex128)
28 TEXT ·AxpyUnitary(SB), NOSPLIT, $0
29 MOVQ x_base+16(FP), SI // SI = &x
30 MOVQ y_base+40(FP), DI // DI = &y
31 MOVQ x_len+24(FP), CX // CX = min( len(x), len(y) )
33 CMOVQLE y_len+48(FP), CX
34 CMPQ CX, $0 // if CX == 0 { return }
36 PXOR X0, X0 // Clear work registers and cache-align loop
38 MOVUPS alpha+0(FP), X0 // X0 = { imag(a), real(a) }
40 SHUFPD $0x1, X1, X1 // X1 = { real(a), imag(a) }
42 MOVAPS X0, X10 // Copy X0 and X1 for pipelining
45 ANDQ $3, CX // CX = n % 4
46 SHRQ $2, BX // BX = floor( n / 4 )
47 JZ caxy_tail // if BX == 0 { goto caxy_tail }
50 MOVUPS (SI)(AX*8), X2 // X_i = { imag(x[i]), real(x[i]) }
51 MOVUPS 16(SI)(AX*8), X4
52 MOVUPS 32(SI)(AX*8), X6
53 MOVUPS 48(SI)(AX*8), X8
55 // X_(i+1) = { real(x[i], real(x[i]) }
61 // X_i = { imag(x[i]), imag(x[i]) }
67 // X_i = { real(a) * imag(x[i]), imag(a) * imag(x[i]) }
68 // X_(i+1) = { imag(a) * real(x[i]), real(a) * real(x[i]) }
79 // imag(result[i]): imag(a)*real(x[i]) + real(a)*imag(x[i]),
80 // real(result[i]): real(a)*real(x[i]) - imag(a)*imag(x[i])
87 // X_(i+1) = { imag(result[i]) + imag(y[i]), real(result[i]) + real(y[i]) }
89 ADDPD 16(DI)(AX*8), X5
90 ADDPD 32(DI)(AX*8), X7
91 ADDPD 48(DI)(AX*8), X9
92 MOVUPS X3, (DI)(AX*8) // y[i] = X_(i+1)
93 MOVUPS X5, 16(DI)(AX*8)
94 MOVUPS X7, 32(DI)(AX*8)
95 MOVUPS X9, 48(DI)(AX*8)
98 JNZ caxy_loop // } while --BX > 0
99 CMPQ CX, $0 // if CX == 0 { return }
103 MOVUPS (SI)(AX*8), X2 // X_i = { imag(x[i]), real(x[i]) }
104 MOVDDUP_X2_X3 // X_(i+1) = { real(x[i], real(x[i]) }
105 SHUFPD $0x3, X2, X2 // X_i = { imag(x[i]), imag(x[i]) }
106 MULPD X1, X2 // X_i = { real(a) * imag(x[i]), imag(a) * imag(x[i]) }
107 MULPD X0, X3 // X_(i+1) = { imag(a) * real(x[i]), real(a) * real(x[i]) }
110 // imag(result[i]): imag(a)*real(x[i]) + real(a)*imag(x[i]),
111 // real(result[i]): real(a)*real(x[i]) - imag(a)*imag(x[i])
115 // X_(i+1) = { imag(result[i]) + imag(y[i]), real(result[i]) + real(y[i]) }
117 MOVUPS X3, (DI)(AX*8) // y[i] = X_(i+1)
118 ADDQ $2, AX // i += 2
119 LOOP caxy_tail // } while --CX > 0