О степенной функции, OpenMP и прочих граблях

Есть вот такое вот известное выражение, переложеное из википедии:

#define To_sRGB(q) (((q)<=0.0031308f)? 12.92f*(q) : (1+0.055f)*powf(q,1/2.4f)-0.055f)
В том смысле, что у sRGB внизу - линейный участок, а дальше степенная функция.

Если его напрямую запрограммировать (вот прямо с powf()), то в зависимости от компилятора получается от 7 Mpix/sec (Visual C++ 2010, один поток) до 100Mpix/sec (Intel Parallel Studio XE, включен OpenMP). Пиксель - это 3 компонента по 4 байта (float), то бишь гигабайт с небольшим в секунду - это максимум.

Как-то медленно.

В-принципе, бывает быстрое приближенное возведение в степень для ограниченного диапазона значений. Более того, можно просто приблизить sRGB-кривую полиномом не очень высокой степени и получить вполне приемлемый для практики результат. Но хотелось заодно решить задачу быстрого наложения произвольной кривой.

Табличный лукап

Если просто домножить входное значение (диапазон которого 0..1) на размер таблицы, округлить и взять значение из таблицы, то получается недостаточно точно, не для того мы заводили плавающую точку, чтобы сразу округлять. Скорость для этого решения я мерять не стал.

Линейная интерполяция

Можно брать из таблицы два значения и линейно интеполировать, прямо как велит нам Klaus Post. Для таблицы из 512 записей получается быстро, но неточно (200-210 Mpix/sec при максимальной ошибке около 0.025), а для таблицы в 64k записей - медленнее и тоже не очень точно (150Mpix/sec, максимальная ошибка 3-4x10-4). Растить таблицу больше 64к записей (512kb памяти с учетом удвоения ради выравнивания, как велит нам Клаус) не хочется совсем.

Кубические сплайны

Кубические сплайны оказались отличным компромиссом: для таблицы в 4k записей (64 килобайта памяти) получаем 80 мегапикселей в секунду для реалиазации на C++, 120 Mpix/sec для SSE-варианта и ошибку не более чем 3.5x10-6. Увеличение размера таблицы не нужно, уменьшение тоже на скорость уже почти не влияет, а на точность - влияет очень сильно.

OpenMP

Все цифры выше (кроме одной, в самом начале, где это оговорено) - это в один поток. OpenMP любит нас, мы ее, режем входной поток на куски по несколько десятков мегабайт, #pragma omp parallel for, запускаем, получаем 360-400Mpix/sec для SSE-реализации сплайнов на компиляторе Intel и четырех ядрах.

Правда вылезают интересные грабли, ради которых, собственно, пост и написан (для собственной памяти):

  1. Режим округления у SSE сбивается при инициализации потока. То есть _MM_SET_ROUNDING_MODE надо звать уже внутри OpenMP-блока. Убил бы нахрен!
  2. У Visual Studio 2010 с OpenMP все не так радужно. OpenMP есть, считает правильно, только вот OpenMP-версия не в 3-4 раза быстрее однопоточной, а в полтора-два раза медленнее.

Заметим, что 4-4.8 гигабайта в секунду - это быстрее, чем залить картинку в видеокарту, посчитать (допустим, мгновенно) и вылить оттуда результат. Я к тому, что на внешнем сопроцессоре мелкими кусками обрабатывать смысла нет, туда надо нести сразу много кода.

Да, а "в 4 раза быстрее простой реализации" - очень похоже на константу, не первый раз ее уже вижу.

Comments

если нужно считать степень для многих пикселей (вектор или матрица), советую попробовать Vector Mathematical Functions из Intel MKL. там правда нет условия, как в примере. т.ч. нужно будет сначала считать правую часть для временного вектора, а потом делать условие и возвращать что нужно.
возможно подобная функция есть готовая в Intel IPP, т.к. там есть примитивы для фото и видео.

