О точности вычислений с плавающей точкой

Недавно я тестировал производительность разных компьютеров на предмет быстроты и стабильности операций с вещественной переменной в цикле (ну кому не хотелось бы найти стабильную и производительную систему для рутинных вычислений?). Запускались 3 тестовые программы: (S) вычисление суммы целых чисел в интервале [-N,N] (при суммировании в цикле целые числа преобразовывались в вещественные); (E) вычисление числа E = (1+1/N)**N; (PI) вычисление числа PI = \int_0^1 dx 4/(1+x**2). Во всех тестах для разных N измерялось время вычисления соответствующей величины. Тест S должен показывать насколько быстро и точно компьютер получит сокращение всех членов прямым суммированием без группировки членов и при каких значениях N он начнет ошибаться. Тест E должен измерять производительность относительно операции произведения, а тест PI - это комбинированный тест уже близкий к тому, что делается в реальных программах.

Я написал простые программки на Фортране и запускал их на разных компьютерах, включая разные узлы кластера, меняя N и разрядность вещественных переменных. Результаты этого упражнения мне показались забавными и, рискну предположить, небезЫнтересными. В частности, обнаружилось, что вычисления с вещественными переменными расширенной точности (80 бит) не только точнее аналогичных вычислений со стандартными переменными двойной точности (64 бит), что естественно, но и во многих случаях еще и быстрее, что оказалось довольно неожиданным. Подробности обсуждаются ниже.

Точность вычислений: 32-бит vs. 64-бит vs. 80-бит

Обсудим сначала результаты теста S (файлы можно найти здесь). Я использовал компилятор gfortran и запускал тестовый код на компьютерах с Linux декларируя разную длину вещественной переменной S. Очевидно, мы должны получить S=0 для любых значений N. Результаты расчетов на одном 64-битном процессоре показаны ниже. Для наглядности, на картинке ниже нарисованы результаты для lg(1+|S|) в зависимости от lg N. Врезка слева вверху показывает, что получается при распараллеливании цикла на 40 процессорах (см. обсуждение ниже).

В Фортране вещественная переменная по умолчанию [real=real(4)] имеет длину 4 байта = 32 бит. Какая точность у 32-битной real(4) переменной? Число значимых цифр 32-битной переменной, казалось бы, можно оценить как lg(2**32)~9.5. Однако, 8 бит занимают порядок и знак и фактическая значимая длина 32-битной вещественной переменной (мантисса) составляет 23 бита, т.е. только 6 цифр. Для переменных этого типа код правильно вычисляет S до N=2**12 и начинает ошибаться начиная с N=2**13. Очевидно, сумма в цикле сначала собирается из членов одного знака и достигает значения S~N**2/2 и только после этого к S начинают добавляться члены другого знака. При N=2**13 старшие значения S в цикле уже не помещаются в 32-битную разрядность, код начинает округлять S, ошибки округления накапливаются, и в результате мы не получаем точного сокращения всех членов.

В случае 64-битной вещественной переменной двойной точности типа real(8), мантисса имеет длину 53 бит, т.е. только около 16 цифр являются значимыми. Компьютер правильно вычисляет S до N=2**26, а начиная с N=2**27 значения S в цикле не помещаются в 64-битную разрядность, и код начинает ошибаться.

В статье википедии про двойную точность есть ссылка на формат чисел с расширенной двойной точностью длиной 80 бит. Оказывается, что процессоры x86 и x86-64 поддерживают этот формат на уровне железа. Даже более того, 80-битовый формат вещественных чисел является естественным для математического сопроцессора x87, который занимается арифметикой. Сопроцессор проводит внутренние (промежуточные) вычисления в 80-битовом формате, для уменьшения ошибок округления и повышения точности конечного результата. Если ему, например, подсовывают 64-битную переменную, то он ее преобразовывает в 80-битный формат. Причем он должен делать это дважды: на входе и на выходе. Было интересно проверить, что получится если явно декларировать S в 80-битовом формате. Мы имеем все основания ожидать, что точность вычислений повысится.

