Comments 33
Провозившись с реализацией fft на fpga, я было подумал какой это гимор, на процах ведь проще. Однако после этой статьи задумался - даже описание моей реализации займет меньше текста)
Здесь 90% статьи — это анализ кода на Си и ассемблере. При первом прочтении разделы Пишем векторный Reverse, Пишем векторный Stage и Пишем LastStage можно прокрутить.
Без этих разделов получится статья на 20 минут — вошли и вышли.
А если сравнить с тривиальной реализацией?
vfft.c
#include <math.h>
typedef double vfft_flt;
typedef struct vfft_t {
int p,n; // n=2^p
vfft_flt* w; // [n]
// w[2*k ]=cos(2*PI*k/n)
// w[2*k+1]=sin(2*PI*k/n)
void *ctx;
void (*scale)(void* ctx,vfft_flt a); // x[]*=a
void (*swap )(void* ctx,int i,int j); // swap(x[i],x[j])
void (*twist)(void* ctx,int i,int j,vfft_flt *w);
// complex *x, t=x[i]-x[j]; x[i]=x[i]+x[j]; x[j]=t*w;
} vfft_t;
static void vfft_twist_back(vfft_t *fft) {
int n=fft->n,p=fft->p,a=n,c,d,f,g,i,j;
vfft_flt *w=fft->w;
for(c=1;c<=p;c++) { d=a; a>>=1;
for(f=0;f<a;f++) {
for(g=0;g<n;g+=d) { i=g+f; j=i+a;
fft->twist(fft->ctx,i,j,w+(f<<c));
}
}
}
}
static void vfft_twist(vfft_t *fft) {
int n=fft->n,p=fft->p,a=n,c,d,f,g,i,j;
vfft_flt *w=fft->w, t[2];
for(c=1;c<=p;c++) { d=a; a>>=1;
for(f=0;f<a;f++) {
for(g=0;g<n;g+=d) { i=g+f; j=i+a;
t[0]=w[f<<c]; t[1]=-w[(f<<c)+1];
fft->twist(fft->ctx,i,j,t);
}
}
}
}
static void vfft_reorder(vfft_t *fft) {
int p=fft->n>>1,q=fft->n-1,j=0,i,k;
for(i=0;i<q;i++,j+=k) {
if (i<j) fft->swap(fft->ctx,i,j);
k=p; while(k<=j) { j-=k; k>>=1; }
}
}
void vfft_mkw(vfft_t *fft) {
int i, n2=fft->n/2; vfft_flt *w=fft->w, k=(vfft_flt)M_PI/n2;
w[0]=1; w[1]=0;
for(i=1;i<n2;i++) {
w[2*i ]=cos(k*i);
w[2*i+1]=sin(k*i);
}
}
void vfft_run(vfft_t *fft,int forward) {
if (forward) vfft_twist(fft); else vfft_twist_back(fft);
vfft_reorder(fft);
if (forward) fft->scale(fft->ctx,(vfft_flt)1/fft->n);
}
//---------------------------------------------------------------------------
typedef struct vfft_vector_xy_t { int n; vfft_flt *x, *y; } vfft_vector_xy_t;
void vfft_vector_xy_twist(void* ctx,int i,int j,vfft_flt *w) {
vfft_vector_xy_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, *y=p->y, t[2];
t[0]=x[i]-x[j];
t[1]=y[i]-y[j];
x[i]=x[i]+x[j];
y[i]=y[i]+y[j];
x[j]=t[0]*w[0]-t[1]*w[1];
y[j]=t[1]*w[0]+t[0]*w[1];
}
void vfft_vector_xy_swap(void* ctx,int i,int j) {
vfft_vector_xy_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, *y=p->y, t;
t=x[i];x[i]=x[j];x[j]=t;
t=y[i];y[i]=y[j];y[j]=t;
}
void vfft_vector_xy_scale(void* ctx,vfft_flt a) {
vfft_vector_xy_t *p; *(void**)&p=ctx; int i,n=p->n;
vfft_flt *x=p->x, *y=p->y;
for(i=0;i<n;i++) { x[i]*=a; y[i]*=a; }
}
//---------------------------------------------------------------------------
typedef struct vfft_vector2_t { int n; vfft_flt *x; } vfft_vector2_t;
void vfft_vector2_twist(void* ctx,int i,int j,vfft_flt *w) {
vfft_vector2_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, t[2]; i*=2; j*=2;
t[0]=x[i ]-x[j ];
t[1]=x[i+1]-x[j+1];
x[i ]+=x[j ];
x[i+1]+=x[j+1];
x[j ]=t[0]*w[0]-t[1]*w[1];
x[j+1]=t[1]*w[0]+t[0]*w[1];
}
void vfft_vector2_swap(void* ctx,int i,int j) {
vfft_vector2_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, t;
i*=2;j*=2;
t=x[i];x[i]=x[j];x[j]=t;
i++;j++;
t=x[i];x[i]=x[j];x[j]=t;
}
void vfft_vector2_scale(void* ctx,vfft_flt a) {
vfft_vector2_t *p; *(void**)&p=ctx; int i,n=2*p->n;
vfft_flt *x=p->x;
for(i=0;i<n;i++) x[i]*=a;
}
//---------------------------------------------------------------------------
typedef struct vfft_image_t { int n,m ; vfft_flt **x, **y; } vfft_image_t;
void vfft_image_twist(void* ctx,int i,int j,vfft_flt *w) {
vfft_image_t *p; *(void**)&p=ctx; int k,m=p->m; vfft_flt t[2];
vfft_flt *xi=p->x[i], *yi=p->y[i], *xj=p->x[j], *yj=p->y[j];
for(k=0;k<m;k++) {
t[0]=xi[k]-xj[k];
t[1]=yi[k]-yj[k];
xi[k]+=xj[k];
yi[k]+=yj[k];
xj[k]=t[0]*w[0]-t[1]*w[1];
yj[k]=t[1]*w[0]+t[0]*w[1];
}
}
void vfft_image_swap(void* ctx,int i,int j) {
vfft_image_t *p; *(void**)&p=ctx; int k,m=p->m;
vfft_flt *xi=p->x[i], *yi=p->y[i], *xj=p->x[j], *yj=p->y[j], t;
for(k=0;k<m;k++) {
t=xi[k];xi[k]=xj[k];xj[k]=t;
t=yi[k];yi[k]=yj[k];yj[k]=t;
}
}
void vfft_image_scale(void* ctx,vfft_flt a) {
vfft_image_t *p; *(void**)&p=ctx; int i,n=p->n, k,m=p->m;
for(i=0;i<n;i++) {
vfft_flt *xi=p->x[i], *yi=p->y[i];
for(k=0;k<m;k++) {
xi[k]*=a; yi[k]*=a;
}
}
}
//---------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
int main(int argc,char* argv) {
enum { P=4, N=1<<P }; int i;
double x[N],y[N],w[N]; vfft_t vfft[1]; vfft_vector_xy_t v[1];
v->n=N;
v->x=x;
v->y=y;
for(i=0;i<N;i++) { x[i]=cos(i*2*M_PI*3/16+0.7); y[i]=0; }
vfft->p=P;
vfft->n=N;
vfft->w=w;
vfft->ctx=v;
vfft->scale=vfft_vector_xy_scale;
vfft->twist=vfft_vector_xy_twist;
vfft->swap =vfft_vector_xy_swap;
vfft_mkw(vfft);
printf("x: "); for(i=0;i<N;i++) printf(" %6.3f",x[i]); printf("\n");
printf("y: "); for(i=0;i<N;i++) printf(" %6.3f",y[i]); printf("\n");
vfft_run(vfft,1);
printf("freq\n");
printf("x: "); for(i=0;i<N;i++) printf(" %6.3f",x[i]); printf("\n");
printf("y: "); for(i=0;i<N;i++) printf(" %6.3f",y[i]); printf("\n");
vfft_run(vfft,0);
printf("time\n");
printf("x: "); for(i=0;i<N;i++) printf(" %6.3f",x[i]); printf("\n");
printf("y: "); for(i=0;i<N;i++) printf(" %6.3f",y[i]); printf("\n");
return 0;
}
Изменил в начале файла typedef double vfft_flt на typedef float vfft_flt.
Дописал в функцию main замеры скорости, как у меня в _speedtest.
Замеры делаю для функции vfft_twist (проход по всем Stage).
vfft.c
#include <math.h>
typedef /*double*/float vfft_flt;
typedef struct vfft_t {
int p,n; // n=2^p
vfft_flt* w; // [n]
// w[2*k ]=cos(2*PI*k/n)
// w[2*k+1]=sin(2*PI*k/n)
void *ctx;
void (*scale)(void* ctx,vfft_flt a); // x[]*=a
void (*swap )(void* ctx,int i,int j); // swap(x[i],x[j])
void (*twist)(void* ctx,int i,int j,vfft_flt *w);
// complex *x, t=x[i]-x[j]; x[i]=x[i]+x[j]; x[j]=t*w;
} vfft_t;
static void vfft_twist_back(vfft_t *fft) {
int n=fft->n,p=fft->p,a=n,c,d,f,g,i,j;
vfft_flt *w=fft->w;
for(c=1;c<=p;c++) { d=a; a>>=1;
for(f=0;f<a;f++) {
for(g=0;g<n;g+=d) { i=g+f; j=i+a;
fft->twist(fft->ctx,i,j,w+(f<<c));
}
}
}
}
static void vfft_twist(vfft_t *fft) {
int n=fft->n,p=fft->p,a=n,c,d,f,g,i,j;
vfft_flt *w=fft->w, t[2];
for(c=1;c<=p;c++) { d=a; a>>=1;
for(f=0;f<a;f++) {
for(g=0;g<n;g+=d) { i=g+f; j=i+a;
t[0]=w[f<<c]; t[1]=-w[(f<<c)+1];
fft->twist(fft->ctx,i,j,t);
}
}
}
}
static void vfft_reorder(vfft_t *fft) {
int p=fft->n>>1,q=fft->n-1,j=0,i,k;
for(i=0;i<q;i++,j+=k) {
if (i<j) fft->swap(fft->ctx,i,j);
k=p; while(k<=j) { j-=k; k>>=1; }
}
}
void vfft_mkw(vfft_t *fft) {
int i, n2=fft->n/2; vfft_flt *w=fft->w, k=(vfft_flt)M_PI/n2;
w[0]=1; w[1]=0;
for(i=1;i<n2;i++) {
w[2*i ]=cos(k*i);
w[2*i+1]=sin(k*i);
}
}
void vfft_run(vfft_t *fft,int forward) {
if (forward) vfft_twist(fft); else vfft_twist_back(fft);
vfft_reorder(fft);
if (forward) fft->scale(fft->ctx,(vfft_flt)1/fft->n);
}
//---------------------------------------------------------------------------
typedef struct vfft_vector_xy_t { int n; vfft_flt *x, *y; } vfft_vector_xy_t;
void vfft_vector_xy_twist(void* ctx,int i,int j,vfft_flt *w) {
vfft_vector_xy_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, *y=p->y, t[2];
t[0]=x[i]-x[j];
t[1]=y[i]-y[j];
x[i]=x[i]+x[j];
y[i]=y[i]+y[j];
x[j]=t[0]*w[0]-t[1]*w[1];
y[j]=t[1]*w[0]+t[0]*w[1];
}
void vfft_vector_xy_swap(void* ctx,int i,int j) {
vfft_vector_xy_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, *y=p->y, t;
t=x[i];x[i]=x[j];x[j]=t;
t=y[i];y[i]=y[j];y[j]=t;
}
void vfft_vector_xy_scale(void* ctx,vfft_flt a) {
vfft_vector_xy_t *p; *(void**)&p=ctx; int i,n=p->n;
vfft_flt *x=p->x, *y=p->y;
for(i=0;i<n;i++) { x[i]*=a; y[i]*=a; }
}
//---------------------------------------------------------------------------
typedef struct vfft_vector2_t { int n; vfft_flt *x; } vfft_vector2_t;
void vfft_vector2_twist(void* ctx,int i,int j,vfft_flt *w) {
vfft_vector2_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, t[2]; i*=2; j*=2;
t[0]=x[i ]-x[j ];
t[1]=x[i+1]-x[j+1];
x[i ]+=x[j ];
x[i+1]+=x[j+1];
x[j ]=t[0]*w[0]-t[1]*w[1];
x[j+1]=t[1]*w[0]+t[0]*w[1];
}
void vfft_vector2_swap(void* ctx,int i,int j) {
vfft_vector2_t *p; *(void**)&p=ctx;
vfft_flt *x=p->x, t;
i*=2;j*=2;
t=x[i];x[i]=x[j];x[j]=t;
i++;j++;
t=x[i];x[i]=x[j];x[j]=t;
}
void vfft_vector2_scale(void* ctx,vfft_flt a) {
vfft_vector2_t *p; *(void**)&p=ctx; int i,n=2*p->n;
vfft_flt *x=p->x;
for(i=0;i<n;i++) x[i]*=a;
}
//---------------------------------------------------------------------------
typedef struct vfft_image_t { int n,m ; vfft_flt **x, **y; } vfft_image_t;
void vfft_image_twist(void* ctx,int i,int j,vfft_flt *w) {
vfft_image_t *p; *(void**)&p=ctx; int k,m=p->m; vfft_flt t[2];
vfft_flt *xi=p->x[i], *yi=p->y[i], *xj=p->x[j], *yj=p->y[j];
for(k=0;k<m;k++) {
t[0]=xi[k]-xj[k];
t[1]=yi[k]-yj[k];
xi[k]+=xj[k];
yi[k]+=yj[k];
xj[k]=t[0]*w[0]-t[1]*w[1];
yj[k]=t[1]*w[0]+t[0]*w[1];
}
}
void vfft_image_swap(void* ctx,int i,int j) {
vfft_image_t *p; *(void**)&p=ctx; int k,m=p->m;
vfft_flt *xi=p->x[i], *yi=p->y[i], *xj=p->x[j], *yj=p->y[j], t;
for(k=0;k<m;k++) {
t=xi[k];xi[k]=xj[k];xj[k]=t;
t=yi[k];yi[k]=yj[k];yj[k]=t;
}
}
void vfft_image_scale(void* ctx,vfft_flt a) {
vfft_image_t *p; *(void**)&p=ctx; int i,n=p->n, k,m=p->m;
for(i=0;i<n;i++) {
vfft_flt *xi=p->x[i], *yi=p->y[i];
for(k=0;k<m;k++) {
xi[k]*=a; yi[k]*=a;
}
}
}
//---------------------------------------------------------------------------
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
uint64_t get_clock_count()
{
uint64_t dst;
#pragma asm_inline
asm ("rrd %%clkr, %0" : "=r" (dst));
return dst;
}
int main(int argc,char** argv) {
if(argc < 1 + 2 || atoi(argv[1]) <= 0 || atoi(argv[2]) <= 0)
{
printf("Usage: %s <bit_count> <repeats>\n", argv[0]);
printf("Example: %s 12 10000 \n", argv[0]);
printf("<bit_count> and <repeats> must be positive\n");
return 1;
}
int bit_count = atoi(argv[1]);
int repeats = atoi(argv[2]);
int data_count = 1 << bit_count;
if(argc < 1 + 3)
{
printf("Количество элементов массива (float complex): %d\n", data_count);
printf("Количество повторов генерации входных данных: %d\n", repeats);
}
int P = bit_count;
int N = data_count;
/*enum { P=4, N=1<<P };*/ int i;
/*double*/vfft_flt x[N],y[N],w[N]; vfft_t vfft[1]; vfft_vector_xy_t v[1];
v->n=N;
v->x=x;
v->y=y;
for(i=0;i<N;i++) { x[i]=cos(i*2*M_PI*3/16+0.7); y[i]=0; }
vfft->p=P;
vfft->n=N;
vfft->w=w;
vfft->ctx=v;
vfft->scale=vfft_vector_xy_scale;
vfft->twist=vfft_vector_xy_twist;
vfft->swap =vfft_vector_xy_swap;
vfft_mkw(vfft);
// printf("x: "); for(i=0;i<N;i++) printf(" %6.3f",x[i]); printf("\n");
// printf("y: "); for(i=0;i<N;i++) printf(" %6.3f",y[i]); printf("\n");
//
// vfft_run(vfft,1);
// printf("freq\n");
// printf("x: "); for(i=0;i<N;i++) printf(" %6.3f",x[i]); printf("\n");
// printf("y: "); for(i=0;i<N;i++) printf(" %6.3f",y[i]); printf("\n");
//
// vfft_run(vfft,0);
// printf("time\n");
// printf("x: "); for(i=0;i<N;i++) printf(" %6.3f",x[i]); printf("\n");
// printf("y: "); for(i=0;i<N;i++) printf(" %6.3f",y[i]); printf("\n");
// return 0;
puts(" (Алгоритм) (Время, мкс) (Такты) (Такты/мкс) (Байты/такт)");
{
struct timespec t0, t1;
clock_gettime(CLOCK_REALTIME, &t0);
uint64_t ticks0 = get_clock_count();
for(int rep = 0; rep < repeats; ++rep)
vfft_twist(vfft);
uint64_t ticks1 = get_clock_count();
clock_gettime(CLOCK_REALTIME, &t1);
double usec = (t1.tv_sec - t0.tv_sec)*1000000 + (t1.tv_nsec - t0.tv_nsec)/1000;
double ticks = ticks1 - ticks0;
double avg_usec = usec /repeats;
double avg_ticks = ticks/repeats;
double byte_count = data_count*sizeof(vfft_flt)*2;
printf("%3d. %-44s" "%9.0lf %12.0lf %9.2lf %9.2lf\n",
0 + 1, "vfft",
avg_usec, avg_ticks, ticks/usec, byte_count/avg_ticks
);
}
}
Замеры скорости
Руками что ли все?
В год когда Навье-Стокс сдался на милость агентам. Ручное написание кода выглядит мягко говоря странным
Изначальная цель — изучение ассемблера процессора Эльбрус.
Я здесь пишу не столько на Си, сколько на ассемблере.
У Вас очень хороший результат для 16К, но просто обвал для 64К. ИИ убеждает меня , что этого не избежать. Можете ли Вы подтвердить, что разбиение БПФ 64 К на четыре части по 16К + заключительная общая стадия, не даст улучшения скорости для 8СВ?
Я понимаю, как это сделать на классическом Stage, но пока не понимаю, можно ли такое сделать на «адаптированном для Эльбруса».
В классическом Stage обработка происходит иерархически.
В «адаптированном» постоянно всё перемешивается.
Последние две статьи я пишу не на Эльбрус-8СВ (elbrus v5), а на e2c3 (elbrus-v6, 1600 МГц).
КДПВ из статей
Спасибо за уточнение. Немного проштудировал Е2С3. Оказывается там только 2 М кэша L2 для одного ядра. При 64К мы имеем возможность хранить в нем 4 массива. И при этом им будет мешать еще и код программы.
массив данных = (64К чисел) * (8 Байт/число) = 512 КБ
+ второй массив = 512 КБ (алгоритм не «in-place»)
+ коэффициенты = 256 КБ (Radix-2) или 384 КБ (Radix-4)
Вроде как все нормально с количеством массивов. Я посмотрел Ваш исходник. Входной массив после реверса идет в выходной, а далее два массива чередуются. По логике 3 массива констант должны быть расположены друг за другом и образовывать один большой. Я бы создал один большой и поделил его на 3 :) Пока не видно причины почему бы и БПФ 64К не делать так же быстро как и 16К. Размер кэша как бы позволяет. Хорошо бы если бы Вы далее продолжили оптимизацию не малых размеров БПФ, а больших. Все-таки процессор у Вас серьезный и задачи он должен решать соответствующие :)
Обвал скорости на 216 виден уже в Stage.
Правильнее всего будет изучить, почему Stage не влезает в кэш даже на Radix-2 Vector-2, где коэффициентов меньше.
ИИ советует :) для 64К использовать не более 3-х массивов т.к. 4-й канал может понадобиться для иных нужд. Также советует 3 массива констант обьединить в структуру чтобы гарантировать их последовательное расположение.
Проверил stage_radix2_vector2:
Выделил память одним куском.
Положил сначала data_in, потом сразу coef, остальное — ещё дальше.
int data_count2 = data_count/2;
int align = 4*1024*1024;
int size = data_count2 + data_count2 + data_count2/2;
void *mem0 = malloc(size*sizeof(myComplex2) + (align-1));
void *mem = (void*)(((uint64_t)mem0 + (align-1)) & (~(align-1)));
myComplex2 *data_in = mem;
myComplex2 *coef = (void*)data_in + data_count2 *sizeof(myComplex2);
myComplex2 *data_out = (void*)coef + data_count2/2*sizeof(myComplex2);Лучше не стало:
Я вспомнил, что делал тест скоростей чтения для elbrus-v5 и elbrus-v6.
Хотел опубликовать в отдельной статье после завершения работы над FFT.
Задача: сделать xor всех байтов в массиве uint8_t.
От этого теста elbrus-v5 ещё зависал в какой‑то момент (перегрев?).
График для elbrus‑v6
Как видим, после 800КБ скорость начинает падать.
адепты АЙ опять не хотят учиться
Просто из научного интереса, а как SPIRAL справится с оптимизацией FFT для Эльбруса?..
У них есть библиотека FFTX с уже готовым символьным описанием алгоритма, оптимизирующая его в build-time, и есть даже онлайн-генератор кода для FPGA.
Из консольного вывода от разных размеров массива видно, что эталонный вариант всегда выполняется на 40–50 тактов быстрее. Я не смог понять, с чем это связано (особенности расположения кода в памяти?).
Это может быть связано с тем, что компилятор при выводе показывает один код, а при компиляции полной программы сгенерируект совсем другой, так как исполнение этого когда в контексте другого когда сильно меняет всю картину компилятору и он может например вообще выкинуть векторизацию, вставить какой-то рыхлый цикл с apb и он почему то оказывается быстрей чем все попытки векторизовать это все принудительно. Может конечно это и просто моя ошибка была, но в любом случае дамп бинарника никому ненавредит.
ldis ./stage_radix2_vector2_speedtest > stage_radix2_vector2_speedtest.ldis.txt
Положил «stage_radix2_vector2_speedtest.ldis.txt» и «stage_radix2_vector2_speedtest.s» рядом с соответствующим исходником (то есть сюда).
Если вы про мой код, то это было давно и я не помню где эти файлы. Но там проблема скорей всего была в том, что внешняя процедура с её циклом заняли и так много регистров и когда к ним добавлялись еще базированные-конвейризированные из раскрученного цикла векторизованной функции, регистровый файл видать закончился. Компилятор такие вещи при компиляции всей программы видит и недопускает, а если ему отдельные функции суют да еще с прагмами требующими раскрутки, ему понятное дело ничего не остается.
У вас же код обоих функций абсолютно идентичный, наверное просто не стоило все функции заталкивать в одну программу.
В предыдущем комментарии я выложил файл‑дамп бинарника, как было предложено.
В комментарий оно не влезает, даже под спойлером, хабр пишет, что так много текста нельзя.
Нас интересуют функции stage_radix2_vector2_etalon и stage_radix2_vector2.
Выделяю мышкой весь код функции и копирую в текстовый редактор.
Так же делаю со второй функцией.
Дальше Alt+Tab туда‑сюда (или Ctrl+Tab, если это вкладки редактора) — сравниваю код.
Извините, я не выспался. Я посмотрел файлы, сопоставил функции diff -ом посмотрел глазами main там и там, всё вроде одинаково. Единственное что может быть вот эта вот проверка как то влияет:
ldd,0 [ %dr12 + f64,lts0 0x17090 -> stage_radix2_vector2_algorithm ], %dr0
cmpedb,1 %dr0, f64,lts0 0x123a0 -> stage_radix2_vector2_etalon , %pred1
То есть он берет из вашей таблицы вызовов текущий адрес по индексу и проверяет зачем то сравнивает с адресом функции etalon дальше там идет предикатное наперстачничество и какой-то мухлеж это уже разобрать довольно сложно.
Может быть он решил заинлайнить функции из таблицы прямо в main?
В таком случае, я сравнивал не тот код.
Похоже еще веселее, вот три путя:
0155<000000013710> :
ldw,0 [ _f64,_lts1 0x170c0 -> stage_radix2_vector2_algorithm_count ], %b[16]
shld,1,sm %db[9], 0x4, %db[9]
subd,2,sm %dr2, _f16s,_lts0lo 0x10, %db[8]
addd,4,sm %dr2, %db[11], %dr45
fdivd,5,sm %db[14], %db[13], %dr41
disp %ctpr3, M_13be8
0156<000000013738> :
addd,0,sm %dr2, %db[11], %dr46
addd,1,sm %dr7, %db[9], %dr15
subd,2,sm %dr2, _f16s,_lts0lo 0x10, %dr60
addd,3,sm 0x0, _f64,_lts1 0x15ec0, %dr59
addd,4,sm %dr2, %db[11], %dr43
disp %ctpr2, M_138a0
0157<000000013760> :
faddd,3,sm %db[12], %db[12], %dr42
fdivd,5,sm %db[15], %db[13], %dr44
disp %ctpr1, M_13f30
0158<000000013770> :
cmplsb,0 0x0, %b[16], %pred1M_13be8 и M_138a0 практически одинаковые, и ведут примерно вот в такой код:
while(i < stage_radix2_vector2_algorithm_count) {
$clock_gettime@plt(...)
i++
$clock_gettime@plt(...)
print(...)
}обычный цикл с прыжком наверх как на интеле, и тут же обрабатываются элементы массива, вытаскиваются лоадом и вместо перестановки просто кладутся со смещением:
0308<000000013e28> :nop 3
fmuld,3 %dr9, %dr7, %dr5
fmuld,4 %dr36, %dr7, %dr6
0312<000000013e38> :
std,5 %dr5, [ %dr3 + _f16s,_lts0lo 0x30 ]
0313<000000013e48> :
std,5 %dr6, [ %dr3 + _f16s,_lts0lo 0x48 ]И третий бранч как раз ведет туда где косвенно вызываются функции как положено. Так что вот, компилятору на отдельные функции плевать он исходит из того что творится в процедуре, а тут цикличные вызовы с принтами и таймклоками в цикле, и смысла нет видимо заморачиваться.
Я думаю, M_13be8 и M_138a0 ведут к случаю, когда data_count слишком мал и цикл внутри функций stage_ не выполняется ни разу.
Поэтому там остаются только замеры времени.
А M_13f30 и M_14840 ведут уже к случаям нормального data_count.
Одна функция занимает примерно эти адреса (между чтениями регистра %clkr):
0352<000000013f30> M_13f30
...
0428<000000014540> M_14540
Вторая занимает примерно эти адреса (между чтениями регистра %clkr):
0511<000000014840> M_14840
...
0588<000000014e78> M_14e78
Разница в том, что у второй функции есть дыра на позиции 0566<000000014dd8>
Хорошая статья.
Попробуйте NEX CLI - сэкономит время на разборе расписания широких команд.
https://github.com/smworklair/NEX-CLI-x-Elbrus
Статья интересная, но меня еще волнует 1 вопрос: измеряли ли вы подготовку констант и планируете ли вы ускорять их подсчет? Ведь на практике перед выполнением БПФ мы неизбежно вызовем функцию, которая заготовит нам данные. И возникает вопрос, насколько тяжела по вычислениям подготовка констант, чтобы работать с БПФ в ЦОС: обычно это 1-2 раза "чистого" вычисления преобразования. Но с другой стороны для тысячи вызовов это уже будет не так важно
Не измерял.
Ускорение подсчёта происходит за счёт понижения точности.
Либо можно хранить готовые таблицы, которые тоже надо один раз считать откуда‑то.
Для ускорения подсчёта можно использовать те же методы, что при вычислении на ходу.
Не думал про этот вопрос.
У меня примерно такой был план развития (по мере сил):
— сначала сделать вариант с const‑коэффициентом на кусок данных и, возможно, исследовать границу между const‑коэффициентом и обычным способом получения коэффициентов;
— потом сделать вычисление коэффициентов на ходу;
— имея это всё, сделать тестилку, которая будет для разных размеров входных данных выбирать наилучший по скорости.
До этого надо перемежение коэффициентов везде вставить.
С этим неудобство в том, что для разных функций нужны разные способы перемежения.
Надо будет делать несколько массивов для каждого случая.
И потом либо разным функциям разные массивы коэффициентов давать (не совсем честное сравнение), либо копировать их в одно и то же место и уже указатель на это место давать в функцию.
могу предложить для предположений некоторые результаты DSP, о котором я упоминал в комментариях к прошлой статье автора. Библиотечная функция готовит коэффициенты для 1К,4К,16К БПФ за 21000,84000,336000 тактов. Абсолютно линейная зависимость. Алгоритм вычислений простой , может быть векторизован , легко разделен на ядра. Библиотечная функция тратит 5000,23500,120000 тактов на БПФ соответственно. С ростом размера БПФ уже можно думать о вычислении массива коэффициентов, а не табличной заготовке. Для Е2С3 вычисление коэффициентов будет намного быстрее ( в сравнении с DSP) и вес их, в сравнении со временем БПФ, будет поменьше.
Пишу алгоритм FFT на Си для процессора Эльбрус: прямая векторизация