В IPP, насколько я вижу в мануале, гамма по BT.709, а порулить коэффициентами не дают. И, главное, представление изображений - целое, на что я пойтить не могу.

А вся конструкция хороша тем, что все поканальные преобразования запихиваются в одну кривую (из, скажем, 4k точек) и все вместе, вне зависимости от сложности формулы, считается за линейное и очень приемлемое время.
Т.е. и гамму и какую-то кривую и еще что-нибудь (blacks) можно наложить одним действием.

А вся конструкция хороша тем, что все поканальные преобразования запихиваются в одну кривую (из, скажем, 4k точек) и все вместе, вне зависимости от сложности формулы, считается за линейное и очень приемлемое время.

Хе-хе. Я давно ту же мысль думаю -- если много слоёв преобразований компилировать всё в одну кривую, если нет смешивающих каналы преобразований. Если есть -- то становится сложнее, конечно. Потому как делать три матрицы 64kx64kx64k в случае 16-битного цвета уже как-то черезчур.

Даже одна матрица размером 248 - это уже перебор. И даже размером 236 (3-мерный сплайн по 12 бит на сторону) - тоже перебор. Т.е. смешение каналов хорошо бы сводить к линейной операции, размером 3x3, возможно с каким-то поканальным препроцессингом до и после.

У меня есть ощущение, что все (разумные) смешивающие преобразования в RGB будут несмешивающими в Lab. Нет, ясно, что есть такая штука как Channel mixer делающая любой канал линейной комбинацией других (возможно даже в другом цветовом пространстве!) но часто ли это используют люди (кроме хардкорных последователей старого Маргулиса)? Особенно в фото а не в полиграфии?
Я, в общем, думаю свои мысли применительно к уровню выше RAW -- пост-обработка, заменф фотошопу (толстый)/лайтруму (мне там не хватает некоторых фотошоповых штук). Хотя, конечно, интеграция с RAW должна быть плотнейшая.

Что-то давно я такого объема спама в ЖЖ не встречал, просто за гранью добра и зла уже.

О да.
Все воют, теперь и до меня очередь дошла

> OpenMP есть, считает правильно, только вот OpenMP-версия не в 3-4 раза быстрее однопоточной, а в полтора-два раза медленнее.

Огосе! А вы не разбирались в причинах такого поведения? Бага в сырой 2010 или что? Какой-то совсем непонятный результат....

Я на профиле (VTune) вижу, что сильно больше времени начинает уходить на сохранение результата (в параллельной версии в сравнении с последовательной).
Там простой _mm_store_ps(&out[i],result), который транслируется в movaps

Возможно, причина в том, что это самое [i] в разных тредах очень синхронно попадает в один участок кэша у L2-кэша (L1 то у ядер свой). Разбираться - не очень понимаю как, это место в VTune (или, если хотите, в хардверных счетчиках у i7) я пока не понимаю.

Интел в этом месте делает, как уже обсуждали в соседних каментах, movntps, которые кэш не размывают.

Ага, понятно, хитрО...
Ну, т.е. если я правильно вас понял, то получается что причина действительно в плохом коде... Нда...

2010 студия весьма противоречивая получилась. Крайний раз я что-то похожее помню то ли на версии VS5, то ли на первой с .NET, которая после вполне нормальной VS6 была.

Да, местами она странная.

Но конкретно к movaps/movntps претензий нет: т.е. интел додумал в этом месте правильную инструкцию написать, но вообще надо там писать _mm_stream_ps (о чем я теперь знаю из гугла) и сейчас я это дело изучу...

Ооо, понятно. Знакомо)) Одна беда, что все спецификации, на основе которых надо работать, не начитаешься...

Спасибо за апдейт, ценное уточнение.

Наиболее чУдным открытием оказался ресет rounding mode внутри OpenMP-потока.

Вот это да - реальная жесть.

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