Вещественные переменные расширенной точности декларируются как real(10). Фактическая точность этого формата (т.е. число значимых цифр) отличается от наивной оценки lg(2**80)~24 и близка к lg(2**64)~19 знакам. Как и ожидалось, точность вычислений в 80-битной версии существенно повышается, код достойно держится до N=2**32, но начинает ошибаться при больших N. Как видно из графика, переход к 80-битным переменным позволяет расширить допустимую область N на 6 порядков относительно 32-битной версии, и почти на 2 порядка относительно 64-битной версии.

Для теста параллельных вычислений использовался OpenMP. Производительность системы масштабируется в соответствии с числом задействованных процессорных нитей (см. обсуждение ниже). Отметим также, что параллельное выполнение суммирования в цикле также благоприятно сказывается на точности сокращения членов разного знака при больших N. Выполнение параллельного кода на 40 ядрах позволяет добиться "почти" точного результата для 64-битной версии для N=2**28 (сравниваем с N=2**26 для 1 ядра). Это связано, по-видимому, с тем что в парциальных суммах, которые выполняются параллельно, накапливается меньше ошибок округления из-за уменьшения количества членов в сумме S и, как следствие, уменьшение абсолютной величины S (вспоминаем, что S растет как N**2). Для 80-битной версии результаты выглядят еще оптимистичнее. Как уже отмечалось выше, точный результат получается для N=2**32 уже на 1 ядре. На 40 ядрах мы имеем почти точное сокращение S для N=2**34.

Относительная точность вычисления произведения E и интеграла PI для разного типа переменных иллюстрируется на картинках ниже, на которых показана относительная погрешность вычисления в зависимости от N. В частности, для числа Е погрешность y=|E-E26|/E26, где E26 это число Е точное до 26 знаков, и аналогично для числа PI. Чтобы не потерять точность, погрешность y вычислялась в 128-битной арифметике.

Как видно из этих графиков, погрешность вычисления E невозможно сделать ниже чем 10**-5 в случае 32-битных переменных. Аналогичные пределы для 64 и 80 бит составляют 10**-11 и 10**-13, соответственно. Погрешность вычисления PI несколько лучше и примерно соответствует числу значимых цифр в тесте S. Более низкая точность вычисления E, по-видимому, вызвана накоплением ошибки округления 1+1/N в произведении при больших N.

Производительность для разных типов вещественных переменных

Нет ничего удивительного в том, что при переходе от 64-битным к 80-битным переменным повышается точность вычислений. Однако, неожиданным и забавным оказалось то, что код с переменной real(10) может быть быстрее аналогичного кода с переменной real(8)! Заменив только тип переменных, я наблюдал ускорение кода S почти в 2 раза на машине cluster с процессором Intel Xeon 3 GHz (отпрыск 2005 года семьи Irwindale). Такое наблюдение заставляет почесать затылок и попытаться воспроизвести эффект на других машинах, более современных процессорах и заодно проверить не является ли это глюком компиляторного происхождения или программы измерения времени.

Мои усилия суммированы в следующих графиках, которые показывают относительное изменение производительности (т.е. отношение обратных времен затраченных на получение результата) при изменении разрядности вещественных переменных от 64 до 80 в разных тестах. Все машины были с 64-битными процессорами в среде 64-битных версией Linux (за исключением nuke, на котором установлена 32-битная версия Linux). Список компьютеров и их основных/существенных параметров см. в таблице ниже. Время измерялось с помощью функции omp_get_wtime() из OpenMP. Результаты измерений довольно сильно флуктуируют при небольших значениях N (малые времена выполнения кода), но стабилизируются при больших N > 10**7. Для того чтобы исключить случайные флуктуации, один и тот же код выполнялся несколько раз, после чего вычислялось среднее по запускам время. На графиках ниже представлены результаты для средних времен. Каждая панель проказывает результаты собранные на одном компьютере, на котором выполнялись программы S, E и PI.

