Transcript Реализация базовых функций задачи горения

Реализация базовых функций
задачи горения на основе операции
FMA специализированного
векторного сопроцессора
Ивасюк Евгений Вячеславович
Научно-исследовательский институт системных
исследований РАН
Результаты профилирования
Метод ETDRK4 для системы ОДУ, кинетика Мааса и
Варнаца (горение водорода, 9 компонентов, 19 реакций),
газодинамический код с учетом той же кинетики, решение
системы проводится 4-х стадийным методом Роземброка
Время
выполнения
Название
функции
20.93
GetTransport
11.67
__ieee754_exp
6.88
UpdateExrate
6.70
ElemCofs
5.11
4.02
3.53
3.47
2.59
2.53
GetThermo
Примечание
Коэффициенты переноса для газовой смеси на основании молярных
масс или плотностей компонентов
Библиотечная функция вычисления вещественной экспоненты
(библиотека libm-2.17.so)
Добавка данных от элементарной реакции в интенсивность образования
компонент
Элементарный механизм: коэффициенты для прямой и обратной
реакции
ExplTurbStep
Явная стадия турбулентности
GetJacobean
Вычисление якобиана
__pow_finite
Библиотечная функция возведения в степень (библиотека libm-2.17.so)
GetDifComp
Учет диффузии компонентов
Y2X
Вычисление плотности компонентов в каждой ячейке
Направления оптимизации
задачи горения
•
•
Ускорение выполнения матрично-векторных
операций (умножение, сложение):
 Оптимизация библиотеки BLAS
Ускорение вычисления скалярных специальных
функций (экспоненты, деление, кв. корень,
дробная степень):
 Аппаратная поддержка вычисления
функций
 Программно-аппаратная поддержка
вычисления функций
Структура базового вычислительного
элемента векторного сопроцессора
Мантисса А
Мантисса В
Мантисса С
Порядок В
Порядок С
Вычисление порядков
(Exponent calculation)
Кодирование Бута
(Booth Encoding)
Выравнивание
(Alignment shift)
...
Порядок А
Дерево Уоллеса
(Wallace tree)
Дизъюнкция битов выдвинутых
за пределы разрядной сетки
Сумматор
Сумматор3>2
3>2(Сarry
(Сarrysave
saveadder)
adder)
Предсказатель
Предсказатель
суммы
суммы(LZA)
(LZA)
Счетчик
Счетчикнулей
нулей
(LZC)
(LZC)
Детектор
Знака
результата
(sign detector)
Предварительный
Предварительныйсумматор
сумматор(Pre
(PreAdder)
Adder)
Нормализация
Нормализация(Norm
(Normshift)
shift)
Мультиплексор (MUX)
Дизъюнкция
отбрасываемых битов
Окончательный сумматор, округление
(Rounding Adder)
окончательная
нормализация
(Post Normalization)
Мантисса результата
Коррекция порядка результата
(Result exponent correction)
Порядок результата
Структурная схема
вычислительного ядра
Операция наибольшей производительности:
Умножение со сложением и вычитанием
комплексных чисел («бабочка Фурье»)
CMADDSUB
(A + iB) * (C + iD) + (E + iF) = (G + iH)
(A + iB) * (C + iD) - (E + iF) = (K + iL)
Операция матрица-вектор:
Умножение вещественной матрицы на
вещественный вектор со сложением
MVMADD
K
G
=
L
A B
+
H
E
*
C D
F
K = G + A * E + B * F
L = H + C * E + D * F
Выполнение функций GEMM
Результаты измерения эффективности выполнения функции GEMM векторным сопроцессором, %
M×N×K
L2(512К)
ZGEMM
4×4×8
8×4×8
32×4×8
256×4×8
4×8×8
8×8×8
32×8×8
256×8×8
4×8×16
8×8×16
32×8×16
256×8×16
16×16×16
1024×4×8
1024×8×8
4096×8×8
256×32×32
64×64×64
256×64×64
1024×8×16
4096×8×16
1024×16×16
0.34%
0.59%
2.05%
15.72%
0.6%
1.0%
3.32%
25.2%
0.88%
1.37%
4.3%
31.64%
3.9%
62.6%
100.2%
400.2%
103.12%
62,5%
212.5%
125.39%
500.4%
200.8%
23.92%
34.12%
52.24%
60.72%
30.93%
40.59%
55.05%
61.08%
35.88%
45.13%
57.95%
62.79%
55.30%
61.73%
61.83%
37.48%
63.92%
63.11%
62.15%
CGEMM
DGEMM
SGEMM
21.18%
31.07%
47.19%
55.52%
25.83%
35.65%
49.91%
57.15%
47.44%
14.21%
22.38%
42.16%
57.72%
17.69%
26.99%
46.34%
58.76%
42.24%
11.37%
18.73%
38.04%
55.12%
33.33%
56.56%
60.20%
42.41%
61.78%
59.02%
61.96%
61.66%
38.22%
61.86%
57.95%
56.61%
58.11%
58.06%
58.15%
58.47%
54.73%
59.53%
57.98%
58.29%
Структурная схема аппаратного
модуля вычисления экспоненты
Приведение аргумента
Формула приведения:
r  x  k C
k  INT _ RND( x / C )
C C
r  [
, ]
2 2
C  2 m ln(2)
Внутренняя разрядность представления данных:
WRED  N ZER _ WORST _ CASE  WMAN
N ZER _ WORST _ CASE
W MAN
- число старших нулей при
вычитании
- разрядность мантиссы
Наихудший случай приведения
W MAN  53 C  ln(2) / 4096
x _ worst _ case  112.391547
70762435
Обход проблемы наихудшего
случая
Увеличение разрядности константы С:
C  C1  C 2  2n1
Формула приведения:
r  x  k  C1  k  C 2  FMS( FMS( x, k , C1), k , C 2)
 FMSmod ( x, k , C1, C 2)
