О реализации QR-разложения на трехмерном систолическом массиве

Непрерывно возрастающие объемы вычислений обуславливают необходимость повышения производительности вычислительных систем. Одно из направлений в этой деятельности — разработка и применение систолических массивов (СМ). Вопросам разработки СМ различного назначения посвящено немало работ, так как они являются одним из наиболее эффективных средств решения указанной проблемы [1 –8] .

Одной из важнейших задач вычислительной линейной алгебры является решение систем линейных алгебраических уравнений (СЛАУ). Пусть матрица системы

Ax = f невырождена, и N — ее порядок. Для осуществления QR-разложения матрицы A мы будем использовать преобразования (матрицы) вращения.

Введем обозначения A ^(0,0) = A, Q * j — преобразование вращения, аннулирующее i-ю компоненту j -ого столбца матрицы.

Для j = 1, N — 1, i = j, N вычисляем:

(j-1,i-1)               xj

• x = aj , c =   ;          , s =   ;,

• У xj + ^x 2

(j — 1,i — 1)

для k = j,N : x = a k     ⁷, y = Sj^x, a ^ = y.

По завершении этой процедуры получим верхнюю треугольную матрицу R = A ^{( N —1 ,N} ) , связанную с исходной соотношением B = S * A, где

S = S 1 , 2 ...S 1 ,N S l , 3 . . . S l ,N . . . S n — 2 ,N — 1 S n — 2 ,N S n — 1 ,N .

Следует сказать, что в литературе для представления матрицы А в виде произведения A = QR, где Q и R — ортогональная и верхняя треугольная матрицы соответственно, широко применяется термин QR-разложение. Следуя этой традиции, произведение Q A мы будем называть QR-упрощением матрицы А .

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

Второй областью применения QR-разложения матрицы, очевидно, являются задачи решения СЛАУ большой размерности с разреженными матрицами. Для решения таких систем применяются итерационные методы. Пусть x0 — произвольный вектор, r0 = f — Ax0 и xn+1 = xn + yn, n = 0, 1, 2, ... .

В проекционных методах вычисление поправки y ⁿ сводится к последовательному решению трех систем линейных уравнений

Vnwn = rn,  Hn zn = wn,  V* yn = zn, где Vn — псевдоортогональная матрица, rn = f — Axn и Hn — некоторая плотная невырожденная матрица порядка m, причем m невелико.

Из трех выше представленных систем наиболее трудоемким является решение второй системы. Системы линейных уравнений H n z n = w n , n =0,1, 2,..., обозначенного потока могут быть решены с помощью QR-разложения матрицы H _n , реализованного аппаратурно на специально разработанном вычислительном устройстве (систолическая матрица Джентльмена–Кунга [9] ).

Обоснование заявленного утверждения упрощенно раскроем сле- дующим образом. Для реализации преобразования вращения вектора разработано специальное вычислительное устройство CORDIC, которое производится в двух вариантах: R1, аннулирующее вторую компоненту входного вектора, и R2, работающее в режиме «эхо».

На вход R1 поступают компоненты x и у двумерного вектора

x

.

zy

На выходе получают вектор , где z = sign(x) у/ x ² + y² , и параметры преобразования вращения c и s.

На вход R2 поступают компоненты и и v двумерного вектора

uv

и

параметры преобразования вращения c и s. Соответственного на выхо-

де получают вектор вращения c и s.

[u] = [—s 3 [u]

и параметры преобразования

Устройство CORDIC осуществляет вычисление компонент вектора за 2m тактов, где m — количество разрядов, отводимых под мантиссу машинных чисел. Устройства CORDIC получили применение в систолической структуре (матрице Джентльмена–Кунга), применяемой для приведения матрицы к верхнему треугольному виду. Эти устройства в конвейерном варианте способны обрабатывать потоки данных высокой плотности, но размерность систолической структуры Джентльмена–Кунга (она равна 2) препятствует реализации этой способности указанного устройства. Выходом из этого положения может служить переход к применению трехмерных систолических структур.

Следует сказать, что недостатком метода CORDIC является не слишком высокая скорость сходимости (для достижения приемлемой точности требуется 2m итераций). Соответственно этому задержка вычисленного результата на выходе устройства CORDIC составит 2m тактов. Одним из авторов этой статьи разработаны модификация метода CORDIC, позволившая сократить число итераций до m, и соответствующее ей устройство, обеспечивающее возможность получить вычисленный результат за m тактов [10–13]).

