SSE _mm_load_pd работает в то время как _mm_store_pd segfaults

Question

Я пытаюсь изучить веревки SSE-intrinsics в C. У меня есть часть кода, где я загружаю двухкомпонентный вектор двойных данных, добавляю что-то к нему, а затем пытаюсь его сохранить вернуться в память. Все работает: я могу загружать свои данные в регистры SEE, я могу работать с моими данными в этих регистрах SSE, но в тот момент, когда я пытаюсь записать обработанные данные обратно в исходный массив (где я читаю свои данные из первое место!) Я получаю ошибку сегментации.

Может кто угодно посоветуйте мне по этому вопросу - это сводит меня с ума.

double res[2] __attribute__((aligned(16))); 

for(int k=0; k<n; k++){ 
int i=0; 
for(; i+1<n; i+=2) 
    { 
    __m128d cik = _mm_load_pd(&C[i+k*n]); 
    int j = 0; 
    for(; j+1<n; j+=2) 
     { 
     __m128d aTij = _mm_load_pd(&A_T[j+i*n]); 
     __m128d bjk = _mm_load_pd(&B[j+k*n]); 
     __m128d dotpr = _mm_dp_pd(aTij, bjk,2); 
     cik = _mm_add_pd(cik, dotpr); 
     } 
    _mm_store_pd(res, cik); 
    //C[i+k*n] = res[0]; 
    } 
} 

Как я уже говорил выше, все работает в этом коде для где я храню свои результаты обратно в этот одномерный массив «C», где я читал мои данные в первую очередь, за исключением. То есть, когда я раскомментировать знаки перед

//C[i+k*n] = res[0]; 

Я получаю ошибку сегментации.

Как возможно, что я могу читать с C с выровненной версией памяти _mm_load_pd (так что C должен быть выровнен в память!), В то время как запись на него не работает? «C» должен быть выровнен, и, как вы видите, «res» также должен быть выровнен.

Отказ от ответственности: Мой исходный код чтения

_mm_store_pd(&C[i+k*n], cik); 

который также произвел ошибку сегментации, и я начал внедрять «Рес» с явным выравниванием в моей попытки решить эту проблему.

Добавление

A, B, C объявляются следующим образом:

buf = (double*) malloc (3 * nmax * nmax * sizeof(double)); 
double* A = buf + 0; 
double* B = A + nmax*nmax; 
double* C = B + nmax*nmax; 

Попытка решения с posix_memalign

В попытке решить проблему ошибки сегментации при записи на оригинальный одномерный массив, теперь я использую буферы для соответствующих матриц. Тем не менее, это все еще segfauls при попытке написать обратно C_buff!

double res[2] __attribute__((aligned(16))); 

double * A_T; 
posix_memalign((void**)&A_T, 16, n*n*sizeof(double)); 

double * B_buff; 
posix_memalign((void**)&B_buff, 16, n*n*sizeof(double)); 

double * C_buff; 
posix_memalign((void**)&C_buff, 16, n*n*sizeof(double)); 

for(int y=0; y<n; y++) 
    for(int x=0; x<n; x++) 
    A_T[x+y*n] = A[y+x*n]; 

for(int x=0; x<n; x++) 
    for(int y=0; y<n; y++) 
    B_buff[y+x*n] = B[y+x*n]; 

for(int x=0; x<n; x++) 
    for(int y=0; y<n; y++) 
    C_buff[y+x*n] = C[y+x*n]; 

for(int k=0; k<n; k++){ 
    int i=0; 
    for(; i+1<n; i+=2) 
    { 
     __m128d cik = _mm_load_pd(&C_buff[i+k*n]); 
     int j = 0; 
     for(; j+1<n; j+=2) 
     { 
      __m128d aTij = _mm_load_pd(&A_T[j+i*n]); 
      __m128d bjk = _mm_load_pd(&B_buff[j+k*n]); 
      __m128d dotpr = _mm_dp_pd(aTij, bjk,2); 
      cik = _mm_add_pd(cik, dotpr); 
     } 
     _mm_store_pd(&C_buff[i+k*n], cik); 

    //_mm_store_pd(res, cik); 
     //C_buff[i+k*n] = res[0]; 
    //C_buff[i+1+k*n] = res[1]; 
    } 
} 

