1. Overview
SIMD stands for Single Instruction Multiple Data, Single Instruction Multiple Data Stream, a set of instructions that can read multiple operands and package them in large registers. After obtaining multiple operands at one time, they are stored in a large register and then operated, thereby achieving the effect of completing the calculation of multiple objects with one instruction and achieving acceleration. Currently, common compilers have relatively good support for 128-bit SIMD calculations on X86-64 CPUs. Basically, SIMD can be used to perform a simple optimization for most simple calculations. However, for more complex operations, it is still necessary to manually write the corresponding C/C++ or assembly code.
Guidelines for using SIMD instructions on Intel’s official website
https://software.intel.com/sites/landingpage/IntrinsicsGuide/#
2. Simple to use
I'm not very good at using SIMD yet, including testing some examples that may be slower with SIMD, and my CPU is older and only supports AVX2.0 (AVX-512 seems to be supported by current i9 and i7 9800x, but I, who is still using a fourth-generation i7, can't afford this thing). Therefore, for double-precision floating-point operations, only 256 bits of data can be packed at one time, that is, four double-precision floating-point numbers. By the way, one thing to note here is that when using gcc to compile code with immintrin.h header file functions, you need to add corresponding compilation options according to the function category you use (old SSE is supported by default and usually does not require additional options), such as AVX2.0, you need to use -mavx2, FMA needs to use -mfma, etc., otherwise the compilation will report an error. If you don't know what instruction set is used, but you are sure that your computer supports this operation, then just -march=native. In fact, if you use the AVX-512F instruction set, then -mavx512f can also be successfully compiled. However, if your CPU does not support AVX-512F, you will be prompted for an illegal instruction (#Funny) when you run the compiled file.
2.1 Vector addition
Yes, it’s vector addition again. This time we won’t do anything fancy, it’s just vector addition in C.
#include <iostream>
using namespace std;
#define N 20000000
int main(){ double *x,*y,*z; x=(double*)malloc(sizeof(double)*N); y=(double*)malloc(sizeof(double)*N); z=(double*)malloc(sizeof(double)*N); for(int i=0;i<N;i++) z[i]=x[i]+y[i]; free(x); free(y); free(z); return 1;
}
Then we use SIMD to write another one. To use SIMD technology, you need Intrinsics header files. There are many different header files. For specific ones, you can check my previous article on the introduction of the AVX instruction set. However, for general calculations, immintrin.h is basically sufficient.
#include <iostream>
#include <immintrin.h>
#define N 20000000
using namespace std;
int main(){ double *x,*y,*z,*px,*py,*pz; x=(double*)_mm_malloc(sizeof(double)*N,16);//Apply for memory and align the address according to the 4th power of 2 y=(double*)_mm_malloc(sizeof(double)*N,16); z=(double*)_mm_malloc(sizeof(double)*N,16); px=x;py=y;pz=z; __m128d vx,vy,vz;// __m128d is the data type corresponding to the operation of double-precision floating point numbers in the SSE instruction set for(int i=0;i<N/2;i++){ vx=_mm_load_pd(px);//Retrieve two numbers from the memory pointed to by px and put them into vx vy=_mm_load_pd(py);//Take out two numbers from the memory pointed to by py and put them into vy vz=vx+vy;//Calculate vx+vy and put the result into vz. This line also has a corresponding function. However, if compiled with GCC, there is no problem in writing it directly. The test found that neither the VS compiler nor the Intel compiler supports this writing method. _mm_store_pd(pz,vz);//Put the results in vz into the memory pointed to by pz px+=2;//Since two pieces of data were taken out previously, the pointer moves back two places. py+=2; pz+=2; } _mm_free(x);//Release the memory requested by _mm_malloc _mm_free(y); _mm_free(z); return 1;
}
After the program is written, the next step is to test it. After testing, the calculation results of the two programs are correct. Since the calculation results need to be written and assigned, I don’t think it is troublesome. If you are interested, you can assign a value yourself and see if the result is correct. Let’s test the efficiency here. The first is still the test under the O0 condition of g++.
Time spent without SIMD
real 0m0.167s
user 0m0.092s
sys 0m0.075s
Time to use SIMD
real 0m0.175s
user 0m0.122s
sys 0m0.053s
It can be seen that SIMD is slower under O0 optimization, but the time difference between the two is not big. After multiple tests, it seems that the range basically fluctuates up and down. It can be considered that the two efficiencies are equivalent, or SIMD is slightly slower.
Next is the test under O3 conditions
When SIMD is not used
real 0m0.112s
user 0m0.040s
sys 0m0.072s
When using SIMD
real 0m0.113s
user 0m0.044s
sys 0m0.068s
There is basically no difference between the two. I have tried increasing the array dimension, but since the memory is only 4G, the increase is not much. I have also tried repeated addition. The final test result is that regardless of whether compiler optimization is used or not, SIMD has not improved efficiency. In other words, SIMD has no effect in this example and even makes the calculation slower.
Later, I read some posts online and some people encountered a similar situation to me.
Someone mentioned that when performing some simple operations, the compiler will automatically convert to SIMD code through some optimization techniques without you actively using SIMD. Therefore, the code you write using SIMD will not exceed the results of the compiler's automatic optimization at the fastest.
Some people say that the additional overhead of SIMD is relatively large. For this simple addition, the computational efficiency improvement brought by using SIMD cannot cover up the additional performance loss caused by using SIMD. SIMD may only have an acceleration effect when the calculation is more complex.
There is another saying that there are a total of 8 XMM registers, but all 8 registers are not used here, so the full performance is not used, so it has no effect.
By generating assembly code with the -S option of g++, you can see that the compiler automatically optimizes this part of the code. Whether you use the function in the immintrin.h header file or use a loop to calculate directly, the generated assembly code is basically the same. In other words, the compiler has automatically optimized this part of the operation, so it has no effect.
2.2 Complex array multiplication
The multiplication of complex arrays was tested when I was writing a fast Fourier transform implementation. I found that the efficiency was indeed significantly improved after using AVX2.0. Here is a brief demonstration.
The first is the principle of calculation, because complex number calculations are different from real numbers. According to the relationship between the real part and the imaginary part, different calculation methods are often required due to different positions of floating point numbers during calculation. Therefore, some processing and changes are also required when using SIMD acceleration. The AVX2.0 instruction is used here to operate 256 bits of data at a time, that is, 4 double-precision floating-point numbers and 2 double-precision complex numbers. Therefore, two complex arrays z1[2] and z2[2] are used for multiplication in the test here. The specific calculation is as follows:
To make it easier for us to remember:
The two complex numbers in z1 are a0 + i * b0 and a1 + i * b1
The two complex numbers in z2 are c0+i*d0 and c1+i*d1
At this time, the calculation result is expressed as:
(a0 * c0 - b0 * d0) + i * (a0 * d0 + b0 * c0) and (a1 * c1 - b1 * d1) + i * (a1 * d1+b1 * c1)
The calculation results can be broken down into two parts:
{a0 * c0 , a0 * d0 , a1 * c1 , a1 * d1}
{b0 * d0 , b0 * c0 , b1 * d1 , b1 * c1}
Subtract the even subscripts of these two parts and add the odd subscripts to get the final calculation result. It happens that AVX2.0 provides us with this instruction.
Break it down further:
{a0 , a0 , a1 , a1}
{c0 , d0 , c1 , d1}
Multiplied together we get the first part and
{b0 , b0 , b1 , b1}
{d0 , c0 , d1 , c1}
Multiplying together yields the second part
It can be seen that the complex multiplication calculation can be performed by taking out the real part in z1 and multiplying it directly with z2, then taking out the imaginary part in z1 and exchanging the real part and imaginary part of z2 for multiplication, and finally interleaving the addition and subtraction to obtain the final result.
Although the operation seems very cumbersome, fortunately, these operations can actually be completed by finding the corresponding instructions.
In the immintrin.h header file, there is a __m256d _mm256_permute_pd (__m256d a, int imm8) function.
This function receives a __m256d variable and a 4-digit immediate number (decimal literal value of 0-15). Based on the number on each binary digit of imm8, it determines whether the number at the corresponding position in a should remain unchanged or replaced with the numbers before and after it. Finally, a new __m256d variable is obtained and returned. The specific rules are as follows:
IF (imm8[0] == 0) dst[63:0] := a[63:0]
IF (imm8[0] == 1) dst[63:0] := a[127:64]
IF (imm8[1] == 0) dst[127:64] := a[63:0]
IF (imm8[1] == 1) dst[127:64] := a[127:64]
IF (imm8[2] == 0) dst[191:128] := a[191:128]
IF (imm8[2] == 1) dst[191:128] := a[255:192]
IF (imm8[3] == 0) dst[255:192] := a[191:128]
IF (imm8[3] == 1) dst[255:192] := a[255:192]
The basic principle is that two numbers are a group. 00, 01, 10, and 11 respectively represent that the two numbers are exchanged, remain unchanged, or the former one covers the next one, and the latter one covers the previous one.
As mentioned earlier, we need a __m256d variable that only contains the real part of z1 and a __m256d variable that only contains the imaginary part of z1.
Here, setting imm8=15 can complete the replacement of the real part by the imaginary part, and the real part covering the imaginary part can be achieved by setting imm8=0. However, the same effect can also be achieved by using the _mm256_movedup_pd function (here, in SSE3, 4 single steps can be completed directly through _mm_moveldup_ps The operation of covering even subscripts by odd subscripts in precision numbers, but this instruction is only available for operating 4 single-precision floating-point numbers. In other cases, only even-numbered subscripts are supported to cover odd-numbered subscripts, and the reverse is not supported. Is there an advantage in hardware implementation for 4 single-precision floating-point number operations? ) In this way, we can get the corresponding _mm256d variable that only contains the real part and only the imaginary part.
Next, let imm8=10 to reverse the real and imaginary parts of the two complex numbers in z2. In this way, the four _mm256d variables required for calculation are successfully obtained, and then the calculation can be performed.
FMA multiply-accumulate fusion can complete the operation of a * b + c in one go. This instruction has faster speed and higher accuracy than separate calculations. At the same time, _mm256_fmaddsub_pd is also provided in immintrin.h to complete the operation of adding even subscript bits and subtracting odd subscript bits.
The specific implementation is as follows:
void mul(complex<double> *x,complex<double> *y){ __m256d tx,ty,r; tx=_mm256_loadu_pd((double*)(x)); ty=_mm256_loadu_pd((double*)(y)); r=_mm256_permute_pd(ty,5)*_mm256_permute_pd(tx,15); _mm256_storeu_pd((double*)(x),_mm256_fmaddsub_pd(ty,_mm256_movedup_pd(tx),r));
}
Next, test the efficiency
#define N 10000
int main(){ complex<double> x[2],y[2]; x[0]=complex<double>(1,0); x[1]=complex<double>(-1,0); y[0]=complex<double>(0,1); y[1]=complex<double>(0,-1); for(int j=0;j<N;j++){ for(int i=0;i<N;i++){ mul(x,y); } } cout<<x[0]<<'\t'<<x[1]<<endl; return 0;
}
Using g++ -O3 optimization, the time of multiplication directly using the * overloaded in the standard library is:
(1,-0) (-1,-0) real 0m0.714s
user 0m0.714s
sys 0m0.000s
The time of the mul function after using SIMD acceleration is:
(1,-0) (-1,-0) real 0m0.489s
user 0m0.485s
sys 0m0.004s
It can be seen that the efficiency of complex number calculations has been significantly improved after using SIMD acceleration. Therefore, for the example of array addition, due to the simple calculation, the compiler has already done the corresponding SIMD acceleration. Writing the corresponding operations yourself cannot speed up the speed.
2.3 Matrix transpose
I recently tested using SIMD to implement the matrix transposition method and found that the acceleration effect is good.
Let’s talk about the algorithm first. The AVX instruction set can operate 4 double-precision floating point numbers at one time. Therefore, by reading four __m256d data, a 4x4 square matrix is transposed at a time. For square matrices larger than 4x4, you can first complete the transposition of each small block, and then transpose the position of the sub-block. For non-original transposition, just write the transposed sub-block to the corresponding position.
Taking the 4x4 matrix as an example, assume that the original matrix is
a0,a1,a2,a3
b0,b1,b2,b3
c0,c1,c2,c3
d0,d1,d2,d3
After the transposition is completed, it should be
a0,b0,c0,d0
a1,b1,c1,d1
a2,b2,c2,d2
a3,b3,c3,d3
Here we need to use the four operations of _mm256_unpacklo_pd, _mm256_unpackhi_pd, _mm256_permute4x64_pd, and _mm256_blend_pd in the AVX and AVX2 instruction sets.
_mm256_unpacklo_pd accepts two __m256d parameters, and returns the first and third numbers of the two incoming parameters in order to form a new __m256d. For example, when a0, a1, a2, a3 and b0, b1, b2, b3 are passed in, a0, b0, a2, b2 are returned.
_mm256_unpackhi_pd is basically similar to _mm256_unpacklo_pd, except that it obtains the second and fourth numbers of the incoming value and returns them
_mm256_permute4x64_pd accepts a __m256d parameter and an 8-bit binary number. Each pair of 8-bit binary numbers is assigned to the number of the passed __m256d to the number of the returned __m256d according to its position in the corresponding 8-bit binary number. For example, pass in a0, a1 , a2, a3 and 0b11100100 (converted to four decimal numbers are 3,2,1,0), then the returned value is a3, a2, a1, a0. One thing to note here is the data arrangement problem in the register. If we want to invert the entire sequence in the memory, the binary number that should actually be given is 0b00011011
_mm256_blend_pd accepts two __m256d parameters a and b, and a four-digit binary number. Depending on whether each bit in the binary number is 1 or 0, it is determined that the corresponding position of the new __m256d generated is a. The value is still the value of b. For example, if a0, a1, a2, a3 and b0, b1, b2, b3 and 0b0101 are passed in, the generated data should be a0, b1, a2, b3. Similarly, you need to pay attention to the arrangement of the data in the register.
With these four functions, we can design the transposition algorithm. First, we work backwards
To get a0,b0,c0,d0 and a1,b1,c1,d1
We consider using _mm256_unpacklo_pd and _mm256_unpackhi_pd to generate two transposed sequences simultaneously
Then we should have a0, a1, c0, c1 and b0, b1, d0, d1
We found a0, a1, c0, c1 and can generate a0, a1, a2, a3 and c2, c3, c0, c1 through _mm256_blend_pd
A0, a1, a2, a3 are the original sequences, c2, c3, c0, c1 can be obtained by rearranging c0, c1, c2, c3 through _mm256_permute4x64_pd, and b0, b1, d0, d1 can also be obtained by the same method. By rearranging and shuffling in this way, we can achieve a 4x4 matrix transpose. The specific code is as follows:
double* x=(double*)_mm_malloc(sizeof(double)*4*4,32); double* y=(double*)_mm_malloc(sizeof(double)*4*4,32); double *p1=x,*p2=x+4,*p3=x+8,*p4=x+12; double *d1=x,*d2=x+4,*d3=x+8,*d4=x+12, int k=0; for(int i=0;i<4;i++) { for(int j=0;j<4;j++) { x[k++]=i; } } __m256d s1,s2,s3,s4,t1,t2,t3,t4,t5,t6,t7,t8; s1=_mm256_load_pd(p1); s2=_mm256_load_pd(p2); s3=_mm256_load_pd(p3); s4=_mm256_load_pd(p4); t1=_mm256_permute4x64_pd(s1,0b01001110);//First rearrange the four sequences, swap the first two and the last two with a2, a3, a0, a1 t2=_mm256_permute4x64_pd(s2,0b01001110);//b2,b3,b0,b1 t3=_mm256_permute4x64_pd(s3,0b01001110);//c2,c3,c0,c1 t4=_mm256_permute4x64_pd(s4,0b01001110);//d2,d3,d0,d1 t5=_mm256_blend_pd(s1,t3,0b1100);//Merge a swapped sequence with an original sequence to complete the final step of shuffling a0, a1, c0, c1 t6=_mm256_blend_pd(s2,t4,0b1100);//b0,b1,d0,d1 t7=_mm256_blend_pd(t1,s3,0b1100);//a2,a3,c2,c3 t8=_mm256_blend_pd(t2,s4,0b1100);//b2,b3,d2,d3 s1=_mm256_unpacklo_pd(t5,t6);//Shuffle the merged sequence to obtain the transposed sequence a0, b0, c0, d0 s2=_mm256_unpackhi_pd(t5,t6);//a1,b1,c1,d1 s3=_mm256_unpacklo_pd(t7,t8); s4=_mm256_unpackhi_pd(t7,t8); _mm256_store_pd(d1,s1); _mm256_store_pd(d2,s2); _mm256_store_pd(d3,s3); _mm256_store_pd(d4,s4); k=0 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { cout<<y[k++]<<'\t'; } cout<<'\n'; }
In this way, the 4x4 matrix transposition is completed. The next step is to conduct a large-scale matrix transposition to test the performance. Taking the 4096 square matrix as an example, the SIMD transposition code is as follows:
#include <immintrin.h>
#include <iostream>
using namespace std;
#define N 4096
int main()
{ double* x=(double*)_mm_malloc(sizeof(double)*N*N,32); double* y=(double*)_mm_malloc(sizeof(double)*N*N,32); int k=0; for(int i=0;i<16;i++) { for(int j=0;j<16;j++) { x[k++]=i; } } double *p1,*p2,*p3,*p4; double *d1,*d2,*d3,*d4,*t=y; d1=y;d2=y+N;d3=y+2*N;d4=y+3*N; p1=x;p2=x+N;p3=x+2*N;p4=x+3*N; t+=4; for(int i=0;i<N/4;i++) { for(int j=0;j<N/4;j++) { __m256d s1,s2,s3,s4,t1,t2,t3,t4,t5,t6,t7,t8; s1=_mm256_load_pd(p1); s2=_mm256_load_pd(p2); s3=_mm256_load_pd(p3); s4=_mm256_load_pd(p4); t1=_mm256_permute4x64_pd(s1,0b01001110); t2=_mm256_permute4x64_pd(s2,0b01001110); t3=_mm256_permute4x64_pd(s3,0b01001110); t4=_mm256_permute4x64_pd(s4,0b01001110); t5=_mm256_blend_pd(s1,t3,0b1100); t6=_mm256_blend_pd(s2,t4,0b1100); t7=_mm256_blend_pd(t1,s3,0b1100); t8=_mm256_blend_pd(t2,s4,0b1100); s1=_mm256_unpacklo_pd(t5,t6); s2=_mm256_unpackhi_pd(t5,t6); s3=_mm256_unpacklo_pd(t7,t8); s4=_mm256_unpackhi_pd(t7,t8); _mm256_store_pd(d1,s1); _mm256_store_pd(d2,s2); _mm256_store_pd(d3,s3); _mm256_store_pd(d4,s4); p1+=4;p2+=4;p3+=4;p4+=4; d1+=4*N;d2+=4*N;d3+=4*N;d4+=4*N; } p1+=3*N;p2+=3*N;p3+=3*N;p4+=3*N; d1=t;d2=t+N;d3=t+2*N;d4=t+3*N; t+=4; } cout<<y[100]<<'\n'; /* k=0; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { cout<<y[k++]<<'\t'; } cout<<'\n'; } */ return 0;
}
And the original code is as follows:
#include <immintrin.h>
#include <iostream>
using namespace std;
#define N 4096
int main()
{ double* x=(double*)_mm_malloc(sizeof(double)*N*N,32); double* y=(double*)_mm_malloc(sizeof(double)*N*N,32); int k=0; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { x[k++]=i; } } double *s=x,*d=y,*t=y; t++; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { *d=*(x++); d+=N; } d=t; t++; }
cout<<y[100]<<'\n';
/* k=0; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { cout<<y[k++]<<'\t'; } cout<<'\n'; }
*/ return 0;
}
The tests were compiled using -O2 optimization, and the final test results are as follows:
SIMD not used
real 0m0.453s
user 0m0.393s
sys 0m0.061s
Use SIMD
real 0m0.147s
user 0m0.087s
sys 0m0.060s
It can be seen that through SIMD acceleration, the performance of the square matrix transposition of 4096 has been improved by nearly 3 times.
2.4 Linear congruence generation of pseudo-random numbers
Recently I needed to use SIMD to make a random number generator. I thought it would be easy, but it turned out that AVX's integer operations were full of pitfalls. After a day of hard work, I was able to create a double-precision floating-point number generator that generates a uniform distribution of -1 to 1.
Not much to say about linear congruence. Use seed=(seed*a+c) mod m seed=(seed*a+c)\ mod\ mseed=(seed*a+c) mod m to continuously update the seed to obtain uniformly distributed pseudo-random numbers from 0 to m-1. Here we take m=232−1, a=513, c=9973 m=2^{32}-1, a=513, c=9973m=2
32
−1,a=513,c=9973
Because I have been dealing with floating-point numbers before, there is nothing particularly strange about the calculations. I thought that I would just write the integers accordingly. It turned out that after writing it all, the result was very different. First of all, AVX-512 has only begun to support floating-point-like load and store instructions (functions) for integer types, but I do not have a machine that supports AVX-512. To read and write integers with AVX, you need to use mask_load and mask_store. In addition to the addresses and ymm register variables required by normal laod and store, these two functions also require an additional mask parameter to control which locations of integer data are read and written, and which locations are directly set to zero if they are not read or written. This mask is also a ymm register variable. Whether the data is read or written is determined according to whether the highest bit of each position number in the mask parameter is 1 or 0. This means that if the corresponding position is given a negative number, the corresponding position will be read and written, and if 0 or an integer is given, the corresponding position will be set to zero directly. So if you want to read the four seeds 1, 2, 3, and 4, you need to do the following operations.
unsigned t[8]={1,1,2,2,3,3,4,4};
__m256i k=_mm256_maskload_epi32((int*)t,_mm256_set1_epi32(-1));
After reading the seed, the next step is to use linear congruence to generate random numbers. There are only two integer multiplication functions of AVX, _mm256_mul_epu32 and _mm256_mul_epi32. i means signed and u means unsigned. There is no difference in the others. This multiplication does not replace the 8 ints or unsigned in __m256i Int corresponds to multiplication, but treats the data as four unsigned integers, then truncates the high bits, leaving only the lower 32 bits for multiplication, and then returns __m256i, which contains four 64-bit signed or unsigned integers. However, the problem here is not big, that is, 4 less random numbers are generated. The calculated results can be calculated without any modification and continue to calculate new random numbers.
void Rand(__m256i &seed)
{ seed=_mm256_mul_epu32(seed,_mm256_set1_epi32(513)); seed=_mm256_add_epi64(seed,_mm256_set1_epi64x(9973)); seed=_mm256_and_si256(seed,_mm256_set1_epi64x(4294967295));
}
In fact, if you simply need unsigned integer random numbers, this function is enough. But since what I need is a floating point number uniformly distributed from -1 to 1, this result needs to be further processed. However, the function that converts four unsigned ___m256i to __m256d requires AVX512 support, and the function that converts four 64-bit unsigned __m256i into four unsigned __m128i cannot be found at all. However, there is a function in the cast function that forces conversion from __m256i to __m128i. This function directly forces truncation of half of __m256i to generate a new __m128i. However, if the seed is regarded as eight 32-bit integers, the random numbers we generate are at the four positions of 0, 2, 4, and 6. If they are directly truncated, only two numbers will be left. This requires rearranging the data, and then using the permute function. The _mm256_permutevar8x32_epi32 function is used here for rearrangement. This function is very useful. The number to determine which original position is placed in each position is determined by the __m256i variable. Directly give a __mm256_set_ep i32(p1,p2,…,p8) will do (p1,p2,…,p8 are the data numbers that need to be placed at the corresponding positions. For example, if you want to invert, it is 7,6,5,4,3,2,1,0. Pay attention to the problem of big endian and small endian. The specific number needs to be tested, here is just a hint). In this way, the four random numbers can be concentrated into the first 128 bits, and then directly truncated to obtain the __m128i variable containing four 32-bit integers. Then this variable can be converted into a __m256d variable containing 4 doubles through cvt, and another division can generate 4 double-precision floating-point numbers uniformly distributed from -1 to 1 (double precision here is just a data type, and the actual random number accuracy is only 1/232 {1}/{2^{32}}1/2
32
). What needs to be noted here is that this conversion of double-precision floating point numbers is read as signed, so the integers removed when converted to -1 to 1 are half of the m value in the Rand function.
Rand(k);
__m256i rd=_mm256_permutevar8x32_epi32(k,_mm256_set_epi32(0,0,0,0,0,2,4,6));
__m128i rrd=_mm256_castsi256_si128 (rd);
__m256d u=_mm256_cvtepi32_pd(rrd);
_mm256_storeu_pd(x,u/_mm256_set1_pd(4294967295));
The complete code is as follows
#include <iostream>
#include <immintrin.h>
#include <cstdlib>
#include <ctime>
using namespace std;
void Rand(__m256i &seed)
{ seed=_mm256_mul_epu32(seed,_mm256_set1_epi32(513)); seed=_mm256_add_epi64(seed,_mm256_set1_epi64x(9973)); seed=_mm256_and_si256(seed,_mm256_set1_epi64x(4294967295));
}
int main()
{ unsigned t[8]={1,1,2,2,3,3,4,4},tj[16]={}; double x[4]; __m256i k=_mm256_maskload_epi32((int*)t,_mm256_set1_epi32(1)); auto s=clock(); for(int i=0;i<10000;i++) { Rand(k); __m256i rd=_mm256_permutevar8x32_epi32(k,_mm256_set_epi32(0,0,0,0,0,2,4,6)); __m128i rrd=_mm256_castsi256_si128 (rd); __m256d u=_mm256_cvtepi32_pd(rrd); _mm256_storeu_pd(x,u/_mm256_set1_pd(2147483647)); for(int j=0;j<4;j+=1) { cout<<x[j]<<'\t'; //tj[int(x[j]*8+8)]++; } cout<<'\n'; } auto e=clock(); for(int i=0;i<40000;i++) { for(int j=0;j<4;j++)cout<<double(rand())/RAND_MAX<<'\t'; cout<<'\n'; } auto e1=clock(); cout<<"time = "<<double(e-s)/CLOCKS_PER_SEC<<'\t'<<double(e1-e)/CLOCKS_PER_SEC; for(int i=0;i<16;i++) cout<<tj[i]<<'\n'; return 0;
}
Here we use tj to make a statistics to test whether the generated random numbers are uniformly distributed, and then also compare the time difference with using the rand function of the standard library to generate the same number of random numbers.
The results of separate statistics are as follows:
2588
2532
2536
2448
2592
2416
2572
2788
2688
2420
2420
2396
2356
2356
2668
It can be seen that -1 to 1 are evenly divided into 16 parts and 40,000 samples are counted. The number of random numbers falling in each interval is basically around 2,500, which proves that there is no big problem in this calculation process.
Then take a look at the efficiency. Since it does not print to the console, the calculation process will be directly optimized after optimization is turned on. However, printing to the console is extremely time-consuming. Here is a simple test. Redirect the output results directly to the file, so that the output time will be shorter and the impact on the results will be smaller. Using gcc's O3 optimization, the time statistics are as follows:
0.020918s (using AVX) 0.0904440s (standard library rand function)
You can see that the speed of generating random numbers is significantly faster after using AVX.