Выполняя эту работу, мы преследовали цель обеспечить обработку потоков данных высокой плотности при осуществлении решения СЛАУ. Достижение этой цели мы видели в решении задачи разработки систолических структур, использующих в качестве функциональных устройств (ФУ) CORDIC-устройства конвейерного типа или их указанную выше модификацию [10 –12] .

1.    Основные положения

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

Пусть T m и T m +1 — моменты поступления на ФЭ m -го и m +1-го тактовых импульсов. Промежуток времени A T = T m +1 — T m называют m -м тактом.

Символами I ФЭ и O ФЭ обозначают вход и выход ФЭ. Обозначение 1 фэ (т) = а означает, что на вход ФЭ на m -ом такте (обычно в начале промежутка A T ) поступает сигнал (число) а.

Под воздействием сигнала a и m -го тактового импульса ФЭ с течением времени A t переходит в новое устойчивое состояние, которое наблюдают на его выходе. При этом пишут О фэ (т) = b, где b — сигнал (число) на выходе ФЭ на m -ом такте, соответственно At называют временем срабатывания ФЭ. Выполнение неравенства At <  A T называют условием срабатывания ФЭ.

Переходя к формальному изложению, последнее мы можем записать так.

Лемма 1.1 (о срабатывании ФЭ). Пусть 1 фэ ( т ) = а. Если A t <  A T, то О Фэ ( т ) = b.

Выражение 1 Фэ 2 <= О Фэ 1 (О ФЭ 1 => 1 ФЭ 2 ) означает, что выход ФЭ 1 соединен с входом ФЭ ₂ .

Определение 1.1. Если на m-том такте на выходе ФЭ появляется результат b ( т.е. O фэ ( m ) = b ) , то его мы будем называть тактом, порождающим результат b при отображении ФЭ ( обозначение m = = ind фэ (b)).

Лемма 1.2 (о локальном взаимодействии ФЭ). Пусть 1 Фэ 2 < = О Фэ 1 . Если 1 фэ 1 ( m ) = а и О ФЭ 1 ( m ) = b, то 1 ФЭ 2 ( т + 1) = b.

Определение 1.2. Конвейером C мы будем называть функциональное устройство, представляющее собой последовательность ФЭ

C = {F h ,..., F i } : i F p +i <= O F p , p = 1, h - 1.

Соответственно I c = I F 1 — вход, O c = O F h — выход, и h — длина конвейера C.

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

Лемма 1.3 (о необходимом условии функционирования конвейера). Пусть конвейерный вычислитель С с регулярн ым графом состоит из n функциональных элементов, причем A t i , i = 1 ,h — время срабатывания i-го элемента, и справедливы соотношения: а 1 = F 1 ( а ) , а ₂ = F 2 ( a 1 ) , .. ., a h - 1 = F h - i ( a h - 2 ) , b = F h ( a h - 1 ) . Если I c ( m ) = а, и A t <  A T, где A t >  max { A t 1 ,..., A t h }, то O c ( m + h — 1) = a h = b.

Действуя по индукции, очевидно, применяя при этом леммы 1.1 и 1.2, мы докажем утверждение последней леммы.

Если l R ( m ) = а и O r ( m ) = а, то функциональный элемент R называют элементом задержки, а конвейер LD , представляющий собой последовательность элементов задержки R = {R h ,..., R 1 }: I r _{p +1} < = < = O r _p , p = 1 ,h — 1, I ld = I R , O ld = O R h , называют линией задержки, соответственно I _LD – входом, O _LD – выходом и h – длиной линии задержки.

2.    QR-упрощение матрицы на трехмерном СМ

Поскольку назначением предлагаемого ниже систолического массива прямого хода (СМПХ) является приведение матрицы системы Ax = у к треугольному виду (этап 1) прежде всего для реализации последующего обратного хода (этап 2), и при осуществлении этапа 1 происходит трансформация вектора y, как и всех вектор-столбцов матрицы А, мы сформируем расширенную матрицу A y . Эту матрицу, как и ранее, мы будем обозначать символом А, а вектор у считать (N +1)-м столбцом расширенной матрицы А.

• 2.1.    Конфигурация СМПХ

Для реализации QR-упрощения мы применим трехмерную матрицу ФУ R, состав элементов которой, описывается соотношениями:

{ LD V j, если i = k,

R 1 V i, если j = k,

R2 в остальных случаях, где k = 1, N — 1, i = k, N, j = k, N +1.

В целях упрощения (единообразия) дальнейшего изложения мы будем придерживаться следующих соглашений:

lRij,k = 12Rj , Or^ = O1Rij,k ,k = 1, N - 1, i = k, j = k, N + 1.

Функциональные устройства связаны между собой следующим образом.

• (1)    Одноуровневые (горизонтальные) связи:

O3rp, =>I3rp       k = 1, N - 1, i = k +1, N, j = k^N, i,j,k              i,j+1 ,k p = 1, m,    (направление—вправо);

O1Rij,k =>I 1R+ijk, k = 1, N - 1, i = k + 1, N - 1, j = k, N +1,    (направление — вперед («вниз»)).

• (2)    Межуровневые связи:

О2^,к => I2r   , , k = 1, N - 1, i = k + 1, N, j = k +1, N +1,    (направление—вниз).

Входы СМПХ:

1 2 R j , i = T?X j = 1, N + 1.

Выходы СМПХ:

O 1 R _N , ₃ , _k ,   k = 1, N - 1, j = k,N +1,

O2 r.n

,N,N- 1 ,

O 2 r n

,N +1 ,N - 1 ^.

• 2.2.    Основные свойства СМПХ (маршрутизация и расписание движения данных)

Пусть требуется осуществить QR-упрощение матриц

A ¹ , ..., A l , ... .

Описание свойств СМПХ предполагает применение переменных:

• t, t — номера тактов,

• т , T — номера упрощаемых матриц (решаемых систем уравнений).

Для описания свойств (технических характеристик) СМПХ, представляющихся нам основными, определим следующие величины:

A T — продолжительность порождения упрощенной матрицы, определяющаяся по формуле A T = k e — k b + 1, где k b — номер такта, на котором на вход СМПХ поступает первый элемент упрощаемой матрицы, k e — номер такта, на котором порождается последний элемент упрощенной матрицы;

Р — производительность, определяемая как количество элементов упрощаемых матриц систем уравнений, появляющихся на выходе (выходах) СМПХ одновременно.

Высказывание 2.1. Пусть

• ( i )    продолжительность A T такта на СМПХ удовлетворяет неравенству

AT > At+, где At+ — продолжительность времени срабатывания сумматора.

• ( ii )    для всех l = 1, 2, ..., i = 1,N — 1, j = 1,N + 1 выполняются соотношения

• (1)                      1 2 R _i,_j,i ( t ) = a^ j , где

t = h ( i — 1) + l, ^T = t ^{— h (} i ^— 1) ^{— ( j —} 1), причем соотношения (1) определены для всех l,i,j : т >  0.

Тогда

• 1 )    для всех k = 1, N — 1

  O1 r

  
  

  ® = bh ,

  
  

^O2 R n,n +i ,n - i ^{( h ( N + ( N —} 2) + (l ^— 1)) = b^l N,N ,

O2rn n+in-i (h(N + (N — 2)) + l) = bl+N +1, gde t = h(N +(k — 1)) + (l — 1), ге T = t — h(N +(k — 1)) — (j — 1) + 1;

• 2 )    продолжительность процедуры QR-упрощения матрицы A T равна 2 hN + N +1 тактам;

• 3 )    производительность СМПХ Р равна ( N + 3)N/2 слов за такт.

Замечание 2.1. Значения величин пунктов 2) и 3) высказывания 2.1 получены применением леммы 1.3 .

Заключение

