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

Недавно я тестировал производительность разных компьютеров на предмет быстроты и стабильности суммирования вещественной переменной в цикле (ну кому не хотелось бы найти стабильную и производительную систему для рутинных вычислений, с возможностью параллелить код?). Тест проводился на Фортране с помощью программы, в которой вычислялась сумма S по целым числам m в интервале [-N,N], и измерялось время суммирования в зависимости от N. Результат нам известен заранее -- мы должны иметь S=0 при любом N. При суммировании в цикле целые числа m преобразовывались в вещественные. Было интересно проверить насколько быстро и точно компьютер получит 0 прямым суммированием без группировки членов и при каких значениях N он начнет ошибаться. Результаты этого упражнения мне показались забавными и, рискну предположить, небезЫнтересными, и я решил изложить их в этой небольшой заметке.

Вещественные переменные: 32-бит vs. 64-бит vs. 80-бит

Я использовал компилятор gfortran и запускал тестовый код на компьютерах с Linux декларируя разную длину вещественной переменной S. В Фортране вещественная переменная [real=real(4)] по умолчанию имеет длину 4 байта = 32 бита. Какая точность у 32-битной переменной типа real(4)? Число значимых цифр 32-битной переменной, казалось бы, можно оценить как lg(2**32)~9.5. Однако, как нас учит википедия, фактическая значимая длина 32-битной вещественной переменной составляет 23 бита, т.е. только 6 цифр (8 бит занимают экспонента и знак). Для переменных этого типа код правильно вычисляет 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-битовый формат вещественных чисел является естественным для математического сопроцессора 387. Сопроцессор проводит внутренние (промежуточные) вычисления в 80-битовом формате, для уменьшения ошибок округления и повышения точности конечного результата. Если ему, например, подсовывают 64-битную переменную, то он ее преобразовывает в 80-битный формат. Причем он должен делать это дважды: на входе и на выходе. Было интересно проверить, что получится если явно декларировать S в 80-битовом формате. Мы имеем все основания ожидать, что точность вычислений повысится.

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

 p       N = 2**p                S 32-bit                  S 64-bit                   S 80-bit
 12  0.409600000000E+04    0.0000000000000000E+00    0.0000000000000000E+00    0.0000000000000000E+00
 13  0.819200000000E+04   -0.5793000000000000E+04    0.0000000000000000E+00    0.0000000000000000E+00
 14  0.163840000000E+05   -0.1531100000000000E+05    0.0000000000000000E+00    0.0000000000000000E+00
 15  0.327680000000E+05   -0.3213400000000000E+05    0.0000000000000000E+00    0.0000000000000000E+00
 16  0.655360000000E+05   -0.6489600000000000E+05    0.0000000000000000E+00    0.0000000000000000E+00
 17  0.131072000000E+06   -0.4896000000000000E+04    0.0000000000000000E+00    0.0000000000000000E+00
 18  0.262144000000E+06   -0.2595520000000000E+06    0.0000000000000000E+00    0.0000000000000000E+00
 19  0.524288000000E+06   -0.5058400000000000E+06    0.0000000000000000E+00    0.0000000000000000E+00
 20  0.104857600000E+07   -0.9994160000000000E+06    0.0000000000000000E+00    0.0000000000000000E+00
 21  0.209715200000E+07   -0.1769468000000000E+07    0.0000000000000000E+00    0.0000000000000000E+00
 22  0.419430400000E+07   -0.3670014000000000E+07    0.0000000000000000E+00    0.0000000000000000E+00
 23  0.838860800000E+07    0.1000000000000000E+01    0.0000000000000000E+00    0.0000000000000000E+00
 24  0.167772160000E+08    0.1094587700000000E+08    0.0000000000000000E+00    0.0000000000000000E+00
 25  0.335544320000E+08   -0.3355443200000000E+08    0.0000000000000000E+00    0.0000000000000000E+00
 26  0.671088640000E+08    0.7275248458137600E+15    0.0000000000000000E+00    0.0000000000000000E+00
 27  0.134217728000E+09    0.2251799813685248E+16   -0.1158500000000000E+05    0.0000000000000000E+00
 28  0.268435456000E+09    0.4503599627370496E+16   -0.2324719230000000E+09    0.0000000000000000E+00
 29  0.536870912000E+09    0.9007199254740992E+16   -0.5198229530000000E+09    0.0000000000000000E+00
 30  0.107374182400E+10    0.1801439850948198E+17   -0.1065319903000000E+10    0.0000000000000000E+00
 31  0.214748364800E+10    0.3602879701896397E+17   -0.2143284860000000E+10    0.0000000000000000E+00
 32  0.429496729600E+10    0.7205759403792794E+17   -0.4292867744000000E+10    0.0000000000000000E+00
 33  0.858993459200E+10    0.1441151880758559E+18   -0.8588876000000000E+10   -0.6074000997000000E+10
 34  0.171798691840E+11    0.2882303761517117E+18   -0.1717933977200000E+11   -0.1607029610300000E+11

