Выбранный подход к решению: Считываем в регистры 3 массива по 4 элемента, один со сдвигом влево (копируем с i-1 элемента), второй с нулевым сдвигом и третий со сдвигом вправо (копируем с i+1 элемента). Делаем вычисления. Копируем обратно. Программа работает, но sse версия работает много медленнее чем обычное решение.
Подозреваю, что проблема в использовании _mm_set_ps для копирования в регистры и организации обратного копирования из регистров.
Вопрос: Как можно улучшить данный код с целью ускорения работы?
Более конкретный вопрос: Как можно организовать копирование начиная с любого i-го (не кратного 4) элемента массива в регистр и наоборот? Без использования _mm_set_ps и костыля (for (int j=0;j<4;j=j++) Y3[i+j]=d2Ydx2[j];)
CODE
// 1st_SSE.cpp
//
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include <iostream>
using std::cout;
using std::cin;
using std::endl;
int main()
{
const int N=16*1024+1;
float *d2Ydx2=(float*)malloc(4*sizeof(float));
__m128 Ysse,YRsse,YMsse,YLsse;
__m128 rrs=_mm_set1_ps(-2.0f);
__m128 dX2s=_mm_set1_ps((N*N)/(6.28319f*6.28319f));
__m128 *d2Ydx2SSE=(__m128*) d2Ydx2;
float *Y1=(float*)malloc(N*sizeof(float));
float *Y=(float*)malloc(N*sizeof(float));
float *Y3=(float*)malloc(N*sizeof(float));
__m128 *Y2=(__m128*) Y;
timespec ts_beg, ts_end;
float rr,dX2;
const int SSEN=(N-1)/4;
dX2=(N*N)/(6.28319f*6.28319f);
rr=-2.0f;
//задаем исходную функцию
for (int i=0;i<=N;i++)
{
Y[i]=sin(i*(6.28319f/N));
Y3[i]=0.f;
Y1[i]=0.f;
}
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_beg);
// решение в лоб
for (int i=1;i<=N-1;i++)
{
Y1[i]=(Y[i-1]+rr*Y[i]+Y[i+1])*dX2;
};
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_end);
float time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec)/1e9;
cout << "Time required = " << time <<endl;
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_beg);
for (int i=1;i<=N-1;i=i+4)
{
Ysse=_mm_set1_ps(0.f);
//копируем в регистры
YLsse=_mm_set_ps(Y[i+2],Y[i+1],Y[i],Y[i-1]);
YMsse=_mm_set_ps(Y[i+3],Y[i+2],Y[i+1],Y[i]);
YRsse=_mm_set_ps(Y[i+4],Y[i+3],Y[i+2],Y[i+1]);
//вычисляем производную
Ysse=_mm_mul_ps(YMsse,rrs);
Ysse=_mm_add_ps(Ysse,YLsse);
Ysse=_mm_add_ps(Ysse,YRsse);
_mm_store_ps(d2Ydx2, _mm_mul_ps(Ysse,dX2s));
//копируем обратно
for (int j=0;j<4;j=j++)
{
Y3[i+j]=d2Ydx2[j];
}
};
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_end);
time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec)/1e9;
cout << "Time required = " << time <<endl;
for (int i=0;i<=N;i=i++)
{
if (Y1[i]!=Y3[i]) cout<<i<<"\t"<<Y1[i]<<"\t"<<Y3[i]<<endl;
};
free(d2Ydx2);
free(Y1);
free(Y);
free(Y3);
return 0;
}
//
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include <iostream>
using std::cout;
using std::cin;
using std::endl;
int main()
{
const int N=16*1024+1;
float *d2Ydx2=(float*)malloc(4*sizeof(float));
__m128 Ysse,YRsse,YMsse,YLsse;
__m128 rrs=_mm_set1_ps(-2.0f);
__m128 dX2s=_mm_set1_ps((N*N)/(6.28319f*6.28319f));
__m128 *d2Ydx2SSE=(__m128*) d2Ydx2;
float *Y1=(float*)malloc(N*sizeof(float));
float *Y=(float*)malloc(N*sizeof(float));
float *Y3=(float*)malloc(N*sizeof(float));
__m128 *Y2=(__m128*) Y;
timespec ts_beg, ts_end;
float rr,dX2;
const int SSEN=(N-1)/4;
dX2=(N*N)/(6.28319f*6.28319f);
rr=-2.0f;
//задаем исходную функцию
for (int i=0;i<=N;i++)
{
Y[i]=sin(i*(6.28319f/N));
Y3[i]=0.f;
Y1[i]=0.f;
}
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_beg);
// решение в лоб
for (int i=1;i<=N-1;i++)
{
Y1[i]=(Y[i-1]+rr*Y[i]+Y[i+1])*dX2;
};
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_end);
float time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec)/1e9;
cout << "Time required = " << time <<endl;
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_beg);
for (int i=1;i<=N-1;i=i+4)
{
Ysse=_mm_set1_ps(0.f);
//копируем в регистры
YLsse=_mm_set_ps(Y[i+2],Y[i+1],Y[i],Y[i-1]);
YMsse=_mm_set_ps(Y[i+3],Y[i+2],Y[i+1],Y[i]);
YRsse=_mm_set_ps(Y[i+4],Y[i+3],Y[i+2],Y[i+1]);
//вычисляем производную
Ysse=_mm_mul_ps(YMsse,rrs);
Ysse=_mm_add_ps(Ysse,YLsse);
Ysse=_mm_add_ps(Ysse,YRsse);
_mm_store_ps(d2Ydx2, _mm_mul_ps(Ysse,dX2s));
//копируем обратно
for (int j=0;j<4;j=j++)
{
Y3[i+j]=d2Ydx2[j];
}
};
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_end);
time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec)/1e9;
cout << "Time required = " << time <<endl;
for (int i=0;i<=N;i=i++)
{
if (Y1[i]!=Y3[i]) cout<<i<<"\t"<<Y1[i]<<"\t"<<Y3[i]<<endl;
};
free(d2Ydx2);
free(Y1);
free(Y);
free(Y3);
return 0;
}