В представленном СМПХ за счет увеличения размерности СМ до трех достигнута максимально возможная плотность потоков движения данных. Она ограничена продолжительностью выполнения операции сложения △t+. Сопоставляя упомянутую во введении матрицу Джентльмена-Кунга с СМПХ, мы видим, что производительность разработанного СМПХ в n х 2h выше, чем у СМ Джентльмена-Кунга. В n раз выше за счет увеличения размерности массива до трех, ив 2h раз —за счет появившейся при этом возможности увеличить плотность потока данных в 2h раз, где m для устройства CORDIC, m/2 для модифицированного CORDICа, m — количество двоичных разрядов, отводимых под мантиссу числа.

Проводя обзор публикаций последних лет по реализации QR-разложения на трехмерных систолических массивах с использованием CORDIC-устройств, мы обнаружили следующее положение дел.

В работе [14] предложена архитектура двумерного систолического массива, который может быть использован в сверхбольших интегральных схемах (СБИС) для аппаратной реализации алгоритма QR-разложения комплексных матриц. Предложенная архитектура является эффективной для реализации в СБИС, так как использует только алгоритм CORDIC и, следовательно, не использует операций умножения. Рассматриваемый систолический массив основан на методе треугольных комплексных вращений и позволяет получить значительный выигрыш в производительности по сравнению с традиционно используемым методом комплексных вращений Гивенса.

В статье [15] исследуется специфическое исключение элементов для решения линейной алгебраической системы вида Az = b, где матрица A задана суммой m парных произведений тёплицевых треугольных (нижней и верхней) матриц порядка n. Предложен алгоритм, требующий лишь (3 m — 2) n ² умножений и такого же числа сложений. Для симметричной матрицы A предложен алгоритм, требующий лишь (2 m — 1) n ² умножений и такого же числа сложений. Алгоритмы обладают регулярной структурой и допускают параллельную реализацию за O( nm ) шагов.

В работах иностранных авторов представлены результаты по практической реализации QR-разложения на двумерных систолических массивах.

В статье [16] дано описание двумерного систолического массива, реализованного на ПЛИС Xilinx Virtex5 и предназначенного для реализации QR-разложения с использованием алгоритма вращения Гивенса.

В работе [17] представлена улучшенная архитектура систолического массива для реализации QR-разложения на основе метода вращения Гивенса (GR) для действительной матрицы 4 х 4. Алгоритм цифрового компьютера вращения координат (CORDIC) принят и модифицирован для ускорения и упрощения процесса GR.

В статье [18] представлен улучшенный аппаратный дизайн с фиксированной точкой QR-разложения, специально оптимизированный для ПЛИС Xilinx. Алгоритм вращения Гивенса реализован с использованием свернутого систолического массива и алгоритма CORDIC.

В статье [19] представлена параллельная архитектура систолического массива QR-разложения на основе алгоритма вращения Гивенса на ПЛИС. Предлагаемая архитектура использует прямое отображение 21 процессорных единиц с фиксированной точкой на основе CORDIC, которые могут вычислять QR-разложение для действительной матрицы 4 х 4.

По трехмерным систолическим массивам мы обнаружили только три работы

В [20] предложена архитектура трёхмерного систолического вычислителя дискретного преобразования Фурье на основе кронекеровского произведения матриц.

Авторы в статье [21] сообщают, что ими обнаружено следующее. Систолический массив как уровень программной виртуализации может привести к чрезвычайно масштабируемой парадигме выполнения. Чтобы продемонстрировать эту масштабируемость, они спроектировали и реализовали виртуальный трехмерный систолический массив для реализации QR-разложения. На машине Cray-XT5 продемонстрировано, что виртуальный систолический массив может обеспечить превосходную параллельную производительность.

В работе [22] рассмотрена реализация HLS архитектуры трехмерного систолического массива для умножения матриц, которая нацелена на определенные характеристики ПЛИС Intel Stratix 10, чтобы создавать проекты, которые достигают высокой пропускной способности с плавающей точкой.

Представленный анализ публикаций последних лет позволяет сделать вывод, что полученные в нашей работе результаты являются новыми.

Подведем итоги проделанной работы. Для приведения матрицы СЛАУ к верхнему треугольному виду разработан трехмерный СМ, ориентированный на применение в системах реального времени с жесткими требованиями ко времени реакции, а также в вычислительных системах, осуществляющих решение СЛАУ большой размерности. Выполненная работа имеет очевидное продолжение — проведение проектных работ для последующего изготовления разработанного СМ в «железе».