Можноли как то ускорить код

[Missir]

mihaild в сообщении #1638068 писал(а):

А не получится поменять порядок - сначала посчитать 4 четверки, а потом в них всех искать максимумы одновременно (это за 3 вызова _mm_max_ps даст результат; т.к. там latency сильно больше, чем throughtput, а вычисления сильно зависимы, то лучше, если получится, считать сразу для 8 четверок)

получится, но данные придётся выгружать в массив, т.к. во время вычислений все регистры заняты под "завязку". Вообще не хотелось заводить для этого массив, да и выгружать все данные. Но всё равно спасибо, может попробую, хотя вариант GAA мне представляется более привлекательным.
Сейчас готовлю пробный вариант. По результатам отпишусь

[GAA]

Замечание к моему предыдущему сообщению. Использование инструкции для работы с целыми pshufd может приводить к штрафу при работе с single. Для работы с single предназначена shufps

Используется синтаксис ASM

  movaps xmm1, xmm0
  shufps xmm1, xmm1, 01001110b
  maxps  xmm0, xmm1
  movaps xmm1, xmm0
  shufps xmm1, xmm1, 10110001b
  maxps  xmm0, xmm1 

[Missir]

GAA
так то всё нормально, но почему то не работает тригонометрия
например не работает _mm_cos_ps(), можно как то их включить?

[worm2]

Код в стартовом сообщении вычисляет некоторую функцию от четырёх аргументов: $F(a, b, p, q)$ (переменные $p$ и $q$ в коде переиспользуются). Есть соблазн это всё ускорить, интерполируя не отдельные операции, типа квадратного корня, а всю функцию $F$ целиком, сначала вычислив её на какой-то не очень объёмной сетке. Если функция $F$ — промежуточная, а на самом деле вычислить нужно $G(F(...), ...)$, возможно, лучше целиком $G$ интерполировать.
Правда, методы многомерной интерполяции довольно сложны сами по себе, если они не линейные. А если линейные, то для достижения нужной точности понадобится чересчур много точек, к гадалке не ходи. К тому же нужно понимать, где в этом четырёхмерном объёме функция сильно меняется и там насыпать точек погуще.
В общем, это сложно и не факт, что получится. Но если получится, то выигрыш в скорости будет действительно радикальный. Может быть, размерность функции $G$ на самом деле меньше четырёх?

[Missir]

worm2
я это уже проверял, при сопоставимых вычислительных затратах точность получается ниже и диапазон ограничен,
в общем хуже чем так как сейчас

[GAA]

Missir в сообщении #1638132 писал(а):

например не работает _mm_cos_ps(), можно как то их включить?

Я этим вопросом не владею. В SSE нет тригонометрических функций. _mm_cos_ps() есть в Intel SVML (Short Vector Math Library).

[Missir]

Сделал без синусов, этот участок стал действительно работать раза в 3 быстрее, правда результаты пока выдаёт неверные
GAA
правильно ли я написал предложенную Вами перестановку с извлечением максимального элемента?

v4 =_mm_shuffle_ps(e4, e4, 0b01001110); e4 = _mm_max_ps(v4, e4);
v4 =_mm_shuffle_ps(e4, e4, 0b10110001); e4 = _mm_max_ss(v4, e4);
_mm_store_ss(&res, e4);

[Missir]

Перепроверил, исправил пару опечаток, но всё равно получается что то не то.
Не которые значения совпадают очень хорошо, примерно до 5 цифр, у других погрешность доходит до 10%.
Если бы ошибка, не должно быть таких точных результатов, тем более их много. Возможно так повлияли переход во float и изменение порядка вычислений

[Missir]

GAA пожалуйста объясните как работает _mm_shuffle_ps
этот код

v4 =_mm_shuffle_ps(e4, e4, 0b01001110); e4 = _mm_max_ps(v4, e4);
v4 =_mm_shuffle_ps(e4, e4, 0b10110001); e4 = _mm_max_ss(v4, e4);
_mm_store_ss(&res, e4);
if (res>max_nrm) max_nrm=res; 