Не, ну я то изучал влияние размера интерполяционных таблиц на точность, как только оно стало попадать не в те строки из-за округления, так сразу стало понятно....

И, кстати, я подозреваю, что вообще все настройки FPU/SSE слетают.

Хм... Тогда они, вероятно, не слетают, а возвращаются к некоторым дефолтным настройкам. Значит наверное должен быть какой-то метод задания дефолтных настроек, это было бы логично..

Гугление показало следущее
- у PGI-компилятора слет состояния MXCSR в дочерних тредах считается багом и был исправлен в 10.5
- ничего другого здравого на тему 'OpenMP MXCSR' на первых страницах выдачи не нагуглилось.

MS просит возвращать как было перед вызовом других функций за исключением вариантов "явного согласия". Ну и defaults описаны в доке.

Подозреваю, что rounding mode вообще мало кто меняет.

Кстати, а ваши проекты сильно привязаны к кроссплатформенности? Новыми фишками VC++ не интересуетесь?)

Я тут вчера ВНЕЗАПНО узнал, что VC2010 уже поддерживает лямбда-выражения и что самое офигенское (ради чего пишу) - они предлагают собственную судя по описанию довольно могучую библиотеку Concurrency Runtime - возможности у неё просто офигенские, включая, например, совершенно простейшее распараллеливание циклов, переброс исключений между рабочими потоками, data flow driven архитектуры алгоритмов и ещё ваще дофига всего. Вот, например, как просто реализуются фьючерсы (х.з. как эта техника называется по русски, не встречал ещё, - суть в фоновом параллельном вычислении необходимых для алгоритма данных) - http://msdn.microsoft.com/en-us/library/dd764564.aspx
Это же просто праздник какой-то! :) Второй день челюсть от пола оторвать не могу))

Да, блин, наплодили добра.

У интела - тоже вот TBB, ABB, IPP и прочие ХЗ-ЧТО-ТАКОЕ

Но к кроссплатформенности в порядке Винда-Мак-Линукс я сильно привязан. А в случае с LibRaw - привязан еще сильнее, причем Linux - первее.

Т.е. на интеловские расширения еще можно смотреть, но и то, только на то что для мака есть, а MS пролетает.

Но за ссылки в любом случае спасибо, почитаю.

Ой.... Болять мои крылья.... (с)

Эти.....люди...Ну вот как их назвать? Фактически испортили OpenMP в VC2010 своей непревзойдённой гениальностью и знанием что для мира лучше: http://stackoverflow.com/questions/4738045/openmp-huge-performance-diffe... или http://social.msdn.microsoft.com/Forums/en/parallelcppnative/thread/3ba7...

Вкратце, чтобы не париццо со сложным реюзингом рабочих потоков, которые могут завершаться и начинаться в не слишком предсказуемые моменты времени, они просто перед убиванием каждого потока OpenMP добавили СПИН-ЛОК, мать их!!!, на 100мс!!! Мать-мать-мать... И в результате распараллеливание у чувака даёт мегазагруз процессора и адское время работы: 2 секунды в vc2005 против 3х минут vc2010. Ну как к этому отнестись?...

OpenMP это, канешна, не упомянутый выше Concurrency Runtime, но вот можно ли теперь этой самой ран-тайм доверять при таких гениальных решениях в смежной области - большой вопрос... Чую, мне теперь предстоит та самая эбля-с-перископом по таки выведению собственного thread pool (или если повезёт - успешной допилки около-бустового threadpool) для своих вычислений...

Может проще интеловский компилятор взять? Он не такой дорогой, окупится на экономии времени.

У интела в этом месте, с виду, ума побольше.

Да, явно придётся как-то изгаляться...

По первой ссылке (stackoverflow) пишут, впрочем, что у интела та же херь. Но регулируемая

ага, спасибо, видел. В том то и дело, что регулируемая...