Нужно отметить, что производительность кода существенно зависит от компилятора и выбора параметров компиляции. Флагов у компилятора gfortran много и изучение их возможностей на предмет оптимизации может легко занять все имеющееся время. Я проводил измерения времени работы кода с двумя наборами оптимизирующих флагов: gfortran -fopenmp -O2 (верхний ряд на рисунке ниже) и gfortran -fopenmp -O2 -march=native (нижний ряд). Такой выбор параметров оптимизации кажется мне естественным: -O2 осуществляет стандартную общую оптимизацию кода и этот уровень оптимизации рекомендован мануалом, а флаг -march=native пытается дополнительно оптимизировать код под конкретную x86 архитектуру процессора.

___________________________________________________________________________________________
System      CPU model                          No.CPUs    Linux          kernel    gfortran
___________________________________________________________________________________________
cluster     Intel Xeon 3.0GHz                      4      SL x86_64      2.6.32-*   4.5.2
cluster42   AMD Opteron 6168 1.9GHz                24     SL x86_64      2.6.32-*   4.4.6
cluster51   AMD Opteron 6276 2.3GHz                64     SL x86_64      2.6.32-*   4.4.7
rut         Intel Core i5-2410M 2.30GHz            4      Debian x86_64  3.2.0-4    4.7.2
nuke        Intel Dual-Core Pentium E5400 2.7GHz   2      Debian i686    3.2.0-4    4.7.2
chromebook  Intel Celeron N2840 2.16GHz            2      Debian x86_64  3.10.18    4.9.2
___________________________________________________________________________________________

Картинки получились достаточно веселыми. Во-первых нужно отметить, что в целом, независимо от процессора, 80-битная версия кода, по крайней мере, не медленнее 64-битной, что уже очень полезно в приложениях, поскольку 80-битная арифметика точнее. Во-вторых, в некоторых тестах 80-битная версия работает заметно быстрее. В этом контексте отметим старика cluster, который оказался чемпионом по разгону сложения в тесте S, но, правда, вел себя довольно индифферентно по отношению к тестам E и PI. Оптимизация под конкретный процессор с помощью -march=native почти всегда приводила к более быстрому коду как для 64- так и для 80-битных версий. Что касается относительной 80/64 производительности, то отметим значительный эффект ускорения вычисления интеграла в тесте PI для процессоров AMD Opteron 6276 (cluster51) и Intel Core i5-2410M (rut) (см. нижний набор панелей на рисунке). Двукратное ускорение машины cluster51 в 80-битной версии в тесте PI кажется удивительным на фоне ее безразличия к другим тестам. Процессор AMD Opteron 6168 (cluster42) также быстрее работает в 80-битной версии теста PI. Отметим также, что производительность в тесте Е оказалась систематически безразличной к разрядности вещественных переменных.

Мой коллега переписал код теста S на C и также наблюдает эффект ускорения для 80-битных переменных с похожим соотношением 80/64 времен. Использовался gcc и уровень оптимизации O2. Это является косвенным подтверждением, что наблюдаемое ускорение 80-битного варианта кода не есть специфический эффект gfortran. Правда, gcc и gfortran тесно связаны, и "ускорение 80" может быть особенностью "гнутого" компилятора. В этом контексте было бы интересно проверить другие компиляторы, в первую очередь Intel Fortran. Однако, Intel больше не раздает ученым бесплатные лицензии на свои компиляторы... Если кто-то возьмется проверить эти наблюдения на С/С++ (и на других компиляторах) я был бы рад увидеть результаты.

Зависимость производительности от типа процессора

Раз уж мы взялись измерять времена, интересно выяснить относительную производительность разных систем. На графике ниже нарисована производительность тестовых кодов на разных процессорах относительно машины cluster (напомню, что под производительностью понимается обратное время затраченное на производство результата).

