Рост производительности современных компьютеров позволяет решать все более и более сложные задачи, что приводит к необходимости выполнять большое количество операций с плавающей точкой, почти каждая из которых сопровождается ошибкой округления. Экзафлопсный уровень скорости вычислений (10 18 операций в секунду), вероятно, будет достигнут в ближайшие десять лет и получить точные результаты на таких вычислительных системах в машинной арифметике IEEE-754 будет непростой задачей [1].

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

1. Плохо обусловленные линейные системы. Многие проблемы нехватки точности связаны с необходимостью решения плохо обусловленных систем большого порядка. Например, в задачах строительной механики возникают системы с десятками миллионов неизвестных и с числом обусловленности 10 9 -10 13 [2]. Также задачи данного класса возникают при изучении кулоновских атомных систем [3] и при исследовании электромагнитного рассеяния [4,5]. В задаче среднеквадратичной аппроксимации непериодических функций неортогональным степенным базисом возникает необходимость решения систем с матрицей Гильберта, которая известна плохой обусловленностью [6].

2. Рекуррентные формулы и большие суммы. Аномальные результаты часто связаны с потерей ассоциативности при суммировании, особенно при его выполнении на параллельной компьютерной системе, где порядок суммирования не может контролироваться [7,8]. Необходимость вычисления больших сумм возникает, в частности, при использовании квадратурных формул для интегрирования полиномов высоких степеней на сетках с малым шагом. К практическим задачам, основанным на вычислении рекуррентных соотношений, чувствительных к ошибкам округления, относится реализация дискретного преобразования Фурье в форме рекурсивного фильтра -алгоритм Гёрцеля [9], а также рекурсивное дискретное косинус-преобразование [10]. Эти задачи имеют большое значение в цифровой обработке сигналов. В теории приближений \ast Работа выполнена при финансовой поддержке РФФИ в рамках научного проекта № 14-07-31075 мол_а.

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org важную роль играют полиномы Чебышева, которые также вычисляются с помощью рекуррентных соотношений.

3. Жесткие системы дифференциальных уравнений. Такие системы возникают в задачах химической кинетики с одновременным присутствием очень медленно и очень быстро протекающих химических реакций [11], при исследовании суточных колебаний озона в атмосфере [12], а также в задачах расчета динамики многозвенных систем [11].

4. Крупномасштабное моделирование. Вычисления, устойчивые на небольших задачах, выполняемые на однопроцессорных системах, могут сопровождаться значительными численными ошибками при масштабировании до уровня массивно-параллельных систем. К крупномасштабным задачам относится моделирование климата [4,7,13,14], моделирование атмосферы сверхновых звезд [7], моделирование процессов, протекающих в ядерных реакторах, расчет дозвукового обтекания летательного аппарата [13].

5. Продолжительное моделирование. Практически любой процесс физического моделирования, выполняемого в течение длительного времени, в конечном итоге отходит от реальности. Это связано с накоплением ошибок округления в дополнение к ошибкам, связанным с дискретизацией по времени и пространству. Показательным примером является задача численного исследования орбитальной эволюции небесных тел [15,16].

6. Экспериментальные математические расчеты. Многие результаты в этой сравнительно молодой области научных исследований (например, формула Бэйли -Боруэйна -Плаффа для вычисления числа \pi ) не могли быть получены без использования вычислений с очень высокой точностью [7,16].

Приведенные области применения высокоточных вычислений характеризуются, как правило, большой размерностью задач. Вместе с тем, даже при малом объеме вычислений в машинной арифметике может быть получен результат, не содержащий ни одной значащей цифры [17][18][19]. Таким образом, проблема обеспечения приемлемой точности вычислений актуальна в настоящее время и с ростом объемов производимых вычислений, ее значение, несомненно, будет возрастать.

Одним из способов повышения точности в настоящее время является 128-битный IEEE формат [20], в котором поле мантиссы расширено до 113 разрядов. Однако аппаратная поддержка этого формата требует значительных затрат [16] и, судя по всему, не предвидится в ближайшей перспективе. Более распространенным вариантом является длинная арифметика -программная обработка чисел, разрядность которых превышает стандартную длину машинного слова вычислительной машины. В настоящее время существует достаточно широкий спектр библиотек длинной арифметики. К наиболее известным относятся следующие. Ознакомиться с представленными критериями эффективности высокоточного программного обеспечения, а также найти другие критерии, можно найти в работе [16], которая и взята за основу изложенного материала.

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Следует отметить, что в настоящее время значительно возрастает роль векторных (SIMD) вычислений, в том числе и на центральных процессорах, длина векторных регистров которых неизменно увеличивается. В связи с этим важным требованием является возможность распараллеливания высокоточных арифметических операций на уровне отдельных цифр многоразрядных мантисс. Но такая задача весьма трудоемка в рамках модели позиционной длинной арифметики, которая предполагает в процессе вычислений учет возможных переносов между соседними блоками цифр многоразрядной мантиссы. В результате этого алгоритмы обработки длинных чисел сильно ветвятся и не распараллеливаются, что противоречит основной концепции SIMD-вычислений. Распараллеливание на уровне отдельных цифр мантиссы наиболее актуально при реализации операций с очень высокой точностью, которые выполняются во много раз медленнее операций машинной точности.

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

3. MF-Library -программная библиотека параллельных высокоточных вычислений

В библиотеке MF-Library для внутреннего представления длинных чисел с плавающей точкой используется модулярно-позиционный формат (MF-формат) [27,28], описание полей которого представлено в таблице 1. В MF-формате для представления длинной мантиссы M используется система счисления с параллельной структурой -система остаточных классов (СОК) [29,30]. В соответствии с этим, мантисса представляется набором остатков residue 1 , residue 2 , . . . , residue n по модулям СОК p 1 , p 2 , . . . , p n , где residue i \equiv M mod p i . Обработка residue i по различным модулям p i выполняется параллельно. Мантисса может принимать значения из диапазона [0, P ), где P -произведение всех p i . Следовательно, варьирование количества модулей позволяет задавать произвольную точность вычислений. Обеспечить достаточно большой динамический диапазон, что важно в научных вычислениях, где часто требуется обработка величин очень большого или очень малого масштаба, позволяет двоичный порядок (exp). Значение числа в MF-формате определяется следующим выражением:

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org MF-формат обладает следующими основными характеристиками точности: машинный эпсилон -\epsilon = 2 - \lfloor \mathrm{ \mathrm{ \mathrm{ 2 (P - 1)\rfloor , unit in the last placeulp(x) = 2 \mathrm{ \mathrm{ \mathrm{ - \lfloor \mathrm{ \mathrm{ \mathrm{ 2 ((P - 1)/M )\rfloor , где M -десятичное значение мантиссы. Абсолютные и относительные ошибки округления ограничены, соответственно, величинами ulp(x) и 2\epsilon при округлении мантиссы усечением [27]. Точность вычислений (в терминах IEEE-754 точность соответствует разрядности мантиссы) определяется половиной длины диапазона изменения мантиссы, т.е. величиной log 2 \surd P -1. Дополнительно в MF-формат включена атрибутивная информация -интервальнопозиционная характеристика (ИПХ) мантиссы, представленная двумя направленно округленными двоичными числами с плавающей точкой -нижней (ic_bot) и верхней (ic_top) границами. ИПХ не участвует в образовании значения числа, но позволяет в значительной степени преодолеть основной недостаток СОК -высокую сложность немодульных операций, таких как сравнение, определение знака, контроль переполнения, масштабирование и пр. ИПХ локализует отношение величины мантиссы M к произведению модулей P так, что ic_bot \leq M/P \leq ic_top. Вопросы использования ИПХ в немодульных вычислениях обсуждаются в [27,31].

Структурно MF-Library включает в себя три слабосвязанных уровня, каждый из которых объединяет однотипные модули (Рис. 1). Связь уровней реализуется посредством API, который также доступен пользователю.

Центральной частью пакета является модуль ядра. В подмодуле MF-формата объявлен основной тип данных, реализованы функции по выделению и освобождению памяти, а также подпрограммы инициализации основных статичных объектов, используемых другими подмодулями. Арифметический подмодуль реализует следующие алгоритмы высокоточных вычислений: сложение, вычитание, умножение, деление, выравнивание порядков, сравнение. Подробное описание этих алгоритмов с оценками эффективности можно найти в [27,28]. Подмодуль IPC реализует алгоритмы вычисления и анализа интервально-позиционных характеристик, необходимые для быстрого выполнения немодульных операций над мантиссами. Подмодуль округления включает в себя подпрограммы модулярного масштабирова-ния степенью двойки, проверки необходимости округления и округления чисел. Подмодуль ввода / вывода обеспечивает следующие возможности: ручной ввод MF-числа, установка MF-числа из форматов double и mpfr_t (многоразрядный тип данных библиотеки MPFR), установка случайного MF-числа, преобразование MF-числа в тип double (с возможной потерей точности), форматированная печать. В подмодуле служебных констант объявлены идентификаторы статичных объектов, вычисляемых при инициализации пакета. В подмодуле основных параметров содержатся управляющие флаги и константы, определяющие все аспекты работы пакета. К их числу относятся: точность вычислений, модули СОК, флаг векторизации, флаг режима отладки, флаг работы в "тихом" режиме, в котором игнорируются некоторые некритичные исключения арифметики с плавающей точкой.

модуле служебных подпрограмм содержатся утилиты для управления режимами округления, генераторы случайных чисел, функции преобразования типов, базовые алгоритмы модулярной арифметики: получение мультипликативной и аддитивной инверсии, алгоритм Евклида, преобразование в позиционную систему, алгоритм генерации оснований для заданной точности вычислений и пр. Для отладки ядра в состав пакета включен модуль низкоуровневых тестов. Для автоматической конфигурации MF-Library и перехода на другой базис СОК (при изменении точности вычислений) используется вспомогательный модуль длинной позиционной арифметики (MPFR+GMP).

Промежуточный уровень реализует итерационные методы высокоточного вычисления математических констант и некоторых классов элементарных функций.

На прикладном уровне реализованы высокоточные матричные, матрично-векторные и векторные операции, входящие в состав BLAS (GEMM, GEMV и пр.), метод сопряженных градиентов для решения больших систем уравнений и конечно-разностный метод решения краевой задачи теплопроводности.

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

Для экспериментальной оценки корректности и эффективности методов каждого уровня пакета MF-Library были разработаны и выполнены соответствующие тесты. Остановимся на рассмотрении результатов некоторых из них 1 .

Расчеты с катастрофической потерей точности. Исследовалась корректность MF-Library при решении задач, которые даже при незначительном объеме вычислений сопровождаются возникновением катастрофических погрешностей. Такие задачи часто используются для верификации высокоточных программных средств. Первая задача заключалась в вычислении следующего полинома восьмой степени (S.M. Rump, 1988 [17]):

при a = 77617.0 и b = 33096.0. В качестве эталонного использовалось решение, полученное с использованием 4096-битной арифметики библиотеки MPFR. Результаты представлены в таблице 2. Более высокая точность MF-Library по сравнению с MPFR (256 бит) объясняется использованием предвычислительной схемы округления [28], которая позволяет использовать для представления мантисс весь динамический диапазон, а не его половину. Вторая исследованная задача -вычисление рекуррентного соотношения Мюллера [18]:

x 0 = 4.00, x 1 = 4.25,

Параметры этого соотношения подобраны таким образом, что при точных вычислениях lim n\rightar\infty x n = 5.0. Однако из-за влияния ошибок округления последовательность \{ x n \} отходит от верного ответа к неподвижной для данного рекуррентного соотношения точке x \ast = 100.0. Номер итерации, с которой начинается такой переход, прямо пропорционален точности вычислений. Результаты эксперимента представлены на Рис. 2. Эффективность векторизации. Исследованы операции сложения (add), вычитания (sub), умножения (mult), сравнения (cmp), сложения с накоплением (aac, x \leftarr x + y), вычитания с накоплением (sac, x \leftarr x -y), умножения с накоплением ( Эффективность векторизации оценивалась по формуле В работах [27,28] представлены результаты других экспериментов по исследованию эффективности библиотеки MF-Library.

С ростом производительности вычислительных систем и размерности решаемых задач возрастает значимость вычислений многократной точности. Уже сейчас имеется целый ряд актуальных приложений, которые требуют оперировать числами с разрядностью, на несколько порядков превосходящей длину машинного слова. Следует полагать, что с течением времени количество таких приложений будет увеличиваться. В связи с этим, проблема обеспечения приемлемой точности отмечается многими исследователями, как одна из наиболее актуальных проблем программного обеспечения экзафлопсных суперкомпьютеров.

Современный этап развития вычислительной техники приводит к новым требованиям, предъявляемым к высокоточному программному обеспечению. К числу основных из таких требований относятся высокая скорость, потоковая безопасность и эффективное использование современных параллельных архитектур. В этой статье рассмотрен новый программный пакет высокоточных вычислений MF-Library, нацеленный на удовлетворение этих требований. В основе пакета лежит модулярно-позиционный способ представления длинных чисел, изначально ориентированный на параллельную обработку. Алгоритмы, заложенные в MF-Library, лишены недостатков аналогов, связанных с необходимостью затратной обработки переносов между соседними цифрами длинной мантиссы, и соответствуют архитектуре современных центральных процессоров, позволяя эффективно использовать как многоядерность, так и возможности SIMD-обработки в пределах ядра. Благодаря потоковой безопасности, MF-Library может эффективно использоваться на параллельных системах с общей памятью, что подтверждается представленными результатами экспериментов.

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org В настоящий момент ведутся работы по расширению функциональности MF-Library. Одновременно с этим начата работа по созданию версии пакета для графических процессоров. Конечной целью исследований является создание унифицированного алгоритмического и программного обеспечения, позволяющего эффективно выполнять высокоскоростные параллельные расчеты с плавающей точкой многократной точности на современных и перспективных высокопроизводительных системах с гибридной архитектурой.