Некоторое время назад я писал про то, как получать воспроизводимые результаты и какие сложности с этим связаны. Также подробно рассказал про модели, позволяющие контролировать работу с числами с плавающей точкой в компиляторе и отдельно уточнил, что, если мы используем какие-либо библиотеки или стандарты, то должны позаботится, чтобы нужные флаги были указаны и для них. И вот совсем недавно я натолкнулся на интересную проблемку, связанную именно с воспроизводимостью результатов при работе с OpenMP.
Что такое воспроизводимость? Да всё просто – мы хотим получать одну и ту же «хорошую цифирь» от запуска к запуску, потому что для нас это важно. Это критично во многих областях, где сейчас активно используются параллельные вычисления.
Итак, как вы помните, для машинных вычислений существенную роль играет порядок суммирования, и, если у нас имеются циклы, распараллеленные с помощью любой технологии, то неизбежно возникнет проблема воспроизводимости результатов, потому что никто не знает в каком порядке будет проводиться суммирование, и на сколько «кусков» будет разбит наш исходный цикл. В частности это проявляется при использовании OpenMP в редукциях.
Рассмотрим простой пример.
!$OMP PARALLEL DO schedule(static)
do i=1,n
a=a+b(i)
enddo
В этом случае при использовании статического планировщика мы просто разобьем всё пространство итераций на равное количество частей и дадим каждому потоку выполнить эти итерации. Скажем, если у нас 1000 итераций, то при работе 4 потоков мы получим по 250 итераций «на брата». Наш массив является примером общих данных для разных потоков, и поэтому нам нужно позаботится о безопасности кода. Вполне рабочий вариант использовать редукцию и вычислять в каждом потоке своё значение, а затем складывать полученные «промежуточные» результаты:
!$OMP PARALLEL DO REDUCTION(+:a) schedule(static)
do i=1,n
a=a+b(i)
enddo
Так вот, даже на таком простом примере получить разброс в значениях можно достаточно просто.
Я поменял число потоков с помощью OMP_SET_NUM_THREADS и получил, что при 2 потоках a= 204.5992, а при 4 уже 204.6005. Способ инициализации массива b(i) и a я опустил.
Интересно то, что говорить о воспроизводимых результатах можно только при соблюдении целого ряда условий. Так вот, архитектура, ОС, версия компилятора, которым собиралось приложение и число потоков должно быть всегда постоянным от запуска к запуску. Если мы изменили число потоков, то результаты будут отличаться, и это абсолютно нормально. Тем не менее, даже при соблюдении всех этих условий результат всё равно может отличаться, и здесь нам должна помочь переменная окружения KMP_DETERMINISTIC_REDUCTION и статический планировщик. Оговорюсь, что её использование не даст нам гарантию совпадения результатов параллельной и последовательной версий приложения, равно как и с другим запуском, при котором использовалось отличное количество потоков. Это важно понимать.
Речь о достаточно узком случае, когда мы действительно ничего не меняли, а результаты не сошлись. И вот самый главный сюрприз заключается в том, что в некоторых случаях и KMP_DETERMINISTIC_REDUCTION не работает, хотя мы и «играли по правилам».
Такой код, который незначительно сложнее первого примера, даёт различные результаты:
!$OMP PARALLEL DO REDUCTION(+:ue) schedule(static)
do is=1,ns
do y=1,ny
do x=1,nx
ue(x,y)=ue(x,y) + ua(x,y,is)
enddo
enddo
enddo
!$OMP END PARALLEL DO
Даже после выставленной переменной KMP_DETERMINISTIC_REDUCTION, ничего не изменилось. Почему? Оказывается, в некоторых случаях компилятор из соображений производительности создаёт свою собственную реализацию цикла с использованием локов, и не заботится о результатах при этом. Эти случаи легко отследить по ассемблеру. В «хороших» вариантах обязательно должен быть вызов функции __kmp_reduce_nowait. А вот для моего примера подобного сделано не было, что несколько подрывает доверие к KMP_DETERMINISTIC_REDUCTION.
Итак, если вы написали код и результаты ваших вычислений скачут от запуска к запуску, не спешите посыпать голову пеплом. Отключите оптимизации и запустите ваше приложение. Проверьте, выравнены ли ваши данные и задайте «строгую» модель работы с числами с плавающей точкой. Для теста могут быть полезны следующие опции компилятора:
ifort -O0 -openmp -fp-model strict -align array32byte
Если и при таком наборе опций результаты вас удивляют, проверьте циклы, распараллеленные с помощью OpenMP и редукций и включите KMP_DETERMINISTIC_REDUCTION. Это может сработать и решить проблему. Если нет, то посмотрите на ассемблер и проверьте наличие вызова __kmp_reduce_nowait. В случае наличия этого вызова – проблема, вероятно, не с OpenMP и не с компилятором, а в код закралась ошибка. Кстати, проблему с KMP_DETERMINISTIC_REDUCTION мы должны решить в скором времени. Но учитывайте эту особенность уже сейчас.
Комментарии (22)
pehat
08.07.2015 11:56+8>Цифирь не сошлась
>Блог компании Intel
Ожидал увидеть здесь древнюю историю про баг в пентиуме.
grafmishurov
08.07.2015 12:47+1У меня сейчас gcc на Celeron (Bay Trail) флоаты для деления и умножения в XMM-регистры кладет и с SSE считает, а не в FPU.
ivorobts Автор
08.07.2015 13:55+1Это правильно, сейчас FPU почти не используется.
VioletGiraffe
08.07.2015 14:27Почему? Неужели так быстрее? Не вижу логики.
0serg
08.07.2015 15:20Классический FPU это стековая архитектура (привет наследию x87), SSE — регистровая
С регистрами часто работать гораздо быстрее чем со стекомgrafmishurov
09.07.2015 22:02В FPU регистровый же стек, а не RAM.
0serg
09.07.2015 22:18Так загружать-выгружать данные то из этого стека куда-то надо. А поскольку старый набор инструкций FPU был создан когда сопроцессор был отдельным чипом, то FST/FIST выгружает данные из стека сугубо в память. И загрузкой та же история — данные можно брать только из памяти. Ну и в целом, как ни крути, неудобно манипулировать со стеком. Даже с учетом FXCH все равно нужно туда-сюда данные крутить, поскольку почти все инструкции работают сугубо с вершиной стека. Кроме того в суперскалярах то что практически все инструкции используют один и тот же регистр порядком убивает производительность, поскольку две операции с FPU одновременно выполнить уже становится проблемой.
ivorobts Автор
08.07.2015 15:33Всё просто — так быстрее!
VioletGiraffe
08.07.2015 15:44-1Понятно, что быстрее. Непонятно, почему векторная подсистема при использовании не по назначению быстрее, чем FPU при использовании по назначению. Там же такие сложности с загрузкой и чтением XMM-регистров — латентность, требования по выравниванию памяти и т. д.
0serg
08.07.2015 16:07+1Нет там сложностей с загрузкой и чтением XMM регистров. Многие SSE команды имеют векторный и скалярный варианты и скалярный имеет меньше ограничений чем векторный
Применительно к чтению XMM, чтение одного значения в SSE регистр — это MOVSS/MOVSD и там нет требований к выравниванию. Чтение вектора требующее выравнивания — это другая инструкция MOVAPS/MOVAPD и у нее есть и невыровненный эквивалент MOVUPS/MOVUPD, просто он более медленный
grafmishurov
08.07.2015 13:01Не, ну тут действительно не будет один и тот же резульатат, много факторов, которые вне контроля программиста: состояния в L-кешах, какие из тредов ядра/процессоры немного раньше или немного позже посчитают. В данном случае важно знать допустимую погрешность при вычислениях, чтобы ее учитывать при разработке алгоритма.
9mm
08.07.2015 12:26-1А не правильнее ли здесь использовать алгоритмическое повышение точности, через Pairwise_summation или Алгоритм Кэхэна?
devlev
08.07.2015 13:23+2Я не силен в фортране но при помощи подручной браузерной консоли удалось получить ситуацию как у вас (и без параллельного программирования) — не выполнения ассоциативности сложения.
var a = 1000; var b = 0.0004; console.log(((a + b) + b) == (a + (b + b)));
var a = 10000000; var b = 0.0004; console.log(((a + b) + b) == (a + (b + b)));
В свое время когда писал магистерскую на CUDA, я тоже наткнулся на эту проблему. Эта была программа по вычислению гауссова размытия изображения и тоже заметил что свертка матрицы на GPU и CPU немного отличались. На тот момент я решил эту проблему. Нужно было чтобы порядок операций над числами на GPU и CPU совпадал. Условно: ((a + b) + (a + b)) если я выполню это на CPU или если я выполню на GPU отдельно сложения в скобках а потом просуммирую результат на CPU то результат будет тот же.
Поправьте, если я не правivorobts Автор
08.07.2015 13:57+3Всё правильно — ассоциативность не выполняется. Но вы должны учитывать, что кроме написанного вами кода, есть ещё компилятор, который может его изменить, и библиотеки, которые могут не учитывать вашу необходимость в «правильном» порядке суммирования. Про это и разговор.
devlev
08.07.2015 15:21+1Вообще как мне кажется основная тема вашей публикации — это как раз не выполнение этой ассоциативности. Вы вроде указали даже команды которые отключают оптимизацию кода компилятора но ни слова не сказали о главной сути проблемы. Решение проблемы не понимая ее причины обычно не приводит к успеху. В обычном кольце всегда выполняет равенство 2 + 2 = 4, а в программировании мы можем сделать и так как у вас на рисунке. Чисто алгоритмически процессоры на GPU и CPU считают одинаково. При решении своей задачи, у меня сверки матриц были идентичными во всех случаях, но стоило заменить последовательность вычисления как результат получался уже немного другой.
ivorobts Автор
08.07.2015 15:37+1В самом начале я давал ссылку на статью, где проблема описана намного шире и детальнее. Здесь речь в основном о частном случае, а именно использовании OpenMP и возможностей по контролю воспроизводимости там.
saluev
Хм, вообще относительная ошибка порядка 1E-7 (как в вашем случае) при работе с float'ами — это нормально. Она такая примерно всегда получается. Просто если у вас вещественные вычисления, то воспроизводимость результатов надо мерить не тупым равенством, а приближённым.
grafmishurov
Это НЕ нормально, урезание количества знаков после запятой работать должно по одному и тому же алгоритму и выдавать всегда один и тот же (floor, round, ceil) результат, тем более, что они тестировали с тем же количеством нитей, выполняющих задачу. Фортран я не знаю, а вот код на C в подобных случаях всегда имело бы смысл сначало транслировать в ассемблерный код (gcc -S), его анализировать, а потом уже в бинарник, ну или для чистоты вообще реверсить с objdump, чтобы совсем точно.
ivorobts Автор
Алгоритм «урезания» заложен в стандарте IEEE 754, и компилятор его поддерживает. Но если вы посмотрите пост, на который я ссылался в самом начале, то увидите пример, когда всё в соответствии с стандартом, а результат зависит от порядка суммирования, потому что критичны случаи, когда складывается большое и очень маленькое число.
grafmishurov
Ну в Вашем случае оптимизация в компиляторе наотсебятила. Я вот тут еще нашел про выравнивание и –fp-model (or /fp) flag в «Using the Intel® Math Kernel Library (Intel® MKL) and Intel® Compilers to Obtain Run-to-Run Numerical Reproducible Results by Todd Rosenquist, Technical Consulting Engineer, Intel® Math Kernal Library and Shane Story, Manager of Intel® MKL Technical Strategy, Intel Corporation»
grafmishurov
П.С. Извиняюсь, в посте по ссылке все это есть. Прочитал пост после того, как предыдущий комментарий отправил.