изредка пропускает максимумы, а этот

_mm_store_ps(&p[0], e4);
if (p[0]>max_nrm) max_nrm=p[0]; 
if (p[1]>max_nrm) max_nrm=p[1]; 
if (p[2]>max_nrm) max_nrm=p[2]; 
if (p[3]>max_nrm) max_nrm=p[3];

работает стабильно, но так ведь не должно быть

[GAA]

Запустив gcc 13.2.0

Используется синтаксис C

#include <xmmintrin.h>
#include <stdio.h>
__m128 x, y;

float* p = (float*) &x;

int main(void)


{
p[0] = 0; p[1] = 1; p[2]=2; p[3] =3;
y = x;
y =_mm_shuffle_ps(y, y, 0b01001110);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
y = x;
y =_mm_shuffle_ps(y, y, 0b10110001);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
printf("----------\n");

p[0] = 3; p[1] = 0; p[2]=1; p[3] =2;
y = x;
y =_mm_shuffle_ps(y, y, 0b01001110);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
y = x;
y =_mm_shuffle_ps(y, y, 0b10110001);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
printf("----------\n");

p[0] = 2; p[1] = 3; p[2]=0; p[3] =1;
y = x;
y =_mm_shuffle_ps(y, y, 0b01001110);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
y = x;
y =_mm_shuffle_ps(y, y, 0b10110001);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
printf("----------\n");

p[0] = 1; p[1] = 2; p[2]=3; p[3] =0;
y = x;
y =_mm_shuffle_ps(y, y, 0b01001110);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
y = x;
y =_mm_shuffle_ps(y, y, 0b10110001);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
printf("----------\n");

return 0;
} 

Получим

2.000000, 3.000000, 2.000000, 3.000000
3.000000, 3.000000, 3.000000, 3.000000
----------
3.000000, 2.000000, 3.000000, 2.000000
3.000000, 3.000000, 3.000000, 3.000000
----------
2.000000, 3.000000, 2.000000, 3.000000
3.000000, 3.000000, 3.000000, 3.000000
----------
3.000000, 2.000000, 3.000000, 2.000000
3.000000, 3.000000, 3.000000, 3.000000
----------

-- Tue 07.05.2024 19:09:36 --

Missir, Ваш вариант также работает

Используется синтаксис C

#include <xmmintrin.h>
#include <stdio.h>
__m128 x, y;

float* p = (float*) &x;

int main(void)


{
p[0] = 1; p[1] = 2; p[2]=3; p[3] =0;
y =_mm_shuffle_ps(x, x, 0b01001110);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
y =_mm_shuffle_ps(x, x, 0b10110001);
x = _mm_max_ps(x, y);
printf("%f, %f, %f, %f\n", p[0], p[1], p[2], p[3]);
printf("----------\n");
return 0;
}

Результат:

1.000000, 2.000000, 3.000000, 0.000000
3.000000, 2.000000, 3.000000, 2.000000
3.000000, 3.000000, 3.000000, 3.000000
----------

Поэтому хотелось бы увидеть для каких данных работает неправильно.

-- Tue 07.05.2024 19:18:53 --

В целом, использование intrinsic создаёт иллюзию написания на языке высокого уровня, но, лично мне это не нравится. Работая непосредственно с инструкциями и регистрами, можно минимизировать, например пересылки данных или сохранение в памяти (переменных/стеке) промежуточных результатов. Я стараюсь так оформлять последовательность инструкций, чтобы по завершению вычислений результат (результаты) уже лежали в тех регистрах (регистре), в которых (в котором) функция должна возвращать (в соответствии с принятыми соглашениями в используемом компиляторе).

[Missir]

GAA в сообщении #1638399 писал(а):

хотелось бы увидеть для каких данных работает неправильно

боюсь, что "выудить" их от туда будет непросто, хотя, попытаться можно
спасибо Вам за примеры, сначала попробую разобраться в них, может и так пойму в чём дело

кстати, я обещал отписаться, код на sse сейчас работает корректно, ускорение на оптимизируемом фрагменте почти в 4 раза,
правда точность резко упала, в основном E-5, но бывает и E-3. В double в основном было E-15, и в самых худших случаях E-10.
Но вообще, и такая точность вполне подходит.

[GAA]

Мне в ЛС написали, что в случае Single можно короче.
Совсем не обязательно во все четырёх single получать (одинаковый) результат.
В следующем коде, результат только в младшем single. Однако перенос а половину регистра :(.

Используется синтаксис ASM

movaps  xmm1, xmm0           // xmm0 = (A B C D)
shufps  xmm1, xmm1, 00011011b// L/T: 1/1            : xmm1 = (D C D A)
maxps   xmm0, xmm1           // L/T: 3/1 или (4/0.5): xmm0 = (max(A, D), max(C, B), max(C, D), max(A, D))
movhlps xmm1, xmm0           // L/T: 1/1            : xmm1 = (    x,        x,      max(A, D), max(C, B))
maxss   xmm0, xmm1           // L/T: 3/1            : xmm0 = (     x,       x,       x,  max(max(A, D), max(C, B)))

[Missir]

GAA
правильно ли я понимаю, что единицы в маске обозначают элементы, которые будут выбраны из 8 элементов в 4-х элементный результат? самая левая единица отмечает самый левый элемент результата, и т.д.?

но тогда в xmm1 должно быть C A C D, а у Вас D C D A

-- 08.05.2024, 00:15 --

У меня были такие маски:
01001110
10110001
первая перестановка даёт:
B A B C
они сравниваются с:
A B C D
вторая перестановка даёт:
B B C C
A C D D
в итоге, в младшем элементе будет максимум из последнего столбца, т.е. из C D C D
правильно я всё понял?

-- 08.05.2024, 00:20 --

GAA в сообщении #1638431 писал(а):

Совсем не обязательно во все четырёх single получать (одинаковый) результат

ну да, сначала сравнение блоков, а потом одиночное, в моём коде выше, к стати, так и сделано, там второй _mm_max_ss

[GAA]

GAA в сообщении #1638431 писал(а):

shufps xmm1, xmm1, 00011011b// L/T: 1/1 : xmm1 = (D C D A)

Опечатка. Должно быть: xmm1 = (D C B A) .

Если источник и приёмник один и тот же регистр, то нулевая двойка в маске указывает куда будет перенесён нулевой Single, первая двойка - куда будет перенесён первый Single, вторая - второй Single и третья - третий Single.
00 - в позицию нулевого Single;
01 - в позицию первого Single;
10 - в позицию второго Single;
11 - в позицию третьего Single.

A B C D 00011011 -> D C B A

Missir в сообщении #1638435 писал(а):

У меня были такие маски: 01001110 <..>

A B C D 01001110 -> C D A B.

----------------

Используется синтаксис C

#include <xmmintrin.h>
#include <stdio.h>
__m128 x, y;

float* p = (float*) &x;

int main(void)
{
p[0] = 3; p[1] = 2; p[2]=1; p[3] =0;
printf("%f, %f, %f, %f\n", p[3], p[2], p[1], p[0]);
x =_mm_shuffle_ps(x, x, 0b01001110);
printf("%f, %f, %f, %f\n", p[3], p[2], p[1], p[0]);
return 0;
} 

0.000000, 1.000000, 2.000000, 3.000000
2.000000, 3.000000, 0.000000, 1.000000

[GAA]

Пытался упростить и сильно заврался, спасибо за указание мне ошибки в ЛС.
Значение двойки не куда, а откуда, а номер двойки - куда.
A B C D 00011011 -> D C B A
с права налево:
11 из третьего Single в нулевой Single;
10 из второго Single в первый Single;
01 из первого Single в 0 второй Single;
00 из нулевого Single в третий Single.

A B C D 01001110 -> C D A B
с права налево:
10 из второго Single в нулевой: xBxx -> xxxB;
11 из третьего Single в первый: Axxx -> xxAx;
00 из нулевого Single во второй: xxxD -> xDxx;
01 из первого Single в третий: xxСx -> Cxxx;
(Спать давно надо)