Как видно из графиков, даже скромные, но более новые системы работают быстрее кластерных узлов. Ноутбук rut оказался в этом тесте в 2 раза быстрее сервера cluster, несмотря на то что частота процессора rut ниже (см. таблицу параметров машин выше). Забавно, что новый хромбук, который был поставлен на эти тесты скорее для забавы, также довольно шустро складываает и перемножает, но отстает в вычислении интеграла PI. Однако, в поддержку х-бука отмечу, что это полностью 64-битная система с, хоть и урезанным, но новым, процессором. В боевом режиме процессор N2840 работает на частоте 3.1GHz (можно подсмотреть в /proc/cpuinfo), что выше чем частота процессора cluster. Также на нем установлен свежий компилятор gfortran v.4.9. Можно предположить, что этот компилятор оптимизирует код лучше, чем старые компиляторы кластерных машин. В этом контексте, было бы интересно протестировать вычислительные возможности процессоров ARM в планшетах и телефонах. Может быть уже пора попытаться использовать телефоны не только по их прямому назначению :-\?

Параллельные вычисления с 80-битными переменными

Описанные выше тесты проводились на одном ядре процессора. Очевидно, что преимущество кластера состоит скорее в большом количестве ядер и в возможности проводить параллельные вычисления, а не в индивидуальной скорости процессора (ну, правда, неплохо иметь быстрое ядро в любом случае). Поэтому, важно проверить, что мы выиграем в производительности и точности кода с 80-битными переменными при попытке распараллелить программу. Для параллельного теста я использовал OpenMP. Результаты тестов суммированы на графике ниже. Здесь показано отношение производительности (т.е. обратных времен) для заданной системы для параллельной и последовательной версий кода.

Эффект параллельности впечатляет! Для теста Е производительность увеличивается пропорционально числу задействованных нитей (threads). Для тестов S и PI наблюдаемый прирост несколько меньше для зеленых и красных машинок. Синяя машинка cluster42 в этом отношении ведет себя более стабильно. Возможно, это связано с тем, что у этой системы число нитей = числу задействованных физических ядер. Зеленые и красные могут выпустить по 2 нити из каждого физического ядра и с этим могут быть связаны наблюдаемые изменения в производительности параллельных вычислений для разных тестов. Как отмечалось выше, параллельное выполнение суммирования в тесте S также благоприятно сказывается на точности сокращения членов разного знака при больших N (см. врезку на рис.1).

Времена в параллельных тестах довольно сильно флуктуируют для узлов кластера. Тем не менее, сравнение времен 64/80 в разных тестах показывает, что 80-битная версия работает по крайней мере не медленнее 64-битной, а часто заметно быстрее. Так, например, в тесте PI выигрыш в производительности при переходе 64 -> 80 на узле cluster51 составляет 40-60% и порядка 20-30% для узла cluster42.

Summary

Using the GNU Fortran compiler (gfortran) we discuss accuracy and performance of numerical calculations on a few different x86_64 Linux systems. To probe the systems we use the three simple test codes: (Test S) calculation of the sum S over integers in the interval [-N,N]; (Test E) calculation of E=(1+1/N)**N by a direct evaluation of the product; (Test PI) calculation of PI=\int_0^1 dx 4/(1+x**2) by splitting the integration interval into N rectangles. We calculate the quantities S, E, and PI as a function of N for different settings of the real floating point variable (32-, 64-, and 80-bit length) and study convergence and performance of the codes. The corresponding timings were measured using the OpenMP function. A detailed comparison of the test results reveal a surprising effect of using the 80-bit floating point variables: not only the 80-bit version of a code is more accurate but also faster than the 64-bit version. The corresponding gain in performance (elapsed time) may be as large as the factor of 2. We discuss possible underlying reason for this observation. We also test this effect in parallel multi-processor calculations with OpenMP.

Source codes, scripts and data files with the results of the tests can be found here.