БЛОЧНО-РЕКУРСИВНОЕ ПАРАЛЛЕЛЬНОЕ ПЕРЕМНОЖЕНИЕ МАТ
33
УДК 681.3
Б. Я. ШТЕЙНБЕРГ
БЛОЧНО-РЕКУРСИВНОЕ ПАРАЛЛЕЛЬНОЕ ПЕРЕМНОЖЕНИЕ МАТРИЦ
Предложен новый параллельный алгоритм перемножения матриц с количеством пересылок, меньшим, чем у существующих. Для быстрого алгоритма перемножения матриц Штрассена рассмотрена сложность по обращениям к памяти и обоснована его неэффективность для матриц практически значимой размерности. Ключевые слова: параллельные алгоритмы, гиперкуб, размещение данных, межпроцессорные пересылки, сложность по обращениям к памяти.
Введение. Новые тенденции в развитии вычислительных архитектур изменяют представления о сложности алгоритмов [1]. Если 40 лет назад время доступа к памяти было не существенно по сравнению со временем выполнения арифметических операций, то сейчас время доступа к памяти на порядок больше, а межпроцессорная пересылка требует существенно больше времени, чем доступ к памяти. Если раньше основной характеристикой сложности вычислительного алгоритма считалось количество умножений, то теперь такой характеристикой должно стать количество обращений к памяти, а для параллельных алгоритмов — количество межпроцессорных пересылок.
В настоящей работе предложен новый параллельный алгоритм перемножения матриц с количеством пересылок, меньшим, чем у существующих. Для быстрого алгоритма перемножения матриц Штрассена рассмотрена сложность по обращениям к памяти и обоснована нецелесообразность его использования для современных компьютеров; предложена процедура блочно-аффинного размещения данных. Следует отметить, что параллельные алгоритмы перемножения матриц на вычислительных системах с распределенной памятью, описанные в работах [2—5], предполагают размещение матриц полосами или блоками — это частные случаи блочно-аффинных размещений.
Приведенные в [2—5] алгоритмы требуют n3 / p параллельных умножений. Объем пересылаемых данных каждого процессора для алгоритмов Фокса и Кэннона [4, 6] составляет n2 / p1/2 . Объем пересылаемых данных каждого процессора в алгоритмах, предложенных в [3], равен n2 / p2/3 . Объем пересылаемых данных каждого процессора в предлагаемом нами
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
34 Б. Я. Штейнберг
блочно-рекурсивном алгоритме равен
n2 log2 p
p . При этом не используются рассылки данных
типа „от одного — всем“. Предполагается использование полнодоступного коммутатора или коммуникационной сети „ n -мерный куб“. Приводится схема отображения пересылок n -мерного куба на кольцевую сеть и кольцевой сети — на решетку, что позволит использовать данный алгоритм на многопроцессорных системах с соответствующей топологией.
Перемножение матриц при блочно-аффинных размещениях. Возможно размещать элементы массива в параллельной памяти произвольным образом, но необходимо, чтобы к этим элементам массива было удобно обращаться. Выбор способа размещения данных должен зависеть от исполняемой программы. Различные способы размещения данных вне связи с исполняемым кодом описаны в работе [7].
Будем предполагать, что имеется p модулей памяти. Пусть X — некоторый
m -мерный массив. Блочно-аффинные по модулю p размещения данных. Пусть натуральные числа
p, d1, d2,..., dm и целые константы s0 , s1, s2 ,..., sm зависят только от m -мерного массива X . Блочно-аффинное по модулю p размещение m -мерного массива X — это такое размеще-
ние, при котором элемент X (i1,i2 ,...,im ) находится в модуле памяти с номером
u = (]s1 / d1[ i1 +⎤⎦ s2 / d2 [i2 + ... +] sm / dm [im + s0 ) mod p .
При описанном блочно-аффинном способе размещения m -мерный массив представля-
ется как массив блоков размерностью d1, d2 ,…, dm , в котором у каждого блока все элементы
находятся в одном модуле памяти. Описанные далее способы размещения данных являются частными случаями блочно-
аффинного по модулю p размещения.
Аффинные по модулю p размещения данных — такие, при которых элемент массива
X (i1,i2 ,...,im ) находится в модуле памяти с номером u = (s1i1 + s2i2 + ... + smim + s0 ) mod p ,
где s0 , s1, s2 ,..., sm — константы, зависящие только от массива X . Число s0 будем называть сдвигом массива X . Это число соответствует номеру модуля памяти, в котором размещается нулевой элемент X (0, 0,..., 0) .
{-1, 0, +1}-размещения — в них константы s1, s2 ,..., sm принадлежат множеству {-1, 0, +1},
к ним, в частности, относятся размещение матрицы по строкам или по столбцам, или по диагоналям, или скошенная форма размещения (хранения) матрицы.
{0, +1}-размещения — в них константы s1, s2 ,..., sm принадлежат множеству {0, +1}.
К такому типу размещений относятся размещения матриц по строкам, по столбцам и в скошенном виде. Размещение по диагоналям к такому типу не относится.
Покоординатные размещения со сдвигами — это такие {0, +1}-размещения, в которых
лишь одна из величин s1, s2 ,..., sm может быть равна единице. Покоординатные размещения — это такие {0, +1}-размещения, в которых лишь одна из
величин s1, s2 ,..., sm может быть равна единице и s0 = 0 .
Размещения матрицы по строкам или по столбцам являются покоординатными. Рассмотрим фрагмент программы, вычисляющий произведение двух матриц по стандартному алгоритму
do i = 1, n
do j = 1, n
do k = 1, n
X(i,j) = X(i,j)+A(i,k)*B(k,j).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Блочно-рекурсивное параллельное перемножение матриц
35
Теорема 1. Не существует таких нетривиальных аффинных размещений массивов X , A и B в распределенной памяти, при которых параллельное произведение матриц выполняется без пересылок.
Доказательство. Предположим, что искомые аффинные размещения массивов существуют. Для каждого значения i , j и k элементы X (i, j) , A(i, k) и B(k, j) должны оказаться в
одном и том же модуле памяти. Это возможно тогда и только тогда, когда существуют такие целые числа s1, s2 , t1, t2 , u1, u2 , что выполняются равенства
s1i + s2 j = t1i + t2k = u1k + u2 j mod p . Это возможно для всех значений i, j, k лишь при выполнении равенств
s1 = s2 = t1 = t2 = u1 = u2 = 0 mod p , т.е. когда все данные находятся в одном и том же модуле памяти. Итак, задача вычисления произведения матриц не может быть распараллелена на вычислительной системе с распределенной памятью без межпроцессорных пересылок данных.
Блочно-рекурсивный параллельный алгоритм перемножения матриц. Будем считать, что количество процессорных элементов p вычислительной системы является степенью числа 4.
Рассмотрим задачу умножения квадратных матриц X = AB размерностью n на вычислительной системе с n процессорными элементами. Матрицы А и В рассмотрим как блочные
матрицы 2×2 с квадратными блоками размером (n / 2) × (n / 2) . Всего в каждой матрице полу-
чается по 4 блока. Обозначим Aij , ( i, j = 1, 2 ), аналогично обозначим блоки матрицы В.
Модуль памяти 0 А11, В11, С11
Модуль памяти 1 А12, В21, С12
Множество всех процессорных элементов разобьем на 4 равных группы с номерами 0, 1, 2, 3 (рис. 1). Разместим элементы обеих матриц по блокам — в каждой группе модулей памяти по одному.
Модуль памяти 2 А21, В12, С21
Модуль памяти 3 А22, В22, С22
— Блок Aij ( i, j = 1, 2 ) матрицы A разместим в
группе модулей памяти с номером (i −1) 2 + j −1.
— Блок Bij ( i, j = 1, 2 ) матрицы B разместим в
группе модулей памяти с номером (i −1) 2 + i −1 .
Рис. 1
Для вычисления произведения матриц в этом слу-
чае требуется два перемножения, две пересылки блоков (рис. 2, а — первое перемножение и
первая группа пересылок, б — второе перемножение и вторая группа) и два сложения (рис. 3).
а)
А11В11 А12В21
В11 В21
Модуль памяти 0 Модуль памяти 1
А21В12 А22В22
В12 В22
С11= А11В11+ А12В21 С12= А12В22+ А11В12
А11В12 А12В22
б) А11В12 А12В21
А21В11 А22В21
А21В12 А22В21
Модуль памяти 2 С21= А21В11+ А22В21
Модуль памяти 3 С22= А21В12+ А22В22
Рис. 2
Рис. 3
Перемножение блоков в каждой группе процессорных элементов (ПЭ) можно выполнять
по такой же схеме рекурсивно. В этом алгоритме возможны log4 p шагов. На каждом шаге
размерность блоков в два раза уменьшается, но количество блоков в четыре раза возрастает.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
36 Б. Я. Штейнберг
На рис. 4, а представлен пример разбиения каждой из четырех групп ПЭ на 4 группы второго
уровня; рис. 4, б — то же разбиение множества процессорных элементов на группы в двоич-
ной кодировке. Нумерация соответствует нумерации узлов четырехмерного куба.
При двойной нумерации цифра перед точкой означает номер группы ПЭ первого уров-
ня, цифра после запятой — номер группы ПЭ второго уровня (см. рис. 4, а).
а) б)
0,2 0,3 1,2 1,3
0000 0001 0100 0101
2,0 2,1 3,0 3,1
0010 0011 0110 0111
2,2 2,3 3,2 3,3
1000 1001 1100 1101
0,0 0,1 1,0 1,1
1010 1011 1110 1111
Рис. 4
Количество параллельных шагов с умножениями (скалярными) в данном алгоритме, как и в других известных параллельных алгоритмах, равно n3 / p . При наличии p процессорных
элементов и при базовом последовательном алгоритме перемножения матриц со сложностью n3 меньшего порядка сложности достичь невозможно.
Для пересылок данных в этом алгоритме удобно использовать полнодоступный коммутатор или ( log2 p )-мерный куб. Всего в алгоритме получается 2 log4 p = log2 p блочных пересылок.
На каждом шаге пересылается n2 элементов по p штук одновременно (скалярных пересы-
лочных шагов n2 / p ). Общее количество скалярных пересылочных шагов n2 log2 p / p . Если p = n , то количество одновременных скалярных пересылок равно n log2 n — это существенно меньше, чем при размещении матриц по строкам или столбцам, и меньше, чем пересылок в алгоритмах Кэннона и Фокса [4].
Перемножение прямоугольных матриц тоже может быть выполнено с помощью блочнорекурсивного алгоритма. Действительно, пусть требуется перемножить матрицу А размером n × m на матрицу В размером m × k . Обозначим q = p1/2 . Разобьем матрицу А на блоки раз-
мером (n / q) × (m / q) , а матрицу В — на (m / q) × (k / q) . Тогда задача перемножения прямо-
угольных матриц будет сведена к задаче перемножения квадратных блочных матриц размером q × q , где q = p1/2 .
Отображение пересылок между сетями различной топологии. Известна задача переноса программ с одних многопроцессорных вычислительных систем на другие. При этом вычислительные системы могут различаться не только количеством вычислительных устройств, но и коммуникационной системой [8, 9]. Тогда возникает необходимость отображения топологии (моделирования пересылок, ориентированных на одну сеть, на другой сети). Например, в работе [6] рассмотрено моделирование на гиперкубе кольцевой сети с помощью кодов Грея.
Рассмотрим задачу моделирования n -мерного куба на кольце с количеством узлов p = 2n . Напомним, что вершины n -мерного куба кодируются булевыми векторами
(x1, x2, …, xn ) длиной n .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Блочно-рекурсивное параллельное перемножение матриц
37
К сожалению, использовать коды Грея [6] для этой задачи не удобно. Действительно, для n = 3 расстояние по кольцу между вершинами, соответствующими (0,0,0) и (1,0,0), равно 1, а между (0,0,1) и (1,0,1) — 4, хотя оба этих ребра n -куба соответствуют изменению первой координаты. Во многих задачах пересылки данных по таким ребрам предполагается выполнять одновременно. В данном примере на кольце эти пересылки будут выполняться в разное время.
Для случая n = 2 задача тривиальна, поскольку двумерный куб представляет собой кольцо с 4 вершинами. Вершинам двумерного куба (0,0), (0,1), (1,1) (1,0) сопоставляются узлы кольцевой сети с номерами 0, 1, 2, 3 соответственно, т.е. ϕ(0, 0) = 0 , ϕ(0,1) = 1 , ϕ(1,1) = 2 ,
ϕ(1, 0) = 3 .
Определим сопоставление вершин n -мерного куба узлам кольцевой сети для n > 2 по следующей формуле:
( )ϕ( x1, x2 ,…, xn ) = x1 + 2x2 + … + 2n−3 xn−2 + 2n−2 ϕ( xn−1, xn ) .
Первое слагаемое в данной формуле — это число, состоящее из первых n − 2 цифр двоичного кода вершины n -мерного куба.
Рассмотрим подробнее указанное сопоставление узлов n -мерного куба узлам кольца. Сопоставление узлов n -мерного куба узлам кольца будем строить по индукции.
В предположении, что построено сопоставление узлов n -мерного куба узлам кольца для n = k, k > 1 , построим его для случая n = k +1 . Пусть узлы кольца пронумерованы от 0 до
(n −1) . Рассмотрим подмножество узлов кольца с четными номерами как кольцевую сеть
(подкольцо), в которой время пересылки между соседними узлами в 2 раза дольше, чем в исходной кольцевой сети. Это подкольцо будем называть 0-подкольцо. Аналогичным образом узлы с нечетными номерами образуют подкольцо, которое будем называть 1-подкольцом.
Пусть (x1, x2,…, xk+1) — произвольная вершина ( k +1 )-мерного куба. Рассмотрим под-
множество вершин n -мерного куба, состоящее из векторов (0, x2,…, xk+1) . Это подмножество порождает подграф, изоморфный k -мерному кубу. Узлы этого k -мерного куба сопоставим
узлам 0-подкольца. Аналогично подмножество вершин (1, x2,…, xk+1) также порождает подграф, изоморфный k -мерному кубу. Узлы этого k -мерного куба сопоставим узлам 1-подкольца.
Рассмотрим в n -мерном кубе множество ребер, соединяющих вершины вида
(x1,…, xi−1, 0, xi+1,…, xk +1) и (x1,…, xi−1,1, xi+1,…, xk +1) . Ребра n -мерного куба, принадлежащие одному такому множеству, будем называть однотипными. Всего получается n типов ребер. При реализации принципа сдваивания на n -мерном кубе используются одновременные пересылки по однотипным ребрам. Поэтому при моделировании n -мерного куба на кольце важно, чтобы такие пересылки эффективно реализовывались.
Теорема 2. При описанном выше сопоставлении вершинам n -мерного куба узлов кольцевой сети расстояние между концевыми узлами всех однотипных ребер одинаково.
Доказательство может быть проведено методом математической индукции. Следствие. Пересылки на кольце, соответствующие пересылкам по однотипным ребрам n -мерного куба, требуют одинакового времени. В работе [6] рассмотрен способ представления на гиперкубе двумерной решетки, число узлов на сторонах которой является степенью двойки. Представим кольцо на решетке, что позволит на решетке моделировать пересылки n -мерного куба. Если сеть имеет гамильтонов цикл (цикл, проходящий через каждую вершину строго один раз), то такой цикл и можно рассматривать как кольцевую сеть. Покажем, как построить гамильтонов цикл на прямоугольной решетке типа „тор“.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
38 Б. Я. Штейнберг
Допустим, что у прямоугольной решетки хотя бы одна из сторон содержит четное количество узлов. Пусть размерность этой решетки m(2k) . Тогда номера узлов, приведенные на
рис. 5 (без окаймления), составляют гамильтонов цикл на этой решетке. Пусть теперь прямоугольная решетка содержит нечетное количество узлов
(2m + 1) × (2k + 1) . В этом случае в решетке выделяется подрешетка размерностью (2m) × (2k) ,
узлы которой нумеруются в указанной выше последовательности (см. рис. 5, с окаймлением). На многомерном торе кольцо можно моделировать по правилу индукции, для двумер-
ного тора отображение построено. Предположим, что оно существует для m -мерного тора и построим его для ( m +1)-мерного тора. Пусть W1W2 …Wm+1 — ( m +1)-мерный тор. Тогда, по предположению индукции, на m -мерном торе W2 …Wm+1 можно смоделировать кольцо, которое обозначим C . Тогда исходный ( m +1)-мерный тор можно представить в виде двумерного тора W1C и на этом двумерном торе смоделировать кольцо.
У многокольцевых соединений [10], как правило, одно из колец является гамильтоновым циклом.
1 2m (2m+1)
m(2k)
m(2k)+1
2 2m–1 (2m+2)
m(2k)–1
m(2k)+2
…
(m–1)
(m+2) 3m–1 3m+2
M
(m+1) 3m 3m+1 … m(2k–1)
m(2k+1)
(m+1)(2k+1) … … … … m(2k+1)+2 m(2k+1)+1
Рис. 5
О сложности обращений к памяти алгоритма Штрассена. Сложность стандартного алгоритма (в котором, согласно определению, строки одной матрицы скалярно умножаются
на столбцы другой) равна n3 . В данном случае n3 — это количество умножений. У алгоритма
Штрассена [5, 11, 12] сложность равна n2,81 (более точно, показатель степени равен log2 7 ).
Этот специфический алгоритм дал импульс развитию теории сложности. Напомним, что сложность алгоритма — это количество операций как функция от количества входных данных, а сложность задачи — это сложность самого быстрого алгоритма. При этом под операциями подразумевались арифметические операции, как правило — умножения, т.е. требующие существенно большего времени, чем сложение. Позднее были найдены алгоритмы перемножения матриц еще менее сложные, чем алгоритм Штрассена. Существуют алгоритмы, представляющие собой комбинацию алгоритма Штрассена и стандартного алгоритма перемножения матриц [13]. Сложность задачи перемножения матриц пока не определена. Есть результаты, показывающие существование таких задач, для которых самого эффективного алгоритма не существует [14] — не исключено, что задача перемножения матриц является такой же.
Рассмотрим реальные преимущества алгоритма Штрассена по сравнению со стандартным. В 1988 г. на кафедре алгебры и дискретной математики механико-математического факультета Ростовского государственного университета (сейчас Южный федеральный университет) был программно реализован алгоритм Штрассена (студенческая работа, руководитель доц. Козак А.В.). Выяснилось, что алгоритм Штрассена начинает давать преимущества по сравнению со стандартным алгоритмом начиная с размерности матриц n = 32 . В 2002 г. опять был реализован алгоритм Штрассена (тоже студенческая работа, руководитель — автор настоящей работы). В этот раз преимущество алгоритма Штрассена начиналось с размерности матриц n = 512 . Попытаемся объяснить этот эффект. Примерно за 15 лет компьютеры стали совсем другими. Изменилось соотношение времени выполнения арифметических опе-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Блочно-рекурсивное параллельное перемножение матриц
39
раций и времени обращения к памяти. А соотношение количества этих операций в алгоритме Штрассена и в стандартном алгоритме различно.
Умножение вещественных чисел современным компьютером может быть в 15 раз быстрее, чем считывание этих чисел из оперативной памяти. Попробуем оценить сложность обращений к памяти алгоритма Штрассена.
Предполагается, что размерность матрицы является степенью двойки: n = 2k . Исходные и результирующая матрицы представляются как блочные 2 × 2 матрицы с блоками размерностью n / 2 . Тогда произведение матриц имеет вид
C=AB C = {Cij }, A = {Aij }, B = {Bij }, i, j = 1, 2 .
Здесь Cij , Aij , Bij (i, j = 1, 2) — блоки матриц С, А и В соответственно.
По алгоритму Штрассена вычисляются семь вспомогательных матриц (размером n / 2 ):
M1 = ( A12 – A22 )( B21 + B22 ), M 2 = ( A11 + A22 )( B11 + B22 ), M3 = ( A11 – A21 ) ( B11 + B12 ), M 4 = ( A11 + A12 ) B22 , M5 = A11 ( B12 – B22 ) , M6 = A22 ( B21 – B11 ), M7 = ( A21 + A22 ) B11.
(1)
После этого элементы результирующей матрицы вычисляются по формулам
C11 = M1 + M 2 – M 4 + M6 , C12 = M 4 + M5 , C21 = M 6 + M 7 , C22 = M 2 – M3 + M5 – M7 .
(2)
В формулах (1) перемножаются матрицы меньшей размерности, что тоже может выполняться рекурсивно по такому же алгоритму.
Обозначим через F (n) количество операций обращения к памяти (чтение или запись) в алгоритме Штрассена.
В формулах (1), (2) матрицы размером n / 2 встречаются 37 раз. При этом формулы (1) содержат умножения матриц размером n / 2 , каждое из которых потребует F (n / 2) обращений к памяти. Итого можно записать рекуррентную формулу
F (n) = 7F (n / 2) + 37(n / 2)2. Учитывая, что F (1) = 3 , из рекуррентной формулы можно получить аналитическую
F (n) = 7F (n / 2) + 37(n / 2)2 = 7{7F (n / 4) + 37(n / 4)2} + 37(n / 2)2 =
= 7k F (1) + 37{7k−1(n / 2k )2 + … + 7(n / 4)2 + (n / 2)2} =
= 7k ⋅ 3 + 37 / 7n2{(7 / 4)k +… + (7 / 4)2 + 7 / 4 } =
= 7k ⋅ 3 + 37 / 7n2 7 / 3((7 / 4)k – 1) =
= 7k ⋅ 3 + 37 / 3n2 ((7 / 4)k – 1).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
40 Б. Я. Штейнберг
Пренебрегая в скобках слагаемым (–1) и учитывая, что k = log2 n , получим, что количество обращений к памяти в алгоритме Штрассена равно
F (n) ≈ 7k 3 + 37 / 3n2 (7 / 4)k = 7k 3 + 37 / 37k = 7k 46 / 3 ≈ n2,8146 / 3 .
Рассмотрим количество обращений к памяти в стандартном алгоритме. Каждый элемент результирующей матрицы вычисляется как скалярное произведение двух векторов (строки первой матрицы на столбец второй). Такое скалярное произведение вычисляется накоплением суммы произведений в регистре и требует n2 чтений и одну запись в память. Итого 2n3 + n2 ≈ 2n3 обращений к памяти.
Асимптотически количество обращений к памяти алгоритма Штрассена меньше, чем у стандартного алгоритма. Проследим, начиная с какой размерности матриц n в алгоритме Штрассена количество обращений к памяти будет меньше, чем в стандартном алгоритме — 2n3 = n2,81 46 / 3 , n ≈ 25 000 .
Численный эксперимент для матриц такой размерности не представляется разумным, поскольку матрицы такой размерности могут быть размещены уже не в оперативной, а во внешней памяти.
Следует отметить, что рекурсивность алгоритма Штрассена тоже может вызвать дополнительные обращения к памяти (обращения к стеку, в котором хранятся параметры вызовов). Глубина рекурсивных вызовов равна k = log2 n . Если столько параметров вызовов не поместится в регистровой или кэш-памяти, то они будут помещаться в оперативной.
Итак, алгоритм Штрассена асимптотически эффективнее, но преимущество перед стандартным алгоритмом начинается с задач, размерность которых превышает практическую значимость.
СПИСОК ЛИТЕРАТУРЫ
1. Пипер Ш., Пол Д., Скотт М. Новая эра в оценке производительности компьютерных систем // Открытые системы. 2007. № 9. С. 52—58.
2. Корнеев В. Д. Параллельное программирование в MPI. М.—Ижевск: Институт компьютерных исследований, 2003. 304 с.
3. Gupta H., Sadayappan P. Communication Efficient Matrix-Multiplication on Hypercubes // Parallel Computing. 1996. N 22. P. 75—99.
4. Гергель В. П. Введение в методы параллельного программирования. Образовательный комплекс. Н. Новгород: ННГУ, 2005.
5. Прангишвили И. В., Виленкин С. Я., Медведев И. Л. Параллельные вычислительные системы с общим управлением. М.: Энергоатомиздат, 1983. 312 с.
6. Гергель В. П., Стронгин Р. Г. Основы параллельных вычислений для многопроцессорных вычислительных систем. Н. Новгород: ННГУ, 2003.
7. Метлицкий Е. А., Каверзнев В. В. Системы параллельной памяти. Теория, проектирование, применение. Л.: ЛГУ, 1989. 240 с.
8. Treleaven P.C. Parallel architecture overview // Parallel Computing. 1988. N 8. P. 59—70.
9. Корнеев В. В. Параллельные вычислительные системы. М.: Нолидж, 1999. 311 с.
10. Подлазов В. С. Свойства мультикольцевых и гиперкубовых коммутаторов на произвольных перестановках // PACO’2001 Тр. Междунар. конф. „Параллельные вычисления и задачи управления“. М.: ИПУ РАН, 2001. С. 152—164.
11. Strassen V. Gaussian Elimination is not Optimal // Numer. Math. 1969. Vol. 13, N 4. P. 354—356.
12. Ахо А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979. 536 с.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Проблемы распараллеливания в некоторых локальных задачах
41
13. Елфимова Л. Д. Быстрый клеточный метод умножения матриц // Кибернетика и системный анализ. 2008. № 3. С. 55—59.
14. Разборов Л. А. P=NP? — проблема: взгляд из 90-х // Математические события XX в. М.: ФАЗИС, 2003. С. 399—416.
Борис Яковлевич Штейнберг
Сведения об авторе — д-р техн. наук; Южный федеральный университет, Ростов-на-Дону;
кафедра алгебры и дискретной математики; зав. кафедрой; E-mail: borsteinb@mail.ru
Рекомендована институтом
Поступила в редакцию 10.03.09 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
УДК 681.3
Б. Я. ШТЕЙНБЕРГ
БЛОЧНО-РЕКУРСИВНОЕ ПАРАЛЛЕЛЬНОЕ ПЕРЕМНОЖЕНИЕ МАТРИЦ
Предложен новый параллельный алгоритм перемножения матриц с количеством пересылок, меньшим, чем у существующих. Для быстрого алгоритма перемножения матриц Штрассена рассмотрена сложность по обращениям к памяти и обоснована его неэффективность для матриц практически значимой размерности. Ключевые слова: параллельные алгоритмы, гиперкуб, размещение данных, межпроцессорные пересылки, сложность по обращениям к памяти.
Введение. Новые тенденции в развитии вычислительных архитектур изменяют представления о сложности алгоритмов [1]. Если 40 лет назад время доступа к памяти было не существенно по сравнению со временем выполнения арифметических операций, то сейчас время доступа к памяти на порядок больше, а межпроцессорная пересылка требует существенно больше времени, чем доступ к памяти. Если раньше основной характеристикой сложности вычислительного алгоритма считалось количество умножений, то теперь такой характеристикой должно стать количество обращений к памяти, а для параллельных алгоритмов — количество межпроцессорных пересылок.
В настоящей работе предложен новый параллельный алгоритм перемножения матриц с количеством пересылок, меньшим, чем у существующих. Для быстрого алгоритма перемножения матриц Штрассена рассмотрена сложность по обращениям к памяти и обоснована нецелесообразность его использования для современных компьютеров; предложена процедура блочно-аффинного размещения данных. Следует отметить, что параллельные алгоритмы перемножения матриц на вычислительных системах с распределенной памятью, описанные в работах [2—5], предполагают размещение матриц полосами или блоками — это частные случаи блочно-аффинных размещений.
Приведенные в [2—5] алгоритмы требуют n3 / p параллельных умножений. Объем пересылаемых данных каждого процессора для алгоритмов Фокса и Кэннона [4, 6] составляет n2 / p1/2 . Объем пересылаемых данных каждого процессора в алгоритмах, предложенных в [3], равен n2 / p2/3 . Объем пересылаемых данных каждого процессора в предлагаемом нами
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
34 Б. Я. Штейнберг
блочно-рекурсивном алгоритме равен
n2 log2 p
p . При этом не используются рассылки данных
типа „от одного — всем“. Предполагается использование полнодоступного коммутатора или коммуникационной сети „ n -мерный куб“. Приводится схема отображения пересылок n -мерного куба на кольцевую сеть и кольцевой сети — на решетку, что позволит использовать данный алгоритм на многопроцессорных системах с соответствующей топологией.
Перемножение матриц при блочно-аффинных размещениях. Возможно размещать элементы массива в параллельной памяти произвольным образом, но необходимо, чтобы к этим элементам массива было удобно обращаться. Выбор способа размещения данных должен зависеть от исполняемой программы. Различные способы размещения данных вне связи с исполняемым кодом описаны в работе [7].
Будем предполагать, что имеется p модулей памяти. Пусть X — некоторый
m -мерный массив. Блочно-аффинные по модулю p размещения данных. Пусть натуральные числа
p, d1, d2,..., dm и целые константы s0 , s1, s2 ,..., sm зависят только от m -мерного массива X . Блочно-аффинное по модулю p размещение m -мерного массива X — это такое размеще-
ние, при котором элемент X (i1,i2 ,...,im ) находится в модуле памяти с номером
u = (]s1 / d1[ i1 +⎤⎦ s2 / d2 [i2 + ... +] sm / dm [im + s0 ) mod p .
При описанном блочно-аффинном способе размещения m -мерный массив представля-
ется как массив блоков размерностью d1, d2 ,…, dm , в котором у каждого блока все элементы
находятся в одном модуле памяти. Описанные далее способы размещения данных являются частными случаями блочно-
аффинного по модулю p размещения.
Аффинные по модулю p размещения данных — такие, при которых элемент массива
X (i1,i2 ,...,im ) находится в модуле памяти с номером u = (s1i1 + s2i2 + ... + smim + s0 ) mod p ,
где s0 , s1, s2 ,..., sm — константы, зависящие только от массива X . Число s0 будем называть сдвигом массива X . Это число соответствует номеру модуля памяти, в котором размещается нулевой элемент X (0, 0,..., 0) .
{-1, 0, +1}-размещения — в них константы s1, s2 ,..., sm принадлежат множеству {-1, 0, +1},
к ним, в частности, относятся размещение матрицы по строкам или по столбцам, или по диагоналям, или скошенная форма размещения (хранения) матрицы.
{0, +1}-размещения — в них константы s1, s2 ,..., sm принадлежат множеству {0, +1}.
К такому типу размещений относятся размещения матриц по строкам, по столбцам и в скошенном виде. Размещение по диагоналям к такому типу не относится.
Покоординатные размещения со сдвигами — это такие {0, +1}-размещения, в которых
лишь одна из величин s1, s2 ,..., sm может быть равна единице. Покоординатные размещения — это такие {0, +1}-размещения, в которых лишь одна из
величин s1, s2 ,..., sm может быть равна единице и s0 = 0 .
Размещения матрицы по строкам или по столбцам являются покоординатными. Рассмотрим фрагмент программы, вычисляющий произведение двух матриц по стандартному алгоритму
do i = 1, n
do j = 1, n
do k = 1, n
X(i,j) = X(i,j)+A(i,k)*B(k,j).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Блочно-рекурсивное параллельное перемножение матриц
35
Теорема 1. Не существует таких нетривиальных аффинных размещений массивов X , A и B в распределенной памяти, при которых параллельное произведение матриц выполняется без пересылок.
Доказательство. Предположим, что искомые аффинные размещения массивов существуют. Для каждого значения i , j и k элементы X (i, j) , A(i, k) и B(k, j) должны оказаться в
одном и том же модуле памяти. Это возможно тогда и только тогда, когда существуют такие целые числа s1, s2 , t1, t2 , u1, u2 , что выполняются равенства
s1i + s2 j = t1i + t2k = u1k + u2 j mod p . Это возможно для всех значений i, j, k лишь при выполнении равенств
s1 = s2 = t1 = t2 = u1 = u2 = 0 mod p , т.е. когда все данные находятся в одном и том же модуле памяти. Итак, задача вычисления произведения матриц не может быть распараллелена на вычислительной системе с распределенной памятью без межпроцессорных пересылок данных.
Блочно-рекурсивный параллельный алгоритм перемножения матриц. Будем считать, что количество процессорных элементов p вычислительной системы является степенью числа 4.
Рассмотрим задачу умножения квадратных матриц X = AB размерностью n на вычислительной системе с n процессорными элементами. Матрицы А и В рассмотрим как блочные
матрицы 2×2 с квадратными блоками размером (n / 2) × (n / 2) . Всего в каждой матрице полу-
чается по 4 блока. Обозначим Aij , ( i, j = 1, 2 ), аналогично обозначим блоки матрицы В.
Модуль памяти 0 А11, В11, С11
Модуль памяти 1 А12, В21, С12
Множество всех процессорных элементов разобьем на 4 равных группы с номерами 0, 1, 2, 3 (рис. 1). Разместим элементы обеих матриц по блокам — в каждой группе модулей памяти по одному.
Модуль памяти 2 А21, В12, С21
Модуль памяти 3 А22, В22, С22
— Блок Aij ( i, j = 1, 2 ) матрицы A разместим в
группе модулей памяти с номером (i −1) 2 + j −1.
— Блок Bij ( i, j = 1, 2 ) матрицы B разместим в
группе модулей памяти с номером (i −1) 2 + i −1 .
Рис. 1
Для вычисления произведения матриц в этом слу-
чае требуется два перемножения, две пересылки блоков (рис. 2, а — первое перемножение и
первая группа пересылок, б — второе перемножение и вторая группа) и два сложения (рис. 3).
а)
А11В11 А12В21
В11 В21
Модуль памяти 0 Модуль памяти 1
А21В12 А22В22
В12 В22
С11= А11В11+ А12В21 С12= А12В22+ А11В12
А11В12 А12В22
б) А11В12 А12В21
А21В11 А22В21
А21В12 А22В21
Модуль памяти 2 С21= А21В11+ А22В21
Модуль памяти 3 С22= А21В12+ А22В22
Рис. 2
Рис. 3
Перемножение блоков в каждой группе процессорных элементов (ПЭ) можно выполнять
по такой же схеме рекурсивно. В этом алгоритме возможны log4 p шагов. На каждом шаге
размерность блоков в два раза уменьшается, но количество блоков в четыре раза возрастает.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
36 Б. Я. Штейнберг
На рис. 4, а представлен пример разбиения каждой из четырех групп ПЭ на 4 группы второго
уровня; рис. 4, б — то же разбиение множества процессорных элементов на группы в двоич-
ной кодировке. Нумерация соответствует нумерации узлов четырехмерного куба.
При двойной нумерации цифра перед точкой означает номер группы ПЭ первого уров-
ня, цифра после запятой — номер группы ПЭ второго уровня (см. рис. 4, а).
а) б)
0,2 0,3 1,2 1,3
0000 0001 0100 0101
2,0 2,1 3,0 3,1
0010 0011 0110 0111
2,2 2,3 3,2 3,3
1000 1001 1100 1101
0,0 0,1 1,0 1,1
1010 1011 1110 1111
Рис. 4
Количество параллельных шагов с умножениями (скалярными) в данном алгоритме, как и в других известных параллельных алгоритмах, равно n3 / p . При наличии p процессорных
элементов и при базовом последовательном алгоритме перемножения матриц со сложностью n3 меньшего порядка сложности достичь невозможно.
Для пересылок данных в этом алгоритме удобно использовать полнодоступный коммутатор или ( log2 p )-мерный куб. Всего в алгоритме получается 2 log4 p = log2 p блочных пересылок.
На каждом шаге пересылается n2 элементов по p штук одновременно (скалярных пересы-
лочных шагов n2 / p ). Общее количество скалярных пересылочных шагов n2 log2 p / p . Если p = n , то количество одновременных скалярных пересылок равно n log2 n — это существенно меньше, чем при размещении матриц по строкам или столбцам, и меньше, чем пересылок в алгоритмах Кэннона и Фокса [4].
Перемножение прямоугольных матриц тоже может быть выполнено с помощью блочнорекурсивного алгоритма. Действительно, пусть требуется перемножить матрицу А размером n × m на матрицу В размером m × k . Обозначим q = p1/2 . Разобьем матрицу А на блоки раз-
мером (n / q) × (m / q) , а матрицу В — на (m / q) × (k / q) . Тогда задача перемножения прямо-
угольных матриц будет сведена к задаче перемножения квадратных блочных матриц размером q × q , где q = p1/2 .
Отображение пересылок между сетями различной топологии. Известна задача переноса программ с одних многопроцессорных вычислительных систем на другие. При этом вычислительные системы могут различаться не только количеством вычислительных устройств, но и коммуникационной системой [8, 9]. Тогда возникает необходимость отображения топологии (моделирования пересылок, ориентированных на одну сеть, на другой сети). Например, в работе [6] рассмотрено моделирование на гиперкубе кольцевой сети с помощью кодов Грея.
Рассмотрим задачу моделирования n -мерного куба на кольце с количеством узлов p = 2n . Напомним, что вершины n -мерного куба кодируются булевыми векторами
(x1, x2, …, xn ) длиной n .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Блочно-рекурсивное параллельное перемножение матриц
37
К сожалению, использовать коды Грея [6] для этой задачи не удобно. Действительно, для n = 3 расстояние по кольцу между вершинами, соответствующими (0,0,0) и (1,0,0), равно 1, а между (0,0,1) и (1,0,1) — 4, хотя оба этих ребра n -куба соответствуют изменению первой координаты. Во многих задачах пересылки данных по таким ребрам предполагается выполнять одновременно. В данном примере на кольце эти пересылки будут выполняться в разное время.
Для случая n = 2 задача тривиальна, поскольку двумерный куб представляет собой кольцо с 4 вершинами. Вершинам двумерного куба (0,0), (0,1), (1,1) (1,0) сопоставляются узлы кольцевой сети с номерами 0, 1, 2, 3 соответственно, т.е. ϕ(0, 0) = 0 , ϕ(0,1) = 1 , ϕ(1,1) = 2 ,
ϕ(1, 0) = 3 .
Определим сопоставление вершин n -мерного куба узлам кольцевой сети для n > 2 по следующей формуле:
( )ϕ( x1, x2 ,…, xn ) = x1 + 2x2 + … + 2n−3 xn−2 + 2n−2 ϕ( xn−1, xn ) .
Первое слагаемое в данной формуле — это число, состоящее из первых n − 2 цифр двоичного кода вершины n -мерного куба.
Рассмотрим подробнее указанное сопоставление узлов n -мерного куба узлам кольца. Сопоставление узлов n -мерного куба узлам кольца будем строить по индукции.
В предположении, что построено сопоставление узлов n -мерного куба узлам кольца для n = k, k > 1 , построим его для случая n = k +1 . Пусть узлы кольца пронумерованы от 0 до
(n −1) . Рассмотрим подмножество узлов кольца с четными номерами как кольцевую сеть
(подкольцо), в которой время пересылки между соседними узлами в 2 раза дольше, чем в исходной кольцевой сети. Это подкольцо будем называть 0-подкольцо. Аналогичным образом узлы с нечетными номерами образуют подкольцо, которое будем называть 1-подкольцом.
Пусть (x1, x2,…, xk+1) — произвольная вершина ( k +1 )-мерного куба. Рассмотрим под-
множество вершин n -мерного куба, состоящее из векторов (0, x2,…, xk+1) . Это подмножество порождает подграф, изоморфный k -мерному кубу. Узлы этого k -мерного куба сопоставим
узлам 0-подкольца. Аналогично подмножество вершин (1, x2,…, xk+1) также порождает подграф, изоморфный k -мерному кубу. Узлы этого k -мерного куба сопоставим узлам 1-подкольца.
Рассмотрим в n -мерном кубе множество ребер, соединяющих вершины вида
(x1,…, xi−1, 0, xi+1,…, xk +1) и (x1,…, xi−1,1, xi+1,…, xk +1) . Ребра n -мерного куба, принадлежащие одному такому множеству, будем называть однотипными. Всего получается n типов ребер. При реализации принципа сдваивания на n -мерном кубе используются одновременные пересылки по однотипным ребрам. Поэтому при моделировании n -мерного куба на кольце важно, чтобы такие пересылки эффективно реализовывались.
Теорема 2. При описанном выше сопоставлении вершинам n -мерного куба узлов кольцевой сети расстояние между концевыми узлами всех однотипных ребер одинаково.
Доказательство может быть проведено методом математической индукции. Следствие. Пересылки на кольце, соответствующие пересылкам по однотипным ребрам n -мерного куба, требуют одинакового времени. В работе [6] рассмотрен способ представления на гиперкубе двумерной решетки, число узлов на сторонах которой является степенью двойки. Представим кольцо на решетке, что позволит на решетке моделировать пересылки n -мерного куба. Если сеть имеет гамильтонов цикл (цикл, проходящий через каждую вершину строго один раз), то такой цикл и можно рассматривать как кольцевую сеть. Покажем, как построить гамильтонов цикл на прямоугольной решетке типа „тор“.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
38 Б. Я. Штейнберг
Допустим, что у прямоугольной решетки хотя бы одна из сторон содержит четное количество узлов. Пусть размерность этой решетки m(2k) . Тогда номера узлов, приведенные на
рис. 5 (без окаймления), составляют гамильтонов цикл на этой решетке. Пусть теперь прямоугольная решетка содержит нечетное количество узлов
(2m + 1) × (2k + 1) . В этом случае в решетке выделяется подрешетка размерностью (2m) × (2k) ,
узлы которой нумеруются в указанной выше последовательности (см. рис. 5, с окаймлением). На многомерном торе кольцо можно моделировать по правилу индукции, для двумер-
ного тора отображение построено. Предположим, что оно существует для m -мерного тора и построим его для ( m +1)-мерного тора. Пусть W1W2 …Wm+1 — ( m +1)-мерный тор. Тогда, по предположению индукции, на m -мерном торе W2 …Wm+1 можно смоделировать кольцо, которое обозначим C . Тогда исходный ( m +1)-мерный тор можно представить в виде двумерного тора W1C и на этом двумерном торе смоделировать кольцо.
У многокольцевых соединений [10], как правило, одно из колец является гамильтоновым циклом.
1 2m (2m+1)
m(2k)
m(2k)+1
2 2m–1 (2m+2)
m(2k)–1
m(2k)+2
…
(m–1)
(m+2) 3m–1 3m+2
M
(m+1) 3m 3m+1 … m(2k–1)
m(2k+1)
(m+1)(2k+1) … … … … m(2k+1)+2 m(2k+1)+1
Рис. 5
О сложности обращений к памяти алгоритма Штрассена. Сложность стандартного алгоритма (в котором, согласно определению, строки одной матрицы скалярно умножаются
на столбцы другой) равна n3 . В данном случае n3 — это количество умножений. У алгоритма
Штрассена [5, 11, 12] сложность равна n2,81 (более точно, показатель степени равен log2 7 ).
Этот специфический алгоритм дал импульс развитию теории сложности. Напомним, что сложность алгоритма — это количество операций как функция от количества входных данных, а сложность задачи — это сложность самого быстрого алгоритма. При этом под операциями подразумевались арифметические операции, как правило — умножения, т.е. требующие существенно большего времени, чем сложение. Позднее были найдены алгоритмы перемножения матриц еще менее сложные, чем алгоритм Штрассена. Существуют алгоритмы, представляющие собой комбинацию алгоритма Штрассена и стандартного алгоритма перемножения матриц [13]. Сложность задачи перемножения матриц пока не определена. Есть результаты, показывающие существование таких задач, для которых самого эффективного алгоритма не существует [14] — не исключено, что задача перемножения матриц является такой же.
Рассмотрим реальные преимущества алгоритма Штрассена по сравнению со стандартным. В 1988 г. на кафедре алгебры и дискретной математики механико-математического факультета Ростовского государственного университета (сейчас Южный федеральный университет) был программно реализован алгоритм Штрассена (студенческая работа, руководитель доц. Козак А.В.). Выяснилось, что алгоритм Штрассена начинает давать преимущества по сравнению со стандартным алгоритмом начиная с размерности матриц n = 32 . В 2002 г. опять был реализован алгоритм Штрассена (тоже студенческая работа, руководитель — автор настоящей работы). В этот раз преимущество алгоритма Штрассена начиналось с размерности матриц n = 512 . Попытаемся объяснить этот эффект. Примерно за 15 лет компьютеры стали совсем другими. Изменилось соотношение времени выполнения арифметических опе-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Блочно-рекурсивное параллельное перемножение матриц
39
раций и времени обращения к памяти. А соотношение количества этих операций в алгоритме Штрассена и в стандартном алгоритме различно.
Умножение вещественных чисел современным компьютером может быть в 15 раз быстрее, чем считывание этих чисел из оперативной памяти. Попробуем оценить сложность обращений к памяти алгоритма Штрассена.
Предполагается, что размерность матрицы является степенью двойки: n = 2k . Исходные и результирующая матрицы представляются как блочные 2 × 2 матрицы с блоками размерностью n / 2 . Тогда произведение матриц имеет вид
C=AB C = {Cij }, A = {Aij }, B = {Bij }, i, j = 1, 2 .
Здесь Cij , Aij , Bij (i, j = 1, 2) — блоки матриц С, А и В соответственно.
По алгоритму Штрассена вычисляются семь вспомогательных матриц (размером n / 2 ):
M1 = ( A12 – A22 )( B21 + B22 ), M 2 = ( A11 + A22 )( B11 + B22 ), M3 = ( A11 – A21 ) ( B11 + B12 ), M 4 = ( A11 + A12 ) B22 , M5 = A11 ( B12 – B22 ) , M6 = A22 ( B21 – B11 ), M7 = ( A21 + A22 ) B11.
(1)
После этого элементы результирующей матрицы вычисляются по формулам
C11 = M1 + M 2 – M 4 + M6 , C12 = M 4 + M5 , C21 = M 6 + M 7 , C22 = M 2 – M3 + M5 – M7 .
(2)
В формулах (1) перемножаются матрицы меньшей размерности, что тоже может выполняться рекурсивно по такому же алгоритму.
Обозначим через F (n) количество операций обращения к памяти (чтение или запись) в алгоритме Штрассена.
В формулах (1), (2) матрицы размером n / 2 встречаются 37 раз. При этом формулы (1) содержат умножения матриц размером n / 2 , каждое из которых потребует F (n / 2) обращений к памяти. Итого можно записать рекуррентную формулу
F (n) = 7F (n / 2) + 37(n / 2)2. Учитывая, что F (1) = 3 , из рекуррентной формулы можно получить аналитическую
F (n) = 7F (n / 2) + 37(n / 2)2 = 7{7F (n / 4) + 37(n / 4)2} + 37(n / 2)2 =
= 7k F (1) + 37{7k−1(n / 2k )2 + … + 7(n / 4)2 + (n / 2)2} =
= 7k ⋅ 3 + 37 / 7n2{(7 / 4)k +… + (7 / 4)2 + 7 / 4 } =
= 7k ⋅ 3 + 37 / 7n2 7 / 3((7 / 4)k – 1) =
= 7k ⋅ 3 + 37 / 3n2 ((7 / 4)k – 1).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
40 Б. Я. Штейнберг
Пренебрегая в скобках слагаемым (–1) и учитывая, что k = log2 n , получим, что количество обращений к памяти в алгоритме Штрассена равно
F (n) ≈ 7k 3 + 37 / 3n2 (7 / 4)k = 7k 3 + 37 / 37k = 7k 46 / 3 ≈ n2,8146 / 3 .
Рассмотрим количество обращений к памяти в стандартном алгоритме. Каждый элемент результирующей матрицы вычисляется как скалярное произведение двух векторов (строки первой матрицы на столбец второй). Такое скалярное произведение вычисляется накоплением суммы произведений в регистре и требует n2 чтений и одну запись в память. Итого 2n3 + n2 ≈ 2n3 обращений к памяти.
Асимптотически количество обращений к памяти алгоритма Штрассена меньше, чем у стандартного алгоритма. Проследим, начиная с какой размерности матриц n в алгоритме Штрассена количество обращений к памяти будет меньше, чем в стандартном алгоритме — 2n3 = n2,81 46 / 3 , n ≈ 25 000 .
Численный эксперимент для матриц такой размерности не представляется разумным, поскольку матрицы такой размерности могут быть размещены уже не в оперативной, а во внешней памяти.
Следует отметить, что рекурсивность алгоритма Штрассена тоже может вызвать дополнительные обращения к памяти (обращения к стеку, в котором хранятся параметры вызовов). Глубина рекурсивных вызовов равна k = log2 n . Если столько параметров вызовов не поместится в регистровой или кэш-памяти, то они будут помещаться в оперативной.
Итак, алгоритм Штрассена асимптотически эффективнее, но преимущество перед стандартным алгоритмом начинается с задач, размерность которых превышает практическую значимость.
СПИСОК ЛИТЕРАТУРЫ
1. Пипер Ш., Пол Д., Скотт М. Новая эра в оценке производительности компьютерных систем // Открытые системы. 2007. № 9. С. 52—58.
2. Корнеев В. Д. Параллельное программирование в MPI. М.—Ижевск: Институт компьютерных исследований, 2003. 304 с.
3. Gupta H., Sadayappan P. Communication Efficient Matrix-Multiplication on Hypercubes // Parallel Computing. 1996. N 22. P. 75—99.
4. Гергель В. П. Введение в методы параллельного программирования. Образовательный комплекс. Н. Новгород: ННГУ, 2005.
5. Прангишвили И. В., Виленкин С. Я., Медведев И. Л. Параллельные вычислительные системы с общим управлением. М.: Энергоатомиздат, 1983. 312 с.
6. Гергель В. П., Стронгин Р. Г. Основы параллельных вычислений для многопроцессорных вычислительных систем. Н. Новгород: ННГУ, 2003.
7. Метлицкий Е. А., Каверзнев В. В. Системы параллельной памяти. Теория, проектирование, применение. Л.: ЛГУ, 1989. 240 с.
8. Treleaven P.C. Parallel architecture overview // Parallel Computing. 1988. N 8. P. 59—70.
9. Корнеев В. В. Параллельные вычислительные системы. М.: Нолидж, 1999. 311 с.
10. Подлазов В. С. Свойства мультикольцевых и гиперкубовых коммутаторов на произвольных перестановках // PACO’2001 Тр. Междунар. конф. „Параллельные вычисления и задачи управления“. М.: ИПУ РАН, 2001. С. 152—164.
11. Strassen V. Gaussian Elimination is not Optimal // Numer. Math. 1969. Vol. 13, N 4. P. 354—356.
12. Ахо А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979. 536 с.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10
Проблемы распараллеливания в некоторых локальных задачах
41
13. Елфимова Л. Д. Быстрый клеточный метод умножения матриц // Кибернетика и системный анализ. 2008. № 3. С. 55—59.
14. Разборов Л. А. P=NP? — проблема: взгляд из 90-х // Математические события XX в. М.: ФАЗИС, 2003. С. 399—416.
Борис Яковлевич Штейнберг
Сведения об авторе — д-р техн. наук; Южный федеральный университет, Ростов-на-Дону;
кафедра алгебры и дискретной математики; зав. кафедрой; E-mail: borsteinb@mail.ru
Рекомендована институтом
Поступила в редакцию 10.03.09 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2009. Т. 52, № 10