I는 희소 배열 a
(거의 제로)가 :희소 배열 압축
unsigned char a[1000000];
및 I는 SIMD 명령어를 사용 a
0이 아닌 요소 인덱스의 배열 b
를 작성하고자 AVX2가있는 인텔 x64 아키텍처 나는 그것을 효율적으로하는 방법에 대한 조언을 찾고있다. 특히, SIMD 레지스터에 연속적으로 배열 된 0이 아닌 요소의 위치를 얻기 위해 SIMD 명령어가 있습니까?
I는 희소 배열 a
(거의 제로)가 :희소 배열 압축
unsigned char a[1000000];
및 I는 SIMD 명령어를 사용 a
0이 아닌 요소 인덱스의 배열 b
를 작성하고자 AVX2가있는 인텔 x64 아키텍처 나는 그것을 효율적으로하는 방법에 대한 조언을 찾고있다. 특히, SIMD 레지스터에 연속적으로 배열 된 0이 아닌 요소의 위치를 얻기 위해 SIMD 명령어가 있습니까?
비록 AVX2 명령어 세트에는 많은 GATHER 명령어가 있지만 성능은 너무 느립니다. 가장 효과적인 방법은 배열을 수동으로 처리하는 것입니다. 당신이 제로가 아닌 요소의 수 (즉, 훨씬 적은 % 1에 비해) 매우 낮은 것으로 예상한다면, 당신은 단순히 제로가 아닌 것에 대해 각 16 바이트 덩어리를 확인할 수 있습니다
:
int mask = _mm_movemask_epi8(_mm_cmpeq_epi8(reg, _mm_setzero_si128());
if (mask != 65535) {
//store zero bits of mask with scalar code
}
하면 좋은 요소의 비율 잘못 예측 된 분기 비용과 'if'내부의 느린 스칼라 코드 비용은 무시할 수 있습니다.
좋은 일반적인 솔루션에 관해서는, 첫째 스트림 압축의 SSE의 구현을 고려한다. 보시다시피
__m128i shuf [65536]; //must be precomputed
char cnt [65536]; //must be precomputed
int compress(const char *src, int len, char *dst) {
char *ptr = dst;
for (int i = 0; i < len; i += 16) {
__m128i reg = _mm_load_si128((__m128i*)&src[i]);
__m128i zeroMask = _mm_cmpeq_epi8(reg, _mm_setzero_si128());
int mask = _mm_movemask_epi8(zeroMask);
__m128i compressed = _mm_shuffle_epi8(reg, shuf[mask]);
_mm_storeu_si128((__m128i*)ptr, compressed);
ptr += cnt[mask]; //alternative: ptr += 16-_mm_popcnt_u32(mask);
}
return ptr - dst;
}
(_mm_shuffle_epi8
+ 룩업 테이블) 경이를 할 수 있습니다 : 그것은 바이트 배열에서 모두 0 요소 (here에서 가져온 아이디어)을 제거합니다. 스트림 압축과 같은 구조적으로 복잡한 코드를 벡터화하는 다른 방법을 모르겠습니다.
이제는 색인을 얻으려는 경우에만 문제가 남습니다. 각 색인은 4 바이트 값으로 저장해야하므로 16 바이트의 입력 청크가 최대 64 바이트의 출력을 생성 할 수 있으며 이는 단일 SSE 레지스터에 맞지 않습니다.
이 문제를 해결하는 한 가지 방법은 정직하게 출력을 64 바이트로 압축 해제하는 것입니다. 따라서 reg
을 코드에서 상수 (0,1,2,3,4, ..., 15)로 바꾼 다음 SSE 레지스터를 4 개의 레지스터로 압축을 풀고 4 개의 i
값이있는 레지스터를 추가합니다. 이것은 더 많은 지시 사항을 취할 것입니다 : 6 개의 포장 풀기 지침, 4 개의 추가 및 3 개의 매장 (이미 하나 있음). 제게는 엄청난 오버 헤드입니다. 특히 0이 아닌 요소의 수가 25 % 미만이라고 예상 할 경우 더욱 그렇습니다.
또는 단일 루프 반복에 의해 처리 된 0이 아닌 바이트의 수를 4로 제한하여 하나의 레지스터가 출력을 위해 항상 충분하도록 할 수 있습니다. 이 방법의
__m128i shufMask [65536]; //must be precomputed
char srcMove [65536]; //must be precomputed
char dstMove [65536]; //must be precomputed
int compress_ids(const char *src, int len, int *dst) {
const char *ptrSrc = src;
int *ptrDst = dst;
__m128i offsets = _mm_setr_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
__m128i base = _mm_setzero_si128();
while (ptrSrc < src + len) {
__m128i reg = _mm_loadu_si128((__m128i*)ptrSrc);
__m128i zeroMask = _mm_cmpeq_epi8(reg, _mm_setzero_si128());
int mask = _mm_movemask_epi8(zeroMask);
__m128i ids8 = _mm_shuffle_epi8(offsets, shufMask[mask]);
__m128i ids32 = _mm_unpacklo_epi16(_mm_unpacklo_epi8(ids8, _mm_setzero_si128()), _mm_setzero_si128());
ids32 = _mm_add_epi32(ids32, base);
_mm_storeu_si128((__m128i*)ptrDst, ids32);
ptrDst += dstMove[mask]; //alternative: ptrDst += min(16-_mm_popcnt_u32(mask), 4);
ptrSrc += srcMove[mask]; //no alternative without LUT
base = _mm_add_epi32(base, _mm_set1_epi32(dstMove[mask]));
}
return ptrDst - dst;
}
한 가지 단점은 라인 ptrDst += dstMove[mask];
이 이전 반복에 실행될 때까지 지금 이후의 각 루프 반복을 시작할 수 있다는 것입니다 : 다음은 샘플 코드입니다. 따라서 중요한 경로가 크게 증가했습니다. 하드웨어 하이퍼 스레딩 또는 수동 에뮬레이션으로이 페널티를 없앨 수 있습니다.
본 것처럼,이 기본 아이디어에는 다양한 변형이 있습니다.이 모든 아이디어는 각기 다른 효율성으로 문제를 해결합니다. LUT의 크기를 줄이려면 LUT의 크기를 줄일 수 있습니다 (다시 말하면 처리량 성능 저하가 있음).
이 접근법은보다 넓은 레지스터 (즉,AVX2 및 AVX-512), 여러 개의 연속 된 명령을 단일 AVX2 또는 AVX-512 명령에 결합하여 처리량을 약간 높일 수 있습니다.
참고 : 사전 계산 LUT에 눈에 띄는 노력이 필요하기 때문에 코드를 테스트하지 않았습니다. nonzeros의 인덱스를 계산하기
LUT 접근 방식이 내 [답변] (http://stackoverflow.com/a)과 어떻게 비교되는지 보는 것이 좋을 것입니다./41958528/2439725)의 비트 조작 명령 (BMI1 및 BMI2)을 기반으로합니다. – wim
다섯 방법 :
세미 벡터화 루프 : 제로와 비교 문자를 가진 SIMD 벡터를로드하고 movemask 적용. 문자 중 하나라도 0이 아니면 (작은 문자가 @stgatilov이라고도 함) 작은 스칼라 루프를 사용하십시오. 이것은 매우 드문 드문 배열에 적합합니다. 아래 코드의 arr2ind_movmsk
함수는 스칼라 루프에 BMI1 명령어 을 사용합니다.
벡터화 된 루프 : Intel Haswell 프로세서 이상은 BMI1 및 BMI2 명령어 세트를 지원합니다. BMI2는 pext
인스트럭션 (Parallel bits extract, see wikipedia link), 이 여기에 유용합니다. 아래 코드에서 arr2ind_pext
을 참조하십시오.
if 문을 사용한 고전적인 스칼라 루프 : arr2ind_if
.
가지가없는 스칼라 루프 : arr2ind_cmov
.
찾아보기 테이블 : @stgatilov은 pdep 및 다른 정수 명령어 대신에 찾아보기 테이블을 사용할 수 있음을 보여줍니다. 이것은 잘 작동하지만 룩업 테이블은 상당히 크기 때문에 L1 캐시에 맞지 않습니다. 여기에서는 테스트하지 않습니다. 토론 here도 참조하십시오.
/*
gcc -O3 -Wall -m64 -mavx2 -fopenmp -march=broadwell -std=c99 -falign-loops=16 sprs_char2ind.c
example: Test different methods with an array a of size 20000 and approximate 25/1024*100%=2.4% nonzeros:
./a.out 20000 25
*/
#include <stdio.h>
#include <immintrin.h>
#include <stdint.h>
#include <omp.h>
#include <string.h>
__attribute__ ((noinline)) int arr2ind_movmsk(const unsigned char * restrict a, int n, int * restrict ind, int * m){
int i, m0, k;
__m256i msk;
m0=0;
for (i=0;i<n;i=i+32){ /* Load 32 bytes and compare with zero: */
msk=_mm256_cmpeq_epi8(_mm256_load_si256((__m256i *)&a[i]),_mm256_setzero_si256());
k=_mm256_movemask_epi8(msk);
k=~k; /* Search for nonzero bits instead of zero bits. */
while (k){
ind[m0]=i+_tzcnt_u32(k); /* Count the number of trailing zero bits in k. */
m0++;
k=_blsr_u32(k); /* Clear the lowest set bit in k. */
}
}
*m=m0;
return 0;
}
__attribute__ ((noinline)) int arr2ind_pext(const unsigned char * restrict a, int n, int * restrict ind, int * m){
int i, m0;
uint64_t cntr_const = 0xFEDCBA;
__m256i shft = _mm256_set_epi64x(0x04,0x00,0x04,0x00);
__m256i vmsk = _mm256_set1_epi8(0x0F);
__m256i cnst16 = _mm256_set1_epi32(16);
__m256i shf_lo = _mm256_set_epi8(0x80,0x80,0x80,0x0B, 0x80,0x80,0x80,0x03, 0x80,0x80,0x80,0x0A, 0x80,0x80,0x80,0x02,
0x80,0x80,0x80,0x09, 0x80,0x80,0x80,0x01, 0x80,0x80,0x80,0x08, 0x80,0x80,0x80,0x00);
__m256i shf_hi = _mm256_set_epi8(0x80,0x80,0x80,0x0F, 0x80,0x80,0x80,0x07, 0x80,0x80,0x80,0x0E, 0x80,0x80,0x80,0x06,
0x80,0x80,0x80,0x0D, 0x80,0x80,0x80,0x05, 0x80,0x80,0x80,0x0C, 0x80,0x80,0x80,0x04);
__m128i pshufbcnst = _mm_set_epi8(0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80, 0x0E,0x0C,0x0A,0x08,0x06,0x04,0x02,0x00);
__m256i i_vec = _mm256_setzero_si256();
m0=0;
for (i=0;i<n;i=i+16){
__m128i v = _mm_load_si128((__m128i *)&a[i]); /* Load 16 bytes. */
__m128i msk = _mm_cmpeq_epi8(v,_mm_setzero_si128()); /* Generate 16x8 bit mask. */
msk = _mm_srli_epi64(msk,4); /* Pack 16x8 bit mask to 16x4 bit mask. */
msk = _mm_shuffle_epi8(msk,pshufbcnst); /* Pack 16x8 bit mask to 16x4 bit mask. */
msk = _mm_xor_si128(msk,_mm_set1_epi32(-1)); /* Invert 16x4 mask. */
uint64_t msk64 = _mm_cvtsi128_si64x(msk); /* _mm_popcnt_u64 and _pext_u64 work on 64-bit general-purpose registers, not on simd registers.*/
int p = _mm_popcnt_u64(msk64)>>2; /* p is the number of nonzeros in 16 bytes of a. */
uint64_t cntr = _pext_u64(cntr_const,msk64); /* parallel bits extract. cntr contains p 4-bit integers. The 16 4-bit integers in cntr_const are shuffled to the p 4-bit integers that we want */
/* The next 7 intrinsics unpack these p 4-bit integers to p 32-bit integers. */
__m256i cntr256 = _mm256_set1_epi64x(cntr);
cntr256 = _mm256_srlv_epi64(cntr256,shft);
cntr256 = _mm256_and_si256(cntr256,vmsk);
__m256i cntr256_lo = _mm256_shuffle_epi8(cntr256,shf_lo);
__m256i cntr256_hi = _mm256_shuffle_epi8(cntr256,shf_hi);
cntr256_lo = _mm256_add_epi32(i_vec,cntr256_lo);
cntr256_hi = _mm256_add_epi32(i_vec,cntr256_hi);
_mm256_storeu_si256((__m256i *)&ind[m0],cntr256_lo); /* Note that the stores of iteration i and i+16 may overlap. */
_mm256_storeu_si256((__m256i *)&ind[m0+8],cntr256_hi); /* Array ind has to be large enough to avoid segfaults. At most 16 integers are written more than strictly necessary */
m0 = m0+p;
i_vec = _mm256_add_epi32(i_vec,cnst16);
}
*m=m0;
return 0;
}
__attribute__ ((noinline)) int arr2ind_if(const unsigned char * restrict a, int n, int * restrict ind, int * m){
int i, m0;
m0=0;
for (i=0;i<n;i++){
if (a[i]!=0){
ind[m0]=i;
m0=m0+1;
}
}
*m=m0;
return 0;
}
__attribute__((noinline)) int arr2ind_cmov(const unsigned char * restrict a, int n, int * restrict ind, int * m){
int i, m0;
m0=0;
for (i=0;i<n;i++){
ind[m0]=i;
m0=(a[i]==0)? m0 : m0+1; /* Compiles to cmov instruction. */
}
*m=m0;
return 0;
}
__attribute__ ((noinline)) int print_nonz(const unsigned char * restrict a, const int * restrict ind, const int m){
int i;
for (i=0;i<m;i++) printf("i=%d, ind[i]=%d a[ind[i]]=%u\n",i,ind[i],a[ind[i]]);
printf("\n"); fflush(stdout);
return 0;
}
__attribute__ ((noinline)) int print_chk(const unsigned char * restrict a, const int * restrict ind, const int m){
int i; /* Compute a hash to compare the results of different methods. */
unsigned int chk=0;
for (i=0;i<m;i++){
chk=((chk<<1)|(chk>>31))^(ind[i]);
}
printf("chk = %10X\n",chk);
return 0;
}
int main(int argc, char **argv){
int n, i, m;
unsigned int j, k, d;
unsigned char *a;
int *ind;
double t0,t1;
int meth, nrep;
char txt[30];
sscanf(argv[1],"%d",&n); /* Length of array a. */
n=n>>5; /* Adjust n to a multiple of 32. */
n=n<<5;
sscanf(argv[2],"%u",&d); /* The approximate fraction of nonzeros in a is: d/1024 */
printf("n=%d, d=%u\n",n,d);
a=_mm_malloc(n*sizeof(char),32);
ind=_mm_malloc(n*sizeof(int),32);
/* Generate a pseudo random array a. */
j=73659343;
for (i=0;i<n;i++){
j=j*653+1;
k=(j & 0x3FF00)>>8; /* k is a pseudo random number between 0 and 1023 */
if (k<d){
a[i] = (j&0xFE)+1; /* Set a[i] to nonzero. */
}else{
a[i] = 0;
}
}
/* for (i=0;i<n;i++){if (a[i]!=0){printf("i=%d, a[i]=%u\n",i,a[i]);}} printf("\n"); */ /* Uncomment this line to print the nonzeros in a. */
char txt0[]="arr2ind_movmsk: ";
char txt1[]="arr2ind_pext: ";
char txt2[]="arr2ind_if: ";
char txt3[]="arr2ind_cmov: ";
nrep=10000; /* Repeat a function nrep times to make relatively accurate timings possible. */
/* With nrep=1000000: ./a.out 10016 4 ; ./a.out 10016 48 ; ./a.out 10016 519 */
/* With nrep=10000: ./a.out 1000000 5 ; ./a.out 1000000 52 ; ./a.out 1000000 513 */
printf("nrep = \%d \n\n",nrep);
arr2ind_movmsk(a,n,ind,&m); /* Make sure that the arrays a and ind are read and/or written at least one time before benchmarking. */
for (meth=0;meth<4;meth++){
t0=omp_get_wtime();
switch (meth){
case 0: for(i=0;i<nrep;i++) arr2ind_movmsk(a,n,ind,&m); strcpy(txt,txt0); break;
case 1: for(i=0;i<nrep;i++) arr2ind_pext(a,n,ind,&m); strcpy(txt,txt1); break;
case 2: for(i=0;i<nrep;i++) arr2ind_if(a,n,ind,&m); strcpy(txt,txt2); break;
case 3: for(i=0;i<nrep;i++) arr2ind_cmov(a,n,ind,&m); strcpy(txt,txt3); break;
default: ;
}
t1=omp_get_wtime();
printf("method = %s ",txt);
/* print_chk(a,ind,m); */
printf(" elapsed time = %6.2f\n",t1-t0);
}
print_nonz(a, ind, 2); /* Do something with the results */
printf("density = %f %% \n\n",((double)m)/((double)n)*100); /* Actual nonzero density of array a. */
/* print_nonz(a, ind, m); */ /* Uncomment this line to print the indices of the nonzeros. */
return 0;
}
/*
With nrep=1000000:
./a.out 10016 4 ; ./a.out 10016 4 ; ./a.out 10016 48 ; ./a.out 10016 48 ; ./a.out 10016 519 ; ./a.out 10016 519
With nrep=10000:
./a.out 1000000 5 ; ./a.out 1000000 5 ; ./a.out 1000000 52 ; ./a.out 1000000 52 ; ./a.out 1000000 513 ; ./a.out 1000000 513
*/
코드는 N의 어레이 크기로 시험 하였다 = 10016 (데이터가 L1 캐시에 맞는) 및 N = 1000000의 다른 제로 밀도 약 0.5 %, 5 % 및 50 %이다. 정확한 타이밍을 위해 함수는 각각 1000000 및 10000 번 호출되었습니다.
Time in seconds, size n=10016, 1e6 function calls. Intel core i5-6500
0.53% 5.1% 50.0%
arr2ind_movmsk: 0.27 0.53 4.89
arr2ind_pext: 1.44 1.59 1.45
arr2ind_if: 5.93 8.95 33.82
arr2ind_cmov: 6.82 6.83 6.82
Time in seconds, size n=1000000, 1e4 function calls.
0.49% 5.1% 50.1%
arr2ind_movmsk: 0.57 2.03 5.37
arr2ind_pext: 1.47 1.47 1.46
arr2ind_if: 5.88 8.98 38.59
arr2ind_cmov: 6.82 6.81 6.81
벡터화 된 루프 빨리 스칼라 루프보다 이러한 예에서
arr2ind_movmsk
의 성능은
a
의 밀도에 많이 달려 있습니다. 밀도가 충분히 작 으면 일 때
arr2ind_pext
보다 빠릅니다. 손익 분기점은 배열 크기
n
에 따라 다릅니다. 'arr2ind_if'함수는 명확하게 50 %가 아닌 밀도에서 분기 예측에 실패했습니다.
직접적으로는 아니지만'pmovmskb'를 일반 레지스터에'pcmpeqb '할 수 있고'bsf' (그리고 두 번째 등등, 잘하면 너무 많지는 않음)로 첫 번째 인덱스를 추출합니다. – harold
당신은 단지 'SIMD'보다 구체적 일 필요가 있습니다 - 당신이 목표로 삼고있는 아키텍처는 무엇입니까?x86, ARM, PowerPC, POWER 및 일부 GPGPU에는 모두 다른 SIMD 확장이 있습니다. 또한 x86에는 MMX, SSE, SSE2, SSE3, SSSE3, SSE4, AVX, AVX2 등의 여러 SIMD 확장이 있습니다 (AVX2에는이 컨텍스트에서 유용 할 수있는 SIMD 명령어가 있음을 참고하십시오). –
@Paul R 죄송합니다. 내 질문을 편집 - AVX2 허용됩니다. –