Я вот всё пытаюсь понять, какими мыслями надо руководствоваться, чтобы принимать такие решения, как вот здесь они приняли с этой херью? Я б такого человека выгнал бы с работы с совершенно волчьей характеристикой, позорище какое-то просто... МС можно и есть за что ругать, но они всё равно одни из лучших производителей софта в мире. Но такие вот решения просто убивают...

Судя по тому, что интел принял то же решение - там есть какая-то идея.

Хотя мне эта идея непонятна совсем, но может быть какая-то типичная ошибка, вроде пропуска синхронизации (а так - оно в подавляющем случае синхронизируется само).

Чудеса какие-то.

Да, не, что там есть какая-то идея - понятно, но я про то, что МС сделали эту хрень ненастраиваемой. Вот где фейл...

Да, я тоже сути проблемы не понимаю, но претензия у меня только к отсутствию настройки)

Кстати, вот тут, чувак просто кипятком писает от счастья:

Earlier I was thinking that only your and other Microsoft people Demos fit perfectly for parallel solutions ! But indeed with some effort it is possible to make an application make the most of the processing power. Most of the work is already being done for us.

Так шта чую что пора уже заканчивать безуспешно гуглить и пора писать бенчмарк для обычного трид-пула и их Concurrency Runtime...

Не, ну OpenMP - тот еще уродец, особенно в сочетании с C++ (исключения, невозможность в #pragma omp описать член класса).
Зато "дешево, надежно и практично". Было. До этих фокусов со спинлоками.

И столь же понятно, что нужен какой-то более другой механизм, которого (общепринятого) нет.
Проблема в том, только, что массовый механизм будет рассчитан на идиотов, что печально.

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

Хе-хее, грабельки они такие грабельки...

Короче, возможно вам будет таки интересно, я сделал простенький бенчмарк на расчёте своих алгоритмов: Concurrency Runtime::parallel_for() vs. boost::threadpool 0.2.5 (с патчем утечки памяти и падения из их багтрекера и моими мелкими правками совместимости с VC2010).

Результаты на 6 ядрах такие:

Linear = 28237
Thread Pool 1 = 5460
Thread Pool 2 = 5428
Thread Pool 3 = 5523
Thread Pool 4 = 5476
Thread Pool 5 = 5413
Thread Pool 6 = 5444
ConcRT Pool 1 = 6911
ConcRT Pool 2 = 6708
ConcRT Pool 3 = 6771
ConcRT Pool 4 = 6832
ConcRT Pool 5 = 6771
ConcRT Pool 6 = 6614

Linear = 29187
Thread Pool 1 = 5367
Thread Pool 2 = 5398
Thread Pool 3 = 5428
Thread Pool 4 = 5398
Thread Pool 5 = 5382
Thread Pool 6 = 5398
ConcRT Pool 1 = 6692
ConcRT Pool 2 = 6771
ConcRT Pool 3 = 6630
ConcRT Pool 4 = 6708
ConcRT Pool 5 = 6849
ConcRT Pool 6 = 6723

Linear - это линейный алго, 1 прогон за тест, и по 6 прогонов на каждый из распараллеливающих алгоритмов. Числа после равенства - время в миллисекундах. Релизный билд.

Разница в 25% видна невооруженным глазом. Плюс бОльшая дисперсия времени в Concurrency Runtime. Плюс там же ещё не пофиксенные утечки памяти...

В общем, Concurrency Runtime юзать, конечно, можно, там побольше всяких вкусняшек, а утечку рано или позно исправят, но когда надо молотить данные без компромиссов - 25% это совсем не айс.

А пошто оно так тормозит - не видно? 25% должны бы быть видны профайлером...

Хорошая мысль, кстати, спасибо. Ща погляжу...

Ну вот если VCшный профилер не врёт, то получается около 13% прожирается только на следующих 6-ти функциях

Function Name Inclusive Samples Exclusive Samples Inclusive Samples % Exclusive Samples %
Concurrency::details::_SpinWait<1>::_SpinOnce(void) 18 332 17 699 7,83 7,56
Concurrency::_Parallel_chunk_helper,0>::operator()(void)const 233 698 10 170 99,88 4,35
Concurrency::details::InternalContextBase::IsSynchronouslyBlocked(void)const 3 578 3 577 1,53 1,53
Concurrency::details::ExternalContextBase::IsSynchronouslyBlocked(void)const 1 372 1 372 0,59 0,59
Concurrency::details::_HyperNonReentrantLock::_Acquire(void) 209 186 0,09 0,08
Concurrency::details::SchedulingRing::GetNextScheduleGroup(int *,int) 175 175 0,07 0,07

Особенно SpinWait нам как бе намекает... Вот накуа оно там???

ну да, на рабочую функцию получается приходится что-то около 85% всейнагрузки

да, забыл ещё сказать, что делается обсчёт довольно большого объёма данных, который скорее всего полностью в кеш не влазит. Соответственно, если потоки стартуют обсчёт практически одновременно, то неисключено, что они по данным идут примерно рядом друг с другом и соответственно, это вызывает более эффективное использование процессорного кеша и меньше фетчей из памяти. А когда работает через конкарренси, оно может стартовать потоки в соответствии со своими какими-то условиями, а за это время уже работающие потоки уйдут по данным далеко и возвращение к ним вызовет переобращение к памяти... Ну, где-то так в качестве гипотезы объяснения...
Ну и плюс эти спин-локи вообще не понятные... Я ещё могу понять, когда они в ядре используются и везде при этом в DDK аршинными буквами пишется: "выпускайте спин-локи как можно скорее!!!", а тут вот вдруг в юзермоде и такая хрень... ну как-то совсем не понятно... Будем надеятся, что они таки увидят этот юзкейс и поправят..

блин, промахнулся.
Итак, на рабочую функцию приходится около 85% всех семплов, соответственно - остальное - шум. Но от всего времени взять 85% - получается всё равно больше, чем время threadpool)))
Видимо ещё генерируется разный код, хотя рабочая функция вызывается всюду из одинаковой лямбды (впрочем, тут могут быть тонкости, связанные с передачей параметров в эту лямбду, которые могут по разному учитываться в профилировке и давать остальное смещение).
В общем, по моему дело пахнет киросином... Я канешна с этим добром по сути только сегодня познакомился, но пока у меня идей нет, куда дальше копать, чтобы это ускорить до приемлемых значений... А хотелось бы, ибо иначе придёццо самому возиться с синхронизацией доступа к результатам обсчёта, а это чот не хочется)
Ща ещё попробую profile guided optimization, может у неё там какая своя магия есть...