FMS – multiply-sub fused, умножение с вычитанием
...
Дерево Уоллеса 2
(Wallace tree)
Мантисса С
Мультиплексоры
...
Выравнивание
(Alignment shift)
Дерево Уоллеса 1
(Wallace tree)
Выравнивание
(Alignment shift)
Предсказатель
Предсказатель
суммы
суммы(LZA)
(LZA)
Детектор знака результата
(sign detector)
Сумматор
Сумматор5>2
3>2(Сarry
(Сarrysave
saveadder)
adder)
Счетчик
Счетчикнулей
нулей
(LZC)
(LZC)
Предварительный
Предварительныйсумматор
сумматор(Pre
(PreAdder)
Adder)
Нормализация
Нормализация(Norm
(Normshift)
shift)
Мультиплексор (MUX)
Окончательный сумматор, округление
(Rounding Adder)
окончательная
нормализация
(Post Normalization)
Коррекция порядка результата
(Result exponent correction)
Мантисса результата
Порядок результата
Порядок С
Мультиплексоры
Мантисса В1
Порядок В1
Мантисса А
Порядок А
Мантисса В2
Кодирование Бута
Порядок В2
Вычисление
сдвига
Порядок В1
Структурная схема модуля FMA c
модернизацией
Вычисление
порядков
(Exponent
calculation)
Вычисление полинома
1. Прямое вычисление полинома по схеме Горнера
P( x)2n  ((...(( a2n  x  a2n 1 )  x  a2n 2 )...)  x  a0 ).
2. Снижение степени разложением на множители
P( x)2n  ((...(( an  x  an1)  x  an2 )...)  x  a0 ) 
xn  ((...(( a2n  x  a2n1)  x  a2n2 )...)  x  an1).
3. Рациональная аппроксимация
P( x) n 
((...(( an  x  an  1)  x  an  2 )...)  x  a0 )
((...((bn  x  bn  1)  x  bn  2 )...)  x  b0 )
.
Структурные схемы потоковых
вычислителей полинома
IN
ARG_RED
...
...
FMA_N
1FMA_N/2
2FMA_N/2
RECONSTR
2FMA_1
1FMA_1
...
2FMA_1
1FMA_1
...
FMA_1
ARG_RED
...
ARG_RED
IN
1FMA_N/2
2FMA_N/2
POWERING
IN
DIV
FMA
RND
RECONSTR
RECONSTR
OUT
RND
OUT
RND
OUT
a
b
c
а : Ttotal=Tred+Tfma*N+Treconstr+Trnd
b : Ttotal=Tred+Tfma*N/2+Tdiv+Treconstr+Trnd
c : Ttotal=Tred+Tfma*N/2+Tfma+Treconstr+Trnd
Зависимость степени полинома от
константы приведения
Зависимость задержки вычисления
полинома от константы приведения
Применение компенсированных
операций
Сложение
function
[ x, y ]  TwoSum(a, b),
x  fl(a  b), y  err
| a || b |: x  fl(a  b),

 y  fl( fl(a  x)  b)
| a || b |: x  fl(a  b),


 y  fl( fl(b  x)  a)
Умножение
function
[ x, y]  TwoMul(a, b),
x  fl(a * b),
y  FMS(a, b, x)  x  a * b
Обычная схема Горнера
r  Horner( p, x)
0
r a
n
n
i  n  1... 0
ri  ri 1  x  ai
P( x)  ((a n  x  a n1 )  x  ...)  x  a0
Компенсированная схема Горнера
r  CompHorner( p, x)
r  a , c  0, i  n  1... 0
n
i n
[ p i ,  i ]  TwoMul(ri 1 , x)
[ri ,  i ]  TwoSum( p i 1 , a i )
c i  c i 1  x   i   i
r  Horner( p, x), c  Horner( p  p , x)
0
0


r  r c
0
0
x
ADD
MUL
ai
Horner(p,x)
  err1 ,  err2
i
i
i
i
Ci
ADD
ri 1
Ci 1
SUB
FPU
err1i
ADD
FMS
векторный
сопроцессор
FMA
Структурная схема вычислителя
полинома компенсированного
типа
err2i
ri
Результаты






Проведена оптимизация библиотеки BLAS под перспективный
векторный сопроцессор
Проведена оценка эффективности выполнения функций библиотеки
(достигает 64%)
Предложена модернизация модуля FMA для обеспечения
эффективного выполнения стадии приведения аргумента
Предложены варианты организации потоковых вычислителей
полиномов по схеме Горнера
Предложен вариант компенсированного вычислителя полинома,
использующий основной блок вещественной арифметики и
сопроцессор
Работа выполнена в рамках проекта РФФИ № 13-07-12062
«Фундаментальные проблемы создания микропроцессоров и
коммуникационных сред супер-ЭВМ эксафлопсного класса,
ориентированных на предсказательное моделирование задачи
горения».