Сделал FIR-фильтрацию с использованием SSE-инструкций

dimon72

Написал-таки. Может кому пригодится.
  
type
TVectSingle = array of single;
////////////////////////////////////////////////////////////////////////////////
function FIRFilSSE(const X: array of single; const C: array of single): TVectSingle;
//функция FIR-фильтра с использованием SSE-команд
//Обязат.условие: кол-во компонентов в векторе С коэф-тов свёртки должно быть кратным 4-м.
var
S,SC,N: integer;
Vector: TVectSingle;
begin
SC:=SizeOf(C);//определим длину вектора коэф-тов (C) в байтах; необх.для опред-я кол-ва смещений адреса по 16
S:=SizeOf(X)+4-SC;//определим кол-во "шагов" свёртки в байтах по вх.вектору
N:=High(C);//
Finalize(Vector);
SetLength(VectorHigh(X)+1-N;//установить длину вых.вектора (Y)
asm
mov ebx,0//инициализ-я переменной цикла @lb2 (задание инкремента смещения адреса по 16)
mov esi,0//инициализ-я переменной цикла @lb1 (задание инкремента смещения адреса по 4)
mov edx, X//в edx-адрес первого эл-та вх.вектора (X)
mov ecx, C //в ecx-адрес первого эл-та вектора коэф-тов (C)
mov eax, Vector//здесь будет вых.вектор (результат фильтрации)
@lb1:
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xorps xmm2, xmm2//здесь результат свертки на X-N шаге; инициализируем нулем
push eax
@lb2:
//****************************************************************************
//здесь расположим команды asm цикла (аналог for i:=.. to .. do ...)
// загрузить вх. векторы в регистры xmm0-1 по соотв.адресам со смещением по 16 байт
mov eax,esi
add eax,ebx
movups xmm0,[edx+eax]//загрузить вх.вектор X в xmm0
movups xmm1,[ecx+ebx]//загрузить вх.вектор C в xmm0
mulps xmm1, xmm0//умножить 4 к-та вектора коэф-тов С на 4 к-та вх. вектор X и рез-т в xmm0
addps xmm2, xmm1// сложить с очередным результатом 4-х компонентной суммации

add ebx,16//задаём очередной "шаг" = 16 байт
cmp ebx,SC//...to [кол-во "шагов" по 16 байт, кот. нужно сделать] do
jnz @lb2
//****************************************************************************
mov ebx,0//контрольный объект цикла @lb2 - в исх. состояние
//посчитаем сумму эл-тов, получ-х в рез-те данного шага свёртки
movups xmm3, xmm2 //Копировать результаты. Порядок в xmm3: f3,f2,f1,f0
shufps xmm3, xmm3, 4Eh //Переставить: f1,f0,f3,f2
addps xmm2, xmm3 //Сложить: f3+f1,f2+f0,f1+f3,f0+f2
movups xmm3, xmm2 //Копировать результаты
shufps xmm3, xmm3, 11h //Переставить: f0+f2,f1+f3,f0+f2,f1+f3
addps xmm2, xmm3 //Сложить: f0+f1+f2+f3,f0+f1+f2+f3,f0+f1+f2+f3,f0+f1+f2+f3
// результат во всех элементах xmm2
//
pop eax
// записать результат в Vector (младшее число из xmm2)
movss [eax+esi], xmm2
//
add esi,4
cmp esi,S
jnz @lb1
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
mov esi,0//контрольный объект цикла @lb1 - в исх. состояние
end;
Result := Vector;
end;
////////////////////////////////////////////////////////////////////////////////

dimon72

Вот вариант функции без ASM'а:
 
////////////////////////////////////////////////////////////////////////////////
function FIRFil(const X: array of single; const C: array of single): TVectSingle;
//функция FIR-фильтра
var
i,z: integer;
N: integer;
begin
// z -номер отсчета, N - длина фильтра, X - вход, Y - выход , С - коэффициенты.
N:=High(C);
Finalize(Result);
SetLength(ResultHigh(X)+1-N;
For z:=0 to High(X)-N do
begin
For i:=0 to N do
Result[z]:=Result[z]+X[z+i]*C[i];
end;
end;
////////////////////////////////////////////////////////////////////////////////

Тестировал обе реализации функций на двумерном массиве 27х292000. Кол-во коэф. прогоняемого фильтра – 528.
Методика тестирования: запускал фильтрацию всего массива одновременно с внешним секундомером (реализованным не в компьютере). Делал 5 замеров, которые усреднял и округлял до целого.
Результаты:
Время работы функции FIRFilSSE: 4 сек.
Время работы функции FIRFil: 32 сек.

SPARTAK3959

А регистры ebx и esi в ассемблерной вставке дядя сохранять будет?

evgen5555

А нах?

Ivan8209

> //Обязат.условие: кол-во компонентов в векторе С коэф-тов свёртки должно быть кратным 4-м.
Религия не позволяет использовать assert или просто отсутствие оного?
---
"Аллах не ведёт людей неверных."

Olenenok

лучше бы ты распараллелил бы прогу на нормальном языке, имхо, было бы больше смысла.

dimon72

Как сумел, так сделал. Критиковать всегда легче.

dimon72

лучше бы ты распараллелил бы прогу на нормальном языке, имхо, было бы больше смысла.
Во-первых, что такое "нормальный язык";
во-вторых, не понимаю, как это "распараллелил"?
Ты сможешь, скажем, без ассемблера, на "нормальном языке" сделать так, чтобы в программе 4 операции выполнялись за один такт процессора? Просто я вообще не программист и даже не технарь по образованию. Мне интересно. Приведи пример.
И, наконец, что значит "было бы больше смысла"? Мне нужно было решить задачу - я это сделал. И, на мой взгляд, получилось удовлетворительно.

Dasar

> во-вторых, не понимаю, как это "распараллелил"?
в современных компьютерах, кол-во процессоров больше, чем 1.

Marinavo_0507

Больше чем 1 - это, как правило, всего 2, причём один из них более чем занят вирусами и антивирусами

Dasar

> причём один из них более чем занят вирусами и антивирусами
вирусы/антивирусы скорее жрут винт, чем проц

Dasar

размер матрицы кстати какой?

dimon72

в современных компьютерах, кол-во процессоров больше, чем 1.
Да я понимаю. Суперскалярность там, дуал-коры и пр. Здесь я хотел уточнить у товарища, что он понимает под "параллелить на нормальном языке". SSE и пр. расширения ведь не просто так придумали, а для ускорения вычислений. Почему бы их не использовать?

dimon72

размер матрицы кстати какой?
Прогонял на массиве 27х292000. Кол-во коэф-тов прогоняемого фильтра - 528.

Dasar

> Почему бы их не использовать?
потому что прога перестала читаться, соответственно сложность дальнейших модификаций резко подскочила.

dimon72

потому что прога перестала читаться, соответственно сложность дальнейших модификаций резко подскочила.
Не знаю, для меня она прекрасно читается. К тому же, это - один из первых моих опытов использования ASM'а. А для тех, кто его знает, думаю, тем более будет нормально восприниматься.
Хорошо, какие есть варианты? Чтобы также быстро работало, без многопроцессорности и пр. фич?

Dasar

ответ уже звучал: взять готовую либу, скорость скорее всего еще в несколько раз подскочит.

dimon72

И так пока что устраивает.

Codcod

Вопрос:
asm
mov ebx,0//инициализ-я переменной цикла @lb2 (задание инкремента смещения адреса по 16)
mov esi,0//инициализ-я переменной цикла @lb1 (задание инкремента смещения адреса по 4)
mov edx, X//в edx-адрес первого эл-та вх.вектора (X)
mov ecx, C //в ecx-адрес первого эл-та вектора коэф-тов (C)
mov eax, Vector//здесь будет вых.вектор (результат суммы векторов)
@lb1:
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xorps xmm2, xmm2//здесь результат свертки на X-N шаге; инициализируем нулем
push eax
@lb2:
//****************************************************************************
//здесь расположим команды asm цикла (аналог for i:=.. to .. do ...)
// загрузить вх. векторы в регистры xmm0-1 по соотв.адресам со смещением по 16 байт
mov eax,esi
add eax,ebx
movups xmm0,[edx+eax]//загрузить вх.вектор X в xmm0
movups xmm1,[ecx+ebx]//загрузить вх.вектор C в xmm0
mulps xmm1, xmm0//умножить 4 к-та вектора коэф-тов С на 4 к-та вх. вектор X и рез-т в xmm0
addps xmm2, xmm1// сложить с очередным результатом 4-х компонентной суммации
add ebx,16//задаём очередной "шаг" = 16 байт
cmp ebx,SC//...to [кол-во "шагов" по 16 байт, кот. нужно сделать] do
jnz @lb2
Зачем весь этот огород с EAX? что бы потерять 1 такт на цикле?

dimon72

Зачем весь этот огород с EAX? что бы потерять 1 такт на цикле?
Согласен. Мне самому не нравится.
Но прямо: movups xmm0,[edx+ebx+esi], не получается. А как по другому можно?

tamusyav

Можно менять edx, а не esi. Правда, выигрыш от этого будет мизерный (проверял).

Codcod

на вскидку (нога прошу не пинать, проверять времени не было)
asm
  mov ebx,S
  mov esi,X
  mov edi,Vector
  mov edx,C
@lb1:
  xorps xmm2,xmm2
  mov eax,SC
@lb2:
  movups xmm0,[edx+eax] //load Coeff
  movups xmm1,[esi+eax] //load Data
  mulps xmm1,xmm0
  addps xmm2,xmm1
  sub eax,16
  jnz $lb2
  movups xmm3, xmm2 //Копировать результаты. Порядок в xmm3: f3,f2,f1,f0
  shufps xmm3, xmm3, 4Eh //Переставить: f1,f0,f3,f2
  addps xmm2, xmm3 //Сложить: f3+f1,f2+f0,f1+f3,f0+f2
  movups xmm3, xmm2 //Копировать результаты
  shufps xmm3, xmm3, 11h //Переставить: f0+f2,f1+f3,f0+f2,f1+f3
  addps xmm2, xmm3 //Сложить: f0+f1+f2+f3,f0+f1+f2+f3,f0+f1+f2+f3,f0+f1+f2+f3
  movss [edi], xmm2
  
  add edi,4
  add esi,4
  sub ebx,4
  jnz @lb1
end;
 что то типа этого, ещё можно поиграться с jmp, потому как ставить jnz в конце цикла не разумно, по крайней мере для P4.
Да ещё рекомендуется сохрянять регистры

dimon72

нога прошу не пинать, проверять времени не было
Что мне лягаться-то? Я вообще дилетант в Asm'е.
Пробовал предложенный тобой код, но, к сожалению, он работает неправильно.
Идею здесь я, кажется, понял: адресоваться при вычислении на каждом шаге свертки с конца к началу. Дело в том, что команда movups загружает значения из массива от указанного в квадратных скобках адреса и далее, по возрастанию, а не наоборот.

dimon72

Да ещё рекомендуется сохрянять регистры
Ты имеешь в виду это?

asm
pushad
.....
popad
end;

dimon72

Попробовал изменить вот так:

asm
pushad//размещение eax, ebx, ecx, edx, esi, edi, esp, ebp в стеке
mov ebx,S//переменная цикла @lb1 (опред-е кол-ва шагов свертки деремир-ем от S по 4)

mov esi,X//в esi-адрес первого эл-та вх.вектора (X)
mov edx,C//в edx-адрес первого эл-та вектора коэф-тов (C)
mov edi,Vector//здесь будет вых.вектор (результат фильтрации)
@lb1:
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xorps xmm2,xmm2//здесь результат свертки на X-N шаге; инициализируем нулем
mov eax,0//инициализ-я переменной цикла @lb2 (задание инкремента смещения адреса по 16)
@lb2:
//****************************************************************************
movups xmm1,[esi+eax] //загрузить вх.вектор (X)
movups xmm0,[edx+eax] //загрузить вектор коэфф. (С)
mulps xmm1,xmm0
addps xmm2,xmm1

add eax,16//задаём очередной "шаг" = 16 байт
cmp eax,SC//...to [кол-во "шагов" по 16 байт, кот. нужно сделать] do
jnz @lb2
//****************************************************************************
//посчитаем сумму эл-тов, получ-х в рез-те данного шага свёртки
movups xmm3, xmm2 //Копировать результаты. Порядок в xmm3: f3,f2,f1,f0
shufps xmm3, xmm3, 4Eh //Переставить: f1,f0,f3,f2
addps xmm2, xmm3 //Сложить: f3+f1,f2+f0,f1+f3,f0+f2
movups xmm3, xmm2 //Копировать результаты
shufps xmm3, xmm3, 11h //Переставить: f0+f2,f1+f3,f0+f2,f1+f3
addps xmm2, xmm3 //Сложить: f0+f1+f2+f3,f0+f1+f2+f3,f0+f1+f2+f3,f0+f1+f2+f3
// результат во всех элементах xmm2
movss [edi], xmm2// записать результат в Vector (младшее число из xmm2)

add edi,4//...следующее число вых.вектора (Vector)
add esi,4//...следующее число вх.вектора (X)
sub ebx,4//...декремент перем.цикла @lb1 (кол-во шагов свёртки)
jnz @lb1
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
popad//извлечение eax, ebx, ecx, edx, esi, edi, esp, ebp из стека
end;

Работает, и команд меньше, чем было.
Спасибо!
что то типа этого, ещё можно поиграться с jmp, потому как ставить jnz в конце цикла не разумно, по крайней мере для P4.
Да ещё рекомендуется сохрянять регистры

jmp, насколько я понимаю, безусловный переход. Но мне же проверять как-то надо, сколько раз внутренний цикл по @lb2 прогонять.

tamusyav

В обратную сторону бегать по массиву - это долго.

Viktori68

Такой ещё код посчитай a.c

#define N 1000
#define M 100

double a[M], b[M], c[N];

int main
{
int i, j;

for (i = 0; i < M - N; i++)
for (j = 0; j < M; j++)
b[i] = b[j] + a [j + i] * c[i];

return 1;
}

gcc -msse2 -mfpmath=sse
 
.file "a.c"
.text
.globl main
.type main, @function
main:
leal 4(%esp %ecx
andl $-16, %esp
pushl -4(%ecx)
pushl %ebp
movl %esp, %ebp
pushl %ecx
subl $16, %esp
movl $0, -12(%ebp)
jmp .L2
.L3:
movl $0, -8(%ebp)
jmp .L4
.L5:
movl -12(%ebp %edx
movl -8(%ebp %eax
movsd b%eax,8 %xmm2
movl -12(%ebp %eax
addl -8(%ebp %eax
movsd a%eax,8 %xmm1
movl -12(%ebp %eax
movsd c%eax,8 %xmm0
mulsd %xmm1, %xmm0
addsd %xmm2, %xmm0
movsd %xmm0, b%edx,8)
addl $1, -8(%ebp)
.L4:
cmpl $99, -8(%ebp)
jle .L5
addl $1, -12(%ebp)
.L2:
cmpl $-900, -12(%ebp)
jl .L3
movl $1, %eax
addl $16, %esp
popl %ecx
popl %ebp
leal -4(%ecx %esp
ret
.size main, .-main
.comm a,800,32
.comm b,800,32
.comm c,8000,32
.ident "GCC: (GNU) 4.1.1 20061011 (Red Hat 4.1.1-30)"
.section .note.GNU-stack,"",@progbits

dimon72

Задался я тут недавно вопросом, как выравнивать данные по 16 байт в Дельфи?
Нужно для того, что SSE-команды для выровненных данных работают быстрее невыровненных.
Пример команд: movaps вместо movups и т.д.
Средства среды Delphi не позволяют осуществлять такое выравнивание, т.е., например,
в Project->Options->Compiler->Record field alignment есть следующие значения:
1,2,4,8, а 16 нету. Директива компилятора {$ALIGN 16} тоже не прокатывает.
посоветовал выделять память и выравнивать вручную:
выделяем выравненные куски + прячем от остального кода работу по выравниванию (решение в лоб)
code:
byte* AllocAlign16(int size)
{
byte * p = new byte[size + 17];
byte *p16 = p + 16 - int)p)%16; //выравниваем
(p16-1)^ = p16-p; //запоминаем сдвиг, чтобы потом можно было восстановить истинный указатель на блок
return p16;
}
void FreeAlign16(byte *p)
{
byte * oldP = p - (p-1)^;
delete oldP;
}
Поскольку я в С не шарю (да и в Дельфи так, не особо то мало что понял из предложенного кода и попробовал
по-своему. Допустим, есть 2 вектора, которые нужно занести в память, начиная с адресов, кратных 16.
Вот что у меня получилось:
 

Type
TArr = array [0..8] of single;
PArr = ^TArr;

var
V0: PArr;
V1: PArr;

implementation

procedure TForm1.FormCreate(Sender: TObject);
var
Address: integer;
begin
/////////////////////////////////////////////////////////////////////////////////
//выравнивание для первого массива
V0:=nil;
V0:=AllocMem(36+16);//выделим память размером большим (+16 байт чем надо
//(для учёта возможного инкремир-я указ-ля при выравнивании)
Address:=Integer(V0);//получим знач-е возвращенного ф-цией new адреса

While (Address mod $10)<>0 do//пока адрес не кратен 16 (остаток от деления<>0)
Address:=Address+4;//инкремируем его по 4 байта...

V0:= Ptr(Address);//...как только кратен 16 - присваиваем его указателю.
/////////////////////////////////////////////////////////////////////////////////
//выравнивание для второго массива
V1:=nil;
V1:=AllocMem(36+16);//выделим память размером чуть большим (+16 байт чем надо
//(для учёта возможного инкремир-я указ-ля при выравнивании)
Address:=Integer(V1);//получим знач-е возвращенного ф-цией new адреса

While (Address mod $10)<>0 do//пока адрес не кратен 16 (остаток от деления<>0)
Address:=Address+4;//инкремируем его по 4 байта...

V1:= Ptr(Address);//...как только кратен 16 - присваиваем его указателю.
/////////////////////////////////////////////////////////////////////////////////
//заполняю входные массивы по скорректированным адресам
V0^[0]:=6;
V0^[1]:=3;
V0^[2]:=7;
V0^[3]:=3;
//
V0^[4]:=1;
V0^[5]:=2;
V0^[6]:=3;
V0^[7]:=4;
V0^[8]:=5;


V1^[0]:=4;
V1^[1]:=5;
V1^[2]:=8;
V1^[3]:=7;
//
V1^[4]:=8;
V1^[5]:=3;
V1^[6]:=9;
V1^[7]:=1;
V1^[8]:=4;


end;

Потом, если дальше работать с V0 и V1, то всё нормально.
Проблема в следующем: если вместо явных указателей V0 и V1 на статические массивы сделать их, скажем, типом array of single, то программа падает.
Что тут не правильно? Динамические же массивы, как я понимаю, и так являются указателями.

SPARTAK3959

Динамические массивы хранят размер массива и количество выделенной памяти по отрицательным смещениям.

dimon72

Во, блин! Ну и как мне переделать, чтоб работало?
Оставить комментарий
Имя или ник:
Комментарий: