Форум Академгородка, Новосибирск > вторая производная с использованием SSE
Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: вторая производная с использованием SSE
Форум Академгородка, Новосибирск > Компьютеры и сети > Программирование
Texnik
Задача: Таблично задана функция, необходимо вычислить вторую производную методом конечных разностей.
Выбранный подход к решению: Считываем в регистры 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;
}


Nemo
Вместо _mm_set_ps надо использовать _mm_loadu_ps с начальным адресом, чтобы сделать одно большое чтение, а не 4 маленьких. Ну на самом деле можно сделать ещё хитрее и обойтись вообще одним выровненным чтением на 4 итерации: надо сначала зачитать выровнено два полных регистра - там будет 3 лишних элемента, операциями shufps и unpackps можно перепаковать данные в нужном порядке и делать так ещё 3 итерации после, потом загрузить очередной регистр и т.д...

Ну и собственно _mm_loadu_ps - это ответ на более конкретный вопрос.
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Русская версия IP.Board © 2001-2024 IPS, Inc.