Answer 1

Даже с __attribute__((aligned(32))), я получаю ту же ошибку (это было 50% шанс неправильного alingment). Тогда я использовал следующую функцию, чтобы получить% 100 шанса выравнивания (а должно быть мощность два):

void * malloc_float_align(size_t n, unsigned int a/*alignment*/, float *& output) 
{ 
    void * adres=NULL; 
    void * adres2=NULL; 
    adres=malloc(n*sizeof(float)+a); 
    size_t adr=(size_t)adres; 
    size_t adr2=adr+a-(adr&(a-1u)); // a valid address for a alignment 
    adres2=(void *) adr2; 
    output=(float *)adres2; 
    return adres;    //pointer to be used in free() 
} 

Тогда использование в основном:

int main() 
{ 


    float * res=NULL; 
    void * origin=malloc_float_align(1024,32u,res); 
    //use res for sse/avx 
    free(origin); // actual allocation is more than 1024 elements 
    return 0; 
} 

Ofcourse это в C++, так что вы можете заставить работать только с изменением стиля параметра функции.

Answer 2

Простой трюк будет выполнить ASSERT и посмотреть, если это вызывает:

ASSERT(((size_t)(&C_buff[i+k*n]) & 0xF) == 0); 

ASSERT будет срабатывать, когда адрес не SSE выровнены. 64-битные сборки должны обеспечивать выравнивание по 16B по умолчанию. Если вы планируете иметь 32-битный код, используйте одну из вышеперечисленных функций align_malloc. Вам нужно будет использовать соответствующий align_free или вы потерпите крах.

Answer 3

При удалении _mm_store_pd(&C_buff[i+k*n], cik); весь цикл оптимизирован и удален. Компилятор выводит, что цикл for for не приводит к какой-либо значимой работе и удаляет ее. Вот почему вы больше не получаете дефектов сегментации.
Я уверен, что ошибка сегментации связана с размером массива. Рассмотрим эту простую программу, основанную на вашем примере:

#include <stdio.h> 
#include "emmintrin.h" 

int main(){ 
int n = 15; 
int y,x,k,i,j; 

double * A; 
posix_memalign((void**)&A, 16, n*n*sizeof(double)); 

double * B; 
posix_memalign((void**)&B, 16, n*n*sizeof(double)); 

double * C; 
posix_memalign((void**)&C, 16, n*n*sizeof(double)); 

for(y=0; y<n; y++) 
    for(x=0; x<n; x++) 
    A[x+y*n] = 0.1; 

for(x=0; x<n; x++) 
    for(y=0; y<n; y++) 
    B[y+x*n] = 0.1; 

for(x=0; x<n; x++) 
    for(y=0; y<n; y++) 
    C[y+x*n] = 0.1; 

for(k=0; k<n; k++){ 
    i=0; 
    for(; i+1<n; i+=2) 
    { 
     __m128d cik = _mm_load_pd(&C[i+k*n]); 
     j = 0; 
     for(; j+1<n; j+=2) 
     { 
      __m128d aTij = _mm_load_pd(&A[j+i*n]); 
      __m128d bjk = _mm_load_pd(&B[j+k*n]); 
      __m128d dotpr = _mm_add_pd(aTij, bjk); 
      cik = _mm_add_pd(cik, dotpr); 
     } 
     _mm_store_pd(&C[i+k*n], cik); 
    } 
} 
printf("C[15]: %f\n", C[15]); 
printf("C[14]: %f\n", C[14]); 

Это дает ошибку сегментации, потому что п нечетно. Теперь измените n = 15 на n = 16, и все будет работать, как ожидалось. Следовательно, добавление ваших массивов к четному числу (или даже лучше к размеру строки кэша -> 64 байта == 8 элементов DP или 16 элементов SP) предотвратит такие проблемы и приведет к лучшей производительности.