Я не знаю что у вас за задача, но если данных много, а доступ к ним (в пределах треда) локален, то зачем вы вообще фигачите близкие куски данных в разных тредах?

Если можно линейно читать (ну грубо, считаем сумму или там гистограмму), то лично я бы пилил на большие куски, может быть с хитро выставленным смещением, чтобы кэш L3 по возможности не размывался (т.е. ассоциативность в разных потоках кажет в разные места).

Ну и фигачить, фигатором. Оно или упрется в memory bandwidth (обсчет быстрее памяти) и тогда кэш не спасает и практически не нужен (для префетча только), или определяющим является именно скорость расчета, тогда значит в CPU упрется (и тогда рулит SSE).

спасибо за мысли, всё так...
Сейчас просто надо хоть как-то сделать, чтобы оно задышало и показало достижение целей, а уж потом "есичо" - оптимизировать буду... А не то я до того счастливого дня могу просто не дожить;) Грамотное распараллеливание потребует изменения алгоритма, а мне сейчас как воздух его понимаемость и изменяемость нужна. В линейном виде оно уж куда понятнее выглядит...
Хотя, вы конечно, правы, лучше это сразу закладывать, но я пока не уверен, что смогу, опыта в распараллеливании маловато и сейчас не самый подходящий момент, чтобы инвестировать время в это...
В любом случае, спасибо за диалог, чертовски приятно в кои-то веки поговорить со знающим своё дело человеком) Надеюсь, и вам какая-нибудь информация оказалась полезна.