Как видно из этой таблицы, переход к 80-битным переменным позволяет расширить допустимую область N на 6 порядков относительно 32-битной версии, и почти на 2 порядка относительно 64-битной версии. Однако, самым забавным оказалось то, что код с переменной real(10) не только точнее, но и быстрее кода с переменной real(8)! Кажется естественным, что при переходе от 64-битным к 80-битным переменным не должно быть лишних преобразований переменных при общении с мат. сопроцессором. Очевидно, это занимает какое-то время, и, видимо, в этом причина ускорения. Однако, неожиданным оказалось ускорение работы программы почти в 2 раза для Intel Xeon на cluster!

Такое наблюдение заставляет почесать себе затылок и попытаться воспроизвести эффект на других машинах и заодно проверить не является ли это глюком компиляторного происхождения или программы измерения времени. Мои усилия суммированы в следующих графиках, которые показывают относительное изменение производительности при изменении разрядности переменных 32 -> 80 и 64 -> 80. Под производительностью понимается обратное время вычисления S. Время измерялось с помощью функции omp_get_wtime() из OpenMP. Почти все машины были с 64-битной версией Linux, кроме nuke, на котором установлен 32-битный Debian (список компьютеров и их основных/существенных параметров см. в таблице ниже). Все результаты получены с флагами компиляции gfortran -fopenmp -O2 (см. также скрипт, который использовался для сбора данных).

80 vs 32 80 vs 64

___________________________________________________________________________________________
System      CPU                                 N.cores    Linux          kernel    gfortran
___________________________________________________________________________________________
cluster     Intel Xeon 3.0GHz                      4      SL x86_64      2.6.32-*   4.5.2
cluster15   AMD Opteron 270                        4      SL x86_64      2.6.32-*   4.4.7
cluster23   AMD Opteron 8214                       8      SL x86_64      2.6.18-*   4.1.2
cluster35   Intel Xeon E5530 2.4GHz                8      SL x86_64      2.6.18-*   4.1.2
cluster42   AMD Opteron 6168 1.9GHz                24     SL x86_64      2.6.32-*   4.4.6
nuke        Intel Dual-Core Pentium E5400 2.7GHz   2      Debian i686    3.2.0-4    4.7.2
rut         Intel Core i5-2410M 2.30GHz            4      Debian x86_64  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-битных процессоров Intel. Интересно также, что 80-битный код работает по крайней мере не медленнее как 64-битного так и 32-битного для всех процессоров, а для cluster и chromebook быстрее. Старик cluster вообще оказался чемпионом в этом тесте, он разгоняется почти в 2 раза в 80-битной версии. Для других тестируемых процессоров Intel наблюдается ускорение порядка 25-30% при переходе 64 -> 80. Для AMD (и 32-битной системы Linux на nuke) не наблюдается существенного выигрыша в производительности от перехода к 80-битной переменной (но точность вычислений повышается).

Результаты измерения времени оказались относительно стабильны для cluster и chromebook и достаточно сильно флуктуируют для других систем. Последовательные запуски одного и того же кода показывают, что флуктуации носят случайный характер для относительно небольших N. Возможно, это связано с конкуренцией других процессов, в том числе системных. При запуске кода на кластере, я выбирал машины свободные от пользовательских процессов (или машины с активными процессами с большим значением параметра nice, как в случае с cluster42. Я напомню, что чем больше значение nice, тем легче у процесса позаимствовать ресурсы.). Также есть подозрение, что время измеряется не очень точно в gfortran 4.1.2.

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

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

CPU performance

В этом контексте, было бы интересно протестировать вычислительные возможности процессоров ARM, которые стоят в большинстве планшетов и телефонов. Может быть уже пора попытаться использовать телефоны не только по их прямому назначению :-?)

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

Очевидно, что производительность кода зависит не только от процессора, но и от компилятора и от параметров компилирования. Графики нарисованы для результатов компиляции с флагом оптимизации O2 и могут выглядеть по-другому с другими флагами gfortran (я не изучал другие компиляторы). Флагов у компилятора много и изучение их возможностей на предмет оптимизации может легко занять все имеющееся время. Я только отмечу флаг -march, который учитывает архитектуру процессора. Желающие могут попробовать оптимизировать код под конкретную машину с помощью параметра -march=native (не работает для v.4.1.2, т.е. для cluster35 etc). Показательным в этом смысле оказался пример машины nuke, на которой установлен 32-битный Debian. Машина nuke была самой медленной в тестах с оптимизацией O2, хотя частоты процессоров у nuke и cluster сравнимы. Первым впечатлением было, что здесь сказывается разрядность ОС. Однако, использование флага -march=native оказалось весьма эффективным и производительность моего тестового кода на 32-битном ядре Linux (процессор у nuke 64-битный) оказалась сравнимой с производительностью "полностью" 64-битного cluster. С другой стороны, включение -march=native для chromebook и rut не дает дополнительного выигрыша в производительности.

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

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

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

OpenMP performance

Эффект параллельности впечатляет, производительность при достаточно больших N увеличивается пропорционально числу задействованных ядер (исключение составили 4х ядерные системы)! Параллельное выполнение суммирования в цикле также благоприятно сказывается на точности сокращения членов разного знака при больших N. Вот результаты для S(N) для 64-битной и 80-битной версий распараллеленых с помощью OpenMP для разного количества процессоров:

64-bit version:
p        N=2**p   S  No.cores=1   S  No.cores=2   S  No.cores=4   S  No.cores=8   S No.cores=24
26     67108864   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000
27    134217728  -11585.0000000  -5792.00000000   0.00000000000   0.00000000000   0.00000000000
28    268435456  -232471923.000  -183344816.000  -88438552.0000  -42486318.0000   1.00000000000
29    536870912  -519822953.000  -461237984.000  -427413984.000  -310103436.000  -137253444.000
30   1073741824  -1065319903.00  -1002421120.00  -985636384.000  -940626304.000  -738108400.000
31   2147483648  -2143284860.00  -2078265344.00  -2069877248.00  -2047507168.00  -1923873512.00
32   4294967296  -4292867744.00  -4226779136.00  -4222556160.00  -4211342720.00  -4150743808.00
33   8589934592  -8588876000.00  -8519720960.00  -8518819840.00  -8513936384.00  -8484304896.00
34  17179869184  -17179339772.0  -17062330368.0  -17075896320.0  -17095215104.0  -17092972544.0
80-bit version:
p        N=2**p   S  No.cores=1   S  No.cores=2   S  No.cores=4   S  No.cores=8   S No.cores=24
26     67108864   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
27    134217728   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
28    268435456   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
29    536870912   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
30   1073741824   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
31   2147483648   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
32   4294967296   0.00000000000   0.00000000000   0.00000000000   0.00000000000   0.00000000000 
33   8589934592  -6074000997.00  -4294967292.00  -1464933360.00   0.00000000000   0.00000000000 
34  17179869184  -16070296103.0  -13588082112.0  -11364843456.0  -7094119694.00  -591823493.000 

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

Сравнение времен затраченных на вычисление S приводит к картине похожей на результаты на 1 процессоре. Процессоры Intel(64-bit) работают быстрее с 80-битной версией кода (для Intel Xeon сервера cluster имеем выигрыш в 2 раза), тогда как производительность AMD довольно безразлична к разрядности S. Но в любом случае точность вычислений существенно повышается.

Summary

Комментарий: исходные коды и скрипты можно найти здесь.

© С.А. Кулагин, ИЯИ ОТФ, sergey.kulagin на gmail.com