company_banner

Сказ о том, как «цифирь» не сошлась



    Некоторое время назад я писал про то, как получать воспроизводимые результаты и какие сложности с этим связаны. Также подробно рассказал про модели, позволяющие контролировать работу с числами с плавающей точкой в компиляторе и отдельно уточнил, что, если мы используем какие-либо библиотеки или стандарты, то должны позаботится, чтобы нужные флаги были указаны и для них. И вот совсем недавно я натолкнулся на интересную проблемку, связанную именно с воспроизводимостью результатов при работе с 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 мы должны решить в скором времени. Но учитывайте эту особенность уже сейчас.
    Intel
    Company

    Comments 22

      –4
      Хм, вообще относительная ошибка порядка 1E-7 (как в вашем случае) при работе с float'ами — это нормально. Она такая примерно всегда получается. Просто если у вас вещественные вычисления, то воспроизводимость результатов надо мерить не тупым равенством, а приближённым.
        +1
        Это НЕ нормально, урезание количества знаков после запятой работать должно по одному и тому же алгоритму и выдавать всегда один и тот же (floor, round, ceil) результат, тем более, что они тестировали с тем же количеством нитей, выполняющих задачу. Фортран я не знаю, а вот код на C в подобных случаях всегда имело бы смысл сначало транслировать в ассемблерный код (gcc -S), его анализировать, а потом уже в бинарник, ну или для чистоты вообще реверсить с objdump, чтобы совсем точно.
          0
          Алгоритм «урезания» заложен в стандарте IEEE 754, и компилятор его поддерживает. Но если вы посмотрите пост, на который я ссылался в самом начале, то увидите пример, когда всё в соответствии с стандартом, а результат зависит от порядка суммирования, потому что критичны случаи, когда складывается большое и очень маленькое число.
            0
            Ну в Вашем случае оптимизация в компиляторе наотсебятила. Я вот тут еще нашел про выравнивание и –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»

            The Intel compiler also
            provides threading via the OpenMP* model. With the latest compiler,
            the reduction stage in threaded parallel sections can be forced to
            provide reproducible results by setting
            KMP_DETERMINISTIC_REDUCTION=yes.
            Let’s consider the case where you would like the resulting executable
            to produce the same results on the following three systems:
            > One 4-core processor supporting SSE 4.2 instructions
            > One 4-core processor supporting AVX instructions
            > Two 4-core processors supporting AVX instructions
            First, to ensure that Intel MKL functions provide reproducible results,
            we must make sure all arrays are aligned on 16-byte boundaries (the
            mkl_alloc() function suits this purpose). Then, ensure that the number
            of threads used on each system remains fixed and does not vary
            during the run. By default, Intel MKL uses as many threads as there
            are cores, so you should fix the number of threads at four using the
            mkl_set_num_threads() function. Finally, because, Intel MKL dispatches
            an optimized code path based on the instruction set available, you
            need to use the new reproducibility control to configure this: mkl_
            cbwr_set(MKL_CBWR_SSE4_2. You will then also need to use an
            appropriate –fp-model (or /fp) flag to ensure the compiler returns
            consistent results. Applications threaded using OpenMP should also
            specify a fixed number of threads (omp_set_num_threads(4)) and set
            KMP_DETERMINISTIC_REDUCTION=yes.
              0
              П.С. Извиняюсь, в посте по ссылке все это есть. Прочитал пост после того, как предыдущий комментарий отправил.
        +8
        >Цифирь не сошлась
        >Блог компании Intel
        Ожидал увидеть здесь древнюю историю про баг в пентиуме.
          +1
          У меня сейчас gcc на Celeron (Bay Trail) флоаты для деления и умножения в XMM-регистры кладет и с SSE считает, а не в FPU.
            +1
            Это правильно, сейчас FPU почти не используется.
              0
              Почему? Неужели так быстрее? Не вижу логики.
                0
                Классический FPU это стековая архитектура (привет наследию x87), SSE — регистровая
                С регистрами часто работать гораздо быстрее чем со стеком
                  0
                  В FPU регистровый же стек, а не RAM.
                    0
                    Так загружать-выгружать данные то из этого стека куда-то надо. А поскольку старый набор инструкций FPU был создан когда сопроцессор был отдельным чипом, то FST/FIST выгружает данные из стека сугубо в память. И загрузкой та же история — данные можно брать только из памяти. Ну и в целом, как ни крути, неудобно манипулировать со стеком. Даже с учетом FXCH все равно нужно туда-сюда данные крутить, поскольку почти все инструкции работают сугубо с вершиной стека. Кроме того в суперскалярах то что практически все инструкции используют один и тот же регистр порядком убивает производительность, поскольку две операции с FPU одновременно выполнить уже становится проблемой.
                  0
                  Всё просто — так быстрее!
                    –1
                    Понятно, что быстрее. Непонятно, почему векторная подсистема при использовании не по назначению быстрее, чем FPU при использовании по назначению. Там же такие сложности с загрузкой и чтением XMM-регистров — латентность, требования по выравниванию памяти и т. д.
                      +1
                      Нет там сложностей с загрузкой и чтением XMM регистров. Многие SSE команды имеют векторный и скалярный варианты и скалярный имеет меньше ограничений чем векторный

                      Применительно к чтению XMM, чтение одного значения в SSE регистр — это MOVSS/MOVSD и там нет требований к выравниванию. Чтение вектора требующее выравнивания — это другая инструкция MOVAPS/MOVAPD и у нее есть и невыровненный эквивалент MOVUPS/MOVUPD, просто он более медленный
              0
              Не, ну тут действительно не будет один и тот же резульатат, много факторов, которые вне контроля программиста: состояния в L-кешах, какие из тредов ядра/процессоры немного раньше или немного позже посчитают. В данном случае важно знать допустимую погрешность при вычислениях, чтобы ее учитывать при разработке алгоритма.
                0
                П.С. Не в ту ветку последний комментарий, обращался к ОПу.
              –1
              А не правильнее ли здесь использовать алгоритмическое повышение точности, через Pairwise_summation или Алгоритм Кэхэна?
                +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)));

                image

                В свое время когда писал магистерскую на CUDA, я тоже наткнулся на эту проблему. Эта была программа по вычислению гауссова размытия изображения и тоже заметил что свертка матрицы на GPU и CPU немного отличались. На тот момент я решил эту проблему. Нужно было чтобы порядок операций над числами на GPU и CPU совпадал. Условно: ((a + b) + (a + b)) если я выполню это на CPU или если я выполню на GPU отдельно сложения в скобках а потом просуммирую результат на CPU то результат будет тот же.

                Поправьте, если я не прав
                  +3
                  Всё правильно — ассоциативность не выполняется. Но вы должны учитывать, что кроме написанного вами кода, есть ещё компилятор, который может его изменить, и библиотеки, которые могут не учитывать вашу необходимость в «правильном» порядке суммирования. Про это и разговор.
                    +1
                    Вообще как мне кажется основная тема вашей публикации — это как раз не выполнение этой ассоциативности. Вы вроде указали даже команды которые отключают оптимизацию кода компилятора но ни слова не сказали о главной сути проблемы. Решение проблемы не понимая ее причины обычно не приводит к успеху. В обычном кольце всегда выполняет равенство 2 + 2 = 4, а в программировании мы можем сделать и так как у вас на рисунке. Чисто алгоритмически процессоры на GPU и CPU считают одинаково. При решении своей задачи, у меня сверки матриц были идентичными во всех случаях, но стоило заменить последовательность вычисления как результат получался уже немного другой.
                      +1
                      В самом начале я давал ссылку на статью, где проблема описана намного шире и детальнее. Здесь речь в основном о частном случае, а именно использовании OpenMP и возможностей по контролю воспроизводимости там.

                Only users with full accounts can post comments. Log in, please.