Если задача пилится на крупные (десятки-сотни тысяч инструкций) куски, независимые по write-данным, то в распараллеливании никакой особой науки нет, берешь и параллелишь.

Если куски более мелкие, то надо смотреть, не сожрет ли оверхед на синхронизацию это самое распараллеливание.

Вот, кстати, и возможный ответ на загадку intel/MS: на спинлоке быстрее синхронизироваться, чем на чем-то еще (каком-то засыпании-просыпании).

А давать какие-то более конкретные советы, не зная специфики задачи - это газифицировать лужи.

кэша у L2-кэша (L1 то у ядер свой). Разбираться - не очень понимаю как, это место в VTune (или, если хотите, в хардверных счетчиках у i7)

У core i7 l2 не shared, l3 - shared.

http://www.tomshardware.com/reviews/Intel-i7-nehalem-cpu,2041-10.html

А что там с инвалидацией? Если она идет не по адресу, а по ассоциативному адресу, то даже с не-shared кэшем ядра будут портить малину брату.

Алексей, а не пробовали просто играться в 2010й студии
с /Oi и /arch:SSE2 и/или http://msdn.microsoft.com/en-us/library/bthd138d(v=VS.100).aspx
Вдруг существует какая-то неочевидная комбинация, которая работает под OpenMP быстрее.

Вы говорите о прямой имплементации, которая на C?
Ну да, в хорошем случае она догонит Intel, если цикл суммет векторизовать, чтобы значит powf считался для 4-х элементов массива сразу.

Впрочем, попробовал. Oi и SSE2 у меня были включены, добавление _set_SSE2_enable(1) ничего не меняет: 7 мегапикселей в секунду для одного потока и 36 для OpenMP.
Это примерно втрое медленнее интеловского варианта т.к. автовекторизации не случилось (втрое а не вчетверо т.к. есть ветвления и векторизатор далек от оптимальности).

С Visual C++ забавнее другое: в один поток SSE-вариант сплайнов быстрее, чем C-шный вариант сплайнов раза в полтора-два (не в четыре, потому что скалярный вариант просто считает, а мой векторный на ~8 операций счета делает еще 8 операций для транспонирования).

А вот с OpenMP - С-шный вариант ускоряется почти вчетверо (и впрямь, чего ему), а SSE - замедляется.

Я на досуге, конечно, попробую оторвать интеллект компилятору, возможно он увидев ассемблерные макросы внутри OpenMP-блока - дуреет и надо их от него попрятать в отдельную функцию в отдельной библиотеке.

Или треды свои самому запускать, по числу ядер....

Разбираться надо, конечно, Visual - гораздо более распространенный компилятор, чем Intel, а у меня опенсорс...

а просто таблица с прямым lookup - не?

насколько я слышал (точно не знаю) контроллеры управления двигателями в автомобилях ничего не вычисляют вообще - работают по таблицам значений f(x, y, z) если и не большего к-ва переменных

У таблицы на 65k записей и просто лукап - получается точность (максимальное отклонение от идеала) порядка 1% (потому что производная в некоторых местах - большая)

Если увеличивать таблицу - скорость будет падать, даже 65к записей - это уже целиком кэш L2.

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

Да, касаемо таблицы, скажем для замены табличным профилям. Там на входе - триада (RGB), в некоторых случаях хочется тетраду (RGBG2). Даже если 10 бит на элемент (что гадски мало даже для линейного лукапа), получается таблица с 2^30-2^40 входов. Как-то многовато.....

Другой вопрос, что точек для оценки/построения табличного профиля берется, в лучшем случае, несколько тысяч - значит это и будут реальные узлы таблицы. А дальше - интерполяция.