Сделал FIR-фильтрацию с использованием SSE-инструкций
Вот вариант функции без ASM'а:
Тестировал обе реализации функций на двумерном массиве 27х292000. Кол-во коэф. прогоняемого фильтра – 528.
Методика тестирования: запускал фильтрацию всего массива одновременно с внешним секундомером (реализованным не в компьютере). Делал 5 замеров, которые усреднял и округлял до целого.
Результаты:
Время работы функции FIRFilSSE: 4 сек.
Время работы функции FIRFil: 32 сек.
////////////////////////////////////////////////////////////////////////////////
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 сек.
А регистры ebx и esi в ассемблерной вставке дядя сохранять будет?
А нах?
> //Обязат.условие: кол-во компонентов в векторе С коэф-тов свёртки должно быть кратным 4-м.
Религия не позволяет использовать assert или просто отсутствие оного?
---
"Аллах не ведёт людей неверных."
Религия не позволяет использовать assert или просто отсутствие оного?
---
"Аллах не ведёт людей неверных."
лучше бы ты распараллелил бы прогу на нормальном языке, имхо, было бы больше смысла.
Как сумел, так сделал. Критиковать всегда легче.
лучше бы ты распараллелил бы прогу на нормальном языке, имхо, было бы больше смысла.Во-первых, что такое "нормальный язык";
во-вторых, не понимаю, как это "распараллелил"?
Ты сможешь, скажем, без ассемблера, на "нормальном языке" сделать так, чтобы в программе 4 операции выполнялись за один такт процессора? Просто я вообще не программист и даже не технарь по образованию. Мне интересно. Приведи пример.
И, наконец, что значит "было бы больше смысла"? Мне нужно было решить задачу - я это сделал. И, на мой взгляд, получилось удовлетворительно.
> во-вторых, не понимаю, как это "распараллелил"?
в современных компьютерах, кол-во процессоров больше, чем 1.
в современных компьютерах, кол-во процессоров больше, чем 1.
Больше чем 1 - это, как правило, всего 2, причём один из них более чем занят вирусами и антивирусами 
> причём один из них более чем занят вирусами и антивирусами
вирусы/антивирусы скорее жрут винт, чем проц
вирусы/антивирусы скорее жрут винт, чем проц
размер матрицы кстати какой?
в современных компьютерах, кол-во процессоров больше, чем 1.Да я понимаю. Суперскалярность там, дуал-коры и пр. Здесь я хотел уточнить у товарища, что он понимает под "параллелить на нормальном языке". SSE и пр. расширения ведь не просто так придумали, а для ускорения вычислений. Почему бы их не использовать?
размер матрицы кстати какой?Прогонял на массиве 27х292000. Кол-во коэф-тов прогоняемого фильтра - 528.
> Почему бы их не использовать?
потому что прога перестала читаться, соответственно сложность дальнейших модификаций резко подскочила.
потому что прога перестала читаться, соответственно сложность дальнейших модификаций резко подскочила.
потому что прога перестала читаться, соответственно сложность дальнейших модификаций резко подскочила.Не знаю, для меня она прекрасно читается. К тому же, это - один из первых моих опытов использования ASM'а. А для тех, кто его знает, думаю, тем более будет нормально восприниматься.
Хорошо, какие есть варианты? Чтобы также быстро работало, без многопроцессорности и пр. фич?
ответ уже звучал: взять готовую либу, скорость скорее всего еще в несколько раз подскочит.
И так пока что устраивает.
Вопрос:
asmЗачем весь этот огород с EAX? что бы потерять 1 такт на цикле?
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 такт на цикле?Согласен. Мне самому не нравится.
Но прямо: movups xmm0,[edx+ebx+esi], не получается. А как по другому можно?
Можно менять edx, а не esi. Правда, выигрыш от этого будет мизерный (проверял).
на вскидку (нога прошу не пинать, проверять времени не было)
Да ещё рекомендуется сохрянять регистры
asmчто то типа этого, ещё можно поиграться с jmp, потому как ставить jnz в конце цикла не разумно, по крайней мере для P4.
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;
Да ещё рекомендуется сохрянять регистры
нога прошу не пинать, проверять времени не былоЧто мне лягаться-то?
Я вообще дилетант в Asm'е.Пробовал предложенный тобой код, но, к сожалению, он работает неправильно.
Идею здесь я, кажется, понял: адресоваться при вычислении на каждом шаге свертки с конца к началу. Дело в том, что команда movups загружает значения из массива от указанного в квадратных скобках адреса и далее, по возрастанию, а не наоборот.
Да ещё рекомендуется сохрянять регистрыТы имеешь в виду это?
asm
pushad
.....
popad
end;
Попробовал изменить вот так:
Работает, и команд меньше, чем было.
Спасибо!
jmp, насколько я понимаю, безусловный переход. Но мне же проверять как-то надо, сколько раз внутренний цикл по @lb2 прогонять.
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 прогонять.
В обратную сторону бегать по массиву - это долго.
Такой ещё код посчитай a.c
gcc -msse2 -mfpmath=sse
#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
Задался я тут недавно вопросом, как выравнивать данные по 16 байт в Дельфи?
Нужно для того, что SSE-команды для выровненных данных работают быстрее невыровненных.
Пример команд: movaps вместо movups и т.д.
Средства среды Delphi не позволяют осуществлять такое выравнивание, т.е., например,
в Project->Options->Compiler->Record field alignment есть следующие значения:
1,2,4,8, а 16 нету. Директива компилятора {$ALIGN 16} тоже не прокатывает.
посоветовал выделять память и выравнивать вручную:
по-своему. Допустим, есть 2 вектора, которые нужно занести в память, начиная с адресов, кратных 16.
Вот что у меня получилось:
Потом, если дальше работать с V0 и V1, то всё нормально.
Проблема в следующем: если вместо явных указателей V0 и V1 на статические массивы сделать их, скажем, типом array of single, то программа падает.
Что тут не правильно? Динамические же массивы, как я понимаю, и так являются указателями.
Нужно для того, что 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, то программа падает.
Что тут не правильно? Динамические же массивы, как я понимаю, и так являются указателями.
Динамические массивы хранят размер массива и количество выделенной памяти по отрицательным смещениям.
Во, блин! 
Ну и как мне переделать, чтоб работало?
Ну и как мне переделать, чтоб работало?Оставить комментарий
dimon72
Написал-таки. Может кому пригодится.