CC BY
49
2013 Вычислительные методы в дискретной математике №3(21)
УДК 519.7
О РЕАЛИЗАЦИИ АЛГОРИТМА КОППЕРСМИТА ДЛЯ ДВОИЧНЫХ МАТРИЧНЫХ ПОСЛЕДОВАТЕЛЬНОСТЕЙ НА ВЫЧИСЛИТЕЛЯХ
КЛАСТЕРНОГО ТИПА
А. С. Рыжов Лаборатория ТВП, г. Москва, Россия E-mail: ryzhovalexander@mail.ru
Рассматривается задача реализации алгоритма Копперсмита, вычисляющего векторные аннулирующие многочлены для матричных последовательностей, на современных 64-разрядных ЭВМ. Рассмотрены вопросы представления данных для случая последовательностей бинарных матриц с точки зрения снижения трудоёмкости алгоритма. Предложены способы эффективного распараллеливания алгоритма для реализации на ЭВМ с многоядерными процессорами, а также для выполнения алгоритма на вычислителях кластерного типа.
Ключевые слова: матричные последовательности, алгоритм Копперсмита.
Введение
Во многих алгоритмах вычислительной алгебры возникает задача нахождения решений сильно разреженной системы однородных линейных уравнений (СОЛУ). Например, при нахождении дискретного логарифма в поле GF(2n) (при p = 2n — І — простом) методом Копперсмита [1] необходимо найти решения такой системы над простым полем Zp.
В настоящий момент для решения данной задачи обычно используется один из двух алгоритмов: алгоритм Ланцоша — Монтгомери [2] или алгоритм Видемана — Копперсмита [3]. Оба алгоритма имеют свои преимущества и недостатки, а также ограничения на используемые вычислители. Алгоритм Ланцоша — Монтгомери реализован, например, в пакете GGNFS. Он ориентирован на выполнение на одном вычислителе кластерного типа с высокоскоростным внутренним соединением узлов. Алгоритм Ви-демана — Копперсмита, напротив, позволяет наиболее трудоёмкие этапы выполнять на нескольких не связанных друг с другом кластерах.
Блочный алгоритм Видемана — Копперсмита является обобщением алгоритма Видемана поиска решения разреженной СОЛУ [4], ориентированным на вычислительную технику, работающую со словами в несколько байтов. Один из этапов алгоритма Виде-мана — нахождение аннулирующего многочлена двоичной последовательности с помощью алгоритма Берлекэмпа — Мэсси, и его непосредственное обобщение на многомерный случай является, по-видимому, трудной задачей. Копперсмит решил эту проблему, сведя задачу к поиску векторных приближений Паде и разработав новый алгоритм для построения таких приближений.
В данной работе рассматриваются вопросы реализации алгоритма Копперсмита построения векторных приближений Паде при решении СОЛУ над полем GF(2). Приводится описание алгоритма Видемана — Копперсмита нахождения решений разреженных СОЛУ, подробное описание метода Копперсмита построения векторных приближений Паде, особенности реализации этого алгоритма на ЭВМ. Отдельно рассмотрены особенности реализации алгоритма Копперсмита на вычислителях кластерного типа.
1. Блочный алгоритм Видемана — Копперсмита решения систем однородных линейных уравнений
Описание алгоритма Видемана — Копперсмита нахождения решений разреженных систем однородных линейных уравнений можно найти, например, в [3]. Приведём основные шаги алгоритма, чтобы показать, где используется алгоритм Копперсмита построения векторных приближений Паде. Пусть Р = ОЕ(2) — поле из двух элементов, А — матрица размеров М х N над полем Р, и пусть необходимо найти нетривиальное (т. е. ненулевое) решение СОЛУ
Ах = 0. (1)
Заметим, что не ставится задача нахождения всех решений системы (1) или приведения матрицы А к специальному виду. Чтобы гарантировать наличие нетривиальных решений, будем считать, что М < N.
Пусть т, п Є N — параметры алгоритма, т ^ п. Пусть, далее, Z Є Рт,м — случайная матрица размера т х М над полем Р, В Є Рм,м — матрица из первых М столбцов матрицы А решаемой СОЛУ (1), У Є Рм,п — матрица, составленная из столбцов матрицы А с номерами М + 1, М + 2,..., М + п (пусть М + п ^ N). Последовательно выполняем следующие три этапа:
1) Вычислим первые Ь + 1 элементов последовательности матриц а = а0аі... размера т х п над полем Р, где а^ = ZBгУ (здесь Ь — параметр метода; см. ниже).
2) Найдём векторный аннулирующий многочлен для последовательности а, т. е.
й
такие векторы /0,... , / Є Рп, что ^ а+ fj = 0 Є Рдля всех і ^ 0.
3=0
й
3) Из последнего равенства получаем, что вектор ^ В3У/ лежит в ядре любой
3=0
й
матрицы ZBг, і = 0,1,... Отсюда с большой вероятностью ^ В3У/ = 0, т. е.
3=0
(й-1
у- В 3 У / •
3=0 3
/0
дополненный нулями до длины N, является решением системы (1) (напомним, У состоит из столбцов матрицы А, идущими после подматрицы В). Если век-
й-1
тор w равен нулю, то ^ В3У/3+1 = 0, и выполняются аналогичные рассуж-
3=0
дения.
Замечание 1. Необходимо отметить, что для последовательности а указанного
вида векторный аннулирующий многочлен, вычисляемый на втором этапе алгоритма,
всегда существует. При т = п = 1 таким многочленом является любой аннулирующий
многочлен матрицы В — например характеристический или минимальный. В общем й
случае равенство ^ В3У/ = 0 можно рассматривать как систему из М однородных
3=0
линейных уравнений от (^ + 1)п неизвестных элементов векторов {/} над полем Р, и при ^ ^ М/п эта система имеет ненулевое решение — искомый векторный аннулирующий многочлен.
Исходный алгоритм Видемана, опубликованный в работе [4], является частным случаем приведённого выше блочного алгоритма Видемана — Копперсмита при т = п = 1. В алгоритме Видемана на втором этапе необходимо найти аннулирующий многочлен
последовательности а = а0а1 . . . элементов поля Р, при этом достаточно вычислить первые 2М + 1 элементов, так как ранг последовательности не превосходит ранг матрицы В. Для поиска аннулирующего многочлена можно использовать, например, алгоритм Берлекэмпа — Мэсси.
Для нахождения векторного аннулирующего многочлена последовательности матриц Копперсмит предложил новый алгоритм, основанный на последовательном вычислении приближений Паде. При этом в результате работы алгоритма может получиться до п различных независимых аннулирующих многочленов, что в результате дает до п линейно независимых решений системы (1).
Для нахождения векторного аннулирующего многочлена достаточно иметь первые Ь +1 элементов последовательности а, где Ь = [М/т] + [М/п] + А, А ~ 100. Это следует из набора линейных соотношений, а также из теоретико-вероятностных соображений [3].
По сравнению с методом Видемана, в алгоритме Видемана — Копперсмита на первом и третьем этапах требуется выполнить существенно меньшее число итераций. Для удобства реализации на современных ЭВМ числа т и п выбирают, как правило, кратными длине машинного слова: например, т = 128, п = 64. Расчет трудоёмкости приведён автором метода в [3]. Алгоритм дает выигрыш только в случае сильной разреженности матрицы системы: для первого и третьего этапов трудоёмкость пропорциональна М2и>, где и> — средний вес строки матрицы В, т. е. среднее число ненулевых элементов в строке. Трудоёмкость второго этапа при использовании алгоритма Коп-персмита равна 0(М2т).
Далее приводится подробное изложение алгоритма Копперсмита и рассматриваются вопросы, связанные с его эффективной реализацией на вычислителях кластерного типа.
2. Алгоритм Копперсмита вычисления векторных аннулирующих многочленов
Описание алгоритма Копперсмита вычисления векторных аннулирующих многочленов приведено автором алгоритма в [3]. Алгоритм сформулирован достаточно громоздко в терминах отдельных элементов коэффициентов степенных рядов над кольцом матриц. Здесь приводится описание алгоритма, в основном схожее с описанием в работе [5] и более удобное для дальнейшего изложения.
Итак, для последовательности матриц над конечным полем (не обязательно СЕ(2)) необходимо найти векторный аннулирующий многочлен. Поскольку алгоритм Коппер-смита, как было сказано ранее, выдаёт сразу несколько решений, будем рассматривать задачу нахождения матричного аннулирующего многочлена: для последовательности а = а0 а1... матриц размера т х п над конечным полем Р необходимо найти матрицы /0,... , /й Є Рп,г размера п х г, одновременно не равные нулю, такие, что
Уі ^ 0 ^ аі+3/ = 0тхг. (2)
3=0
Будем изначально предполагать существование решения.
Замечание 2. Последовательность а = а0а1..., обладающую свойством (2), в общем случае нельзя называть линейной рекуррентной последовательностью, поскольку при т = п множество возможных значений коэффициентов этой последовательности не образует структуру кольца. Нельзя также говорить о многочлене вида /0хй + /1хй-1 + • • • + /й. Тем не менее далее соответствующие объекты будем
называть многочленами и степенными рядами для удобства и наглядности изложения. Например, будем говорить «степенной ряд Н(ж) = ^ а^жг» или «многочлен Р(ж) = /0жй + Лжй-1 + • • • + /й», а также рассматривать их «произведение», если размеры коэффициентов допускают перемножение. При этом какие-либо алгебраические свойства этих объектов использоваться не будут.
Замечание 3. Как было отмечено выше, при решении исходной задачи — нахождения решений СОЛУ (1) —построенная последовательность а = (2ВгУ : г = 0,1,...)
й
заведомо обладает свойством (2), поскольку система ^ В3У/3- = 0, состоящая из гМ
з=о
уравнений от г(1 + 1)п неизвестных (коэффициентов матриц /0,... ,/й) при 1 ^ М/п имеет ненулевое решение.
Рассмотрим произведение ряда Н(ж) = ^ а*жг и многочлена Р(ж) = /0жй + /1 жй-1 + + • • • + /й:
сю й
к
ак-й+3/3 I ж
г=0 3=0 к=0 \ 3=0
(здесь а^ = 0 при г < 0). Ясно, что если выполняется условие (2), то ряд Н(ж)Р(ж) является многочленом степени не выше 1 — 1, то есть его коэффициенты при жй, жй+1,. . . равны нулю. При этом Р(ж) называется правым порождающим многочленом последовательности а. Обратно, если для некоторых многочленов Р(ж) € Рп,г[ж], ^(ж) € Рт,г[ж] выполняется Н(ж)Р(ж) = ^(ж), то коэффициенты многочлена Р(ж) являются решением поставленной задачи.
Если для некоторых многочленов Р (ж), ^(ж) и рядов Н (ж), С (ж) выполняется равенство
Н (ж)Р (ж) = ^(ж) + ж*С (ж), (3)
то многочлены Р и ^ называют приближениями Паде порядка £ для степенного ряда Н(ж). Алгоритм Копперсмита итеративно строит приближения Паде порядка £ = £0,£0 + 1,... Как было оговорено выше, последовательность а обладает правым порождающим многочленом, поэтому при достаточно большом £ получим правый порождающий многочлен этой последовательности.
Приведём описание алгоритма Копперсмита. Через Ет будем обозначать единичную матрицу размера т х т. Положим А(ж) = (Н(ж)| — Ет)тх(т+П), С*(ж) =
Р*(ж)
?*(^ / (т+га)х(т+га)
носильно равенству А(ж)С*(ж) = ж* С* (ж) (здесь С (ж) из равенства (3) заменено на С*(ж), чтобы показать, что этот матричный многочлен определяется номером шага алгоритма £ = 0,1, 2,...). Для столбцов матрицы С*(ж) = (С*д(ж),... , С4,(т+га)(ж)) введём целочисленные векторы = ($*д,... , $*,(т+п)), элементы которых вычисляются в процессе работы алгоритма и в некотором роде соответствуют степеням этих векторных многочленов. (Под степенью столбца матрицы над кольцом многочленов будем понимать максимум степеней многочленов, являющихся элементами этого столбца:
deg (^1(ж),... , дк(ж))Т = тах {deg^(ж)}.) А именно:
ге{1,...,к>
— степень столбца С4,3- (ж) не превосходит соответствующее значение $4,3-;
— при сложении двух столбцов С*,31 (ж) и С*,32 (ж) результату будет соответствовать максимум их значений тах{$4,31, $4,32};
— при умножении столбца на ж соответствующее ему значение увеличивается на 1.
( ) . Тогда А(ж)С*(ж) = Н(ж)Р*(ж) — ф*(ж) и равенство (3) рав-
V ^*(ж) / (т+га) X (т+га)
На каждом шаге алгоритма Копперсмита последовательно вычисляются пары (Gt(x), Ct(x)), для которых выполняются следующие три условия:
A(x)Gt(x) = (х); (4)
rank Gt(0) = m; (5)
m+n
deg Gt,j(x) ^ ^t,j, j = 1, 2,..., m + n, ^ = tm. (6)
j=i
Инициализация значений (Go(x), Co(x)) выполняется следующим образом: G0(x) = = E(m+n)x(m+n), C0(x) = (H(x)|Em). Значения £o = (^o,i,... ,^o,m+n) задаются равными нулю. Выполнение условий (4) - (6) для t = 0 очевидно.
При выполнении очередного шага t необходимо так изменить матрицу Gt (x), чтобы можно было выделить одну степень х. Для этого достаточно n столбцов младшего коэффициента (т. е. Gt (0)) сделать нулевыми, а остальные m столбцов домножить на х.
Формально, на шаге t по паре (Gt(x), Ct(x)) вычисляется пара (Gt+1(x), Gt+1 (х)) следующим образом:
1) Вычислим матрицу перехода Tt размера (m + n) х (m + n), которая представляет собой невырожденные линейные преобразования столбцов, приводящие матрицу Gt(0) к ступенчатому виду: Gt(0)rt = (0|E). При этом к столбцу Gtj(х) можно прибавлять только те столбцы Gtj2 (х), соответствующие значения $tj2 которых не больше (этого можно добиться за счёт перестановки столбцов; при этом значения $tj также переставляются).
2) Вычислим Gt+1 (х) = G'^T'D, где D = diag(1,... , 1,х,... , х) (единицы стоят на первых n элементах диагонали).
3) Вычислим С*+1(х) = Gt^T'D/T.
4) Положим i1+i,j Л ^ ’ ес"и j e{1--"-n>
l^t,j + 1, если j е {n +1,...,m + n}.
Заметим, что вычисление матрицы Tt возможно в силу выполнения условия (5). Деление на х при вычислении С'+1(х) возможно в силу того, что первые n столбцов матрицы С'(0)т' равны нулю, т. е. делятся на х, а оставшиеся m столбцов домножаются на х при умножении на матрицу D.
Проверим выполнение условий (4) - (6). Выполнение условия (4) для t + 1 очевидно (если оно выполняется для t): левая и правая части равенства (4) для t домножены справа на одну и ту же полиномиальную матрицу TtD. Условие (5) для t + 1 следует из того, что последние m столбцов матрицы Gt(0)т являются единичными. Условие (6) выполняется в силу особенностей построения матрицы перехода Tt, а также умножения G'^T' на матрицу D, содержащую ровно m элементов х на диагонали.
Признаком останова алгоритма является наличие в матрице Gt^) достаточного (требуемого) числа столбцов Gtj (х), для которых выполняется неравенство
N
t — >-------+ Д 2, (7)
m
где Д2 — наперёд заданный постоянный параметр; Д2 ^ 10-100. Если необходимо найти n решений, то в силу условия (6) останов выполнится по крайней мере через t ~ N/m + N/n шагов.
Столбцы, для которых выполняется условие (7), представляют собой решение задачи. Поскольку нас интересуют только сами правые порождающие многочлены Ftj (х), нет необходимости вычислять на каждом шаге всю матрицу Gt^) =
Таким образом, на каждом шаге необходимо выполнить следующие действия:
— с помощью невырожденных линейных преобразований столбцов привести матрицу С*(0) к ступенчатому виду С*(0)т* = (0|Е);
— умножить полиномиальную матрицу С* (ж) на матрицу т*Д;
— умножить полиномиальную матрицу С* (ж) на т*Д/ж.
Замечание 4. За счёт выбора матрицы т*Д на каждом шаге ровно т столбцов матрицы С(ж) увеличивают свою степень. Именно это позволяет к определённому моменту набрать достаточное число столбцов, удовлетворяющих условию (7).
Оценим трудоёмкость алгоритма Копперсмита. На каждом шаге наиболее трудоёмким действием является умножение полиномиальных матриц С*(ж) и С* (ж) на скалярную матрицу перехода т*. Степень матрицы С*(ж), очевидно, не превосходит і. Степень матрицы С* (ж) не превосходит Ь — і, где Ь = deg Со (ж) = deg Н(ж) — длина отрезка последовательности, для которой ищется правый порождающий многочлен. Таким образом, при т ^ п трудоёмкость алгоритма имеет порядок 0(Ь2т3). При Ь = |_М/т] + |_М/п] + А, как в алгоритме Видемана — Копперсмита нахождения решений СОЛУ, получаем трудоёмкость порядка 0(М2т).
Рассмотрим некоторые подходы, позволяющие эффективно реализовать алгоритм Копперсмита нахождения векторных аннулирующих многочленов матричной последовательности над полем СЕ(2) на современной 64-разрядной вычислительной технике. Будем считать, что размеры последовательности т и п кратны 64: т = 64т/, п = 64п;.
С точки зрения оптимизации реализации алгоритма на ЭВМ большую роль играет способ представления данных. Естественным выбором является представление, позволяющее упростить наиболее трудоёмкие операции, выполняемые на каждом шаге алгоритма. Кроме того, при использовании вычислителей кластерного типа необходимо минимизировать объём данных, которыми обмениваются узлы кластера на каждом шаге.
Наиболее трудоёмкими операциями, выполняемыми на каждом шаге алгоритма Копперсмита, являются:
— умножение (справа) матрицы С*(ж) на матрицу т*Д/ж;
— умножение (справа) матрицы С* (ж) на матрицу т*Д.
Выделение нулевого коэффициента С*(0) полиномиальной матрицы С* (ж) при любом её естественном представлении не составляет труда. Операция приведения С*(0) при небольших значениях т и п также не является трудоёмкой по сравнению с умножением матриц С* (ж) и С*(ж) на матрицу перехода т*. Вектор степеней £*+1 = (^*+1,1,..., ^*+1,т+п) вычисляется за 0(т + п) операций.
Матрица перехода т* является скалярной матрицей, поэтому умножение полиномиальных матриц С* (ж) и С*(ж) на т* справа сводится к отдельному умножению каждой строки каждого коэффициента этих матриц на т*. Поэтому с точки зрения снижения трудоёмкости этой операции, а также с учётом использования 64-разрядной ЭВМ естественно представлять полиномиальные матрицы С* (ж) и С* (ж) в виде массивов из т и п (соответственно) полиномиальных строк (длины этих массивов соответству-
, а достаточно выполнять действия только с первыми п стро-
ками.
3. Вопросы реализации алгоритма Копперсмита на ЭВМ
3.1. Представление данных
ют степеням матриц С*(ж) и С*(ж)), каждую полиномиальную строку — в виде массива строк-коэффициентов, а каждую строку-коэффициент (длины т + п = 64(т' + п')) — в виде массива из т' + п' 64-разрядных машинных слов.
Умножение на ж выполняется над последними т = 64т/ столбцами матрицы С* (ж) (соответственно деление на ж — над первыми п = 64п' столбцами матрицы С*(ж)), поэтому эту операцию можно выполнить автоматически, записывая результат умножения строки-коэффициента на матрицу перехода т* со сдвигом: последние т/ машинных слов результата умножения строки-коэффициента при жк записываются в массив, соответствующий коэффициенту при жк+1 (для С* (ж) —первые п' машинных слов результата умножения строки-коэффициента при жк записываются в массив, соответствующий коэффициенту при жк-1).
Предложенный способ хранения данных позволяет записывать результаты очередного шага алгоритма на место предыдущих результатов: Сі+1 (ж) вместо С*(ж), а С*+і(ж) вместо С*(ж), что снимает необходимость постоянного выделения дополнительных объёмов оперативной памяти и, в свою очередь, снижает нагрузку на системный диспетчер памяти.
3.2. Быстрое умножение полиномиальных матриц
н а м а т р и ц у п е р е х о д а
Поскольку основной элементарной операцией на каждом шаге алгоритма Коппер-смита является умножение большого числа строк длины т+п на одну и ту же матрицу перехода т* размера (т + п) х (т + п), для снижения трудоёмкости всего алгоритма необходимо прежде всего оптимизировать эту операцию.
Один из эффективных методов умножения строки на матрицу небольших размеров связан с использованием так называемых Іоокир-таблиц. Идея эта используется давно [6, гл. 6]; она достаточно проста и заключается в следующем. Пусть необходимо реализовать умножение строки а длиной 8к битов на двоичную матрицу С размеров 8к х 8к. Будем считать, что строка а хранится в массиве байтов длины к. Представим массив а в виде суммы к байтовых массивов длины к, в которых только один
к
байт может быть ненулевым: а = ^ а^, где а^ = (0,... , 0, *, 0,... , 0) Є (ненуле-
І=1
вой элемент — на і-й позиции). Тогда произведение строки на матрицу раскладывается
к
в сумму к произведений: аС = ^ а^С. Если для матрицы С заранее вычислить все
І=1
возможные произведения вида а^С — а их число равно 256к, — то умножение строки длины 8к битов на матрицу 8к х 8к сводится к сложению к строк длины 8к битов.
В нашем случае роль матрицы С играет матрица перехода т* размером (т + п) х х(т + п) битов и к = 8(т' + п'). Для реализации приведённого алгоритма необходимо составить таблицу из 211(т' + п') строк длины т + п битов. Составление этой таблицы является предварительным этапом и выполняется один раз на каждом шаге алгоритма Копперсмита после вычисления матрицы перехода т*, поэтому вклад этой процедуры в общую трудоёмкость алгоритма ничтожно мал по сравнению с последующим использованием таблицы при умножении строк С*(ж) и С* (ж) на т*. Умножение одной строки на матрицу т* в этом случае выполняется за к = 8(т' + п') побитовых сложений массивов из т' + п' 64-разрядных слов.
3.3. Многопоточная реализация алгоритма Рассмотрим вопросы, связанные с параллельным выполнением алгоритма Копперсмита несколькими потоками на многоядерной/многопроцессорной ЭВМ с общей памятью.
Наиболее трудоёмкими операциями алгоритма являются операции умножения полиномиальных матриц С*(ж) и С*(ж) на скалярную матрицу перехода т*, причём эти операции выполняются независимо для различных строк матриц С*(ж) и С* (ж), а также для матриц-коэффициентов при различных степенях ж. Это предоставляет практически неограниченные возможности для распараллеливания данных операций.
При организации многопоточных вычислений применительно к рассматриваемому случаю необходимо учитывать следующие особенности:
— выбранный способ представления данных (строки матриц С* (ж) и С* (ж) представляются в виде отдельных массивов скалярных строк-коэффициентов при 1, ж, ж2,...);
— число строк матриц С*(ж) и С*(ж) равно, как правило, т = 64т' и п = 64п' соответственно, где т' и п' — натуральные числа;
— число ядер в современных ЭВМ с серийными процессорами составляет порядка 12-16 ядер на 1 процессор.
С учётом сказанного, достаточно реализовать параллельное умножение различных строк матриц С*(ж) и С*(ж), так как их общее количество даже при т' = п' = 1 в разы превышает возможное число вычислительных ядер. В этом случае многопоточная реализация одного шага алгоритма Копперсмита выглядит следующим образом:
1) управляющий поток вычисляет матрицу С4(0), а также матрицу перехода т*;
2) каждый рабочий поток выполняет умножение выделенных ему строк матрицы С*(ж) (С*(ж)) на т*Д/ж (т*Л);
3) управляющий поток вычисляет вектор степеней и дожидается завершения умножения строк матриц С* (ж) и Д4(ж) всеми рабочими потоками.
Важно отметить, что при большой длине Ь входной матричной последовательности трудоёмкость вычисления матрицы перехода т* и вектора степеней $4+1 ничтожно мала по сравнению с умножением матриц С* (ж) и С* (ж) на т*, поэтому «простой» рабочих потоков является незначительным.
При использовании большего числа вычислительных ядер — скажем, порядка 1000 — вычисления отдельных строк матриц С*+1(ж) и С*+1(ж) можно также распределять по нескольким ядрам, разделяя между ними вычисления по диапазону степеней ж. Например, при умножении С*(ж) на матрицу перехода т* каждый коэффициент каждой строки С* (ж) умножается на т* независимо от остальных: если разделить С* (ж) по строкам С*(ж) = (С*,1(ж),..., С*,п(ж))Т, а каждую строку С^(ж) по степеням ж: *
С*,г(ж) = ^2 С*,г,^-ж3, то умножение выполняется отдельно для каждого коэффициента з=о
*
при различных степенях ж: С*,г(ж)т* = Е(С*,г,3-т*)ж3. В этом случае необходимо лишь
3=0
учитывать перекрытия на границах диапазонов степеней ж, выделенных различным вычислительным ядрам, поскольку одновременно с умножением на т* выполняется умножение последних т столбцов С*(ж)т* на ж, что достигается за счёт записи результата умножения на т* со сдвигом на одну степень ж; аналогичным образом выполняется деление первых п столбцов С*(ж)т* на ж (см. п. 3.1).
Практически единственным препятствием для эффективного использования большого количества вычислительных ядер при реализации алгоритма Копперсмита яв-
ляется большое число обращений к оперативной памяти. Действительно, с вычислительной точки зрения на каждом шаге выполняются простейшие операции, связанные с умножением строки на матрицу, а при использовании lookup-таблиц это умножение сводится к побитному сложению предварительно вычисленных массивов машинных слов. В то же время выполняется полное чтение и перезапись больших по объёму матриц Gt(x) и Ct(x). Другими словами, основными операциями являются чтение из памяти и запись в память. Поскольку современные модули памяти не могут обрабатывать запросы от нескольких ядер одновременно, распараллеливание алгоритма на большое число ядер становится неэффективным.
Одним из способов решения этой проблемы является использование вычислителей, в которых отдельные модули памяти соответствуют различным процессорам. В этом случае ядра одного процессора обращаются к своей доли оперативной памяти и выполнение алгоритма происходит более эффективно. Частным случаем является использование вычислителей кластерного типа.
3.4. Реализация алгоритма на вычислителях
к л а с т е р н о г о т и п а
При организации высокопроизводительных вычислений обычно используются вычислители кластерного типа, или кластеры. Кластер — это разновидность распределённой вычислительной системы, состоящей из нескольких связанных между собой компьютеров (узлов) и используемой как единый ресурс. В качестве связующего элемента выступают Ethernet, InfiniBand или любые другие относительно недорогие сети. Распределение вычислений по отдельным узлам выполняет один или несколько выделенных управляющих узлов.
В отличие от обычных многопоточных реализаций, при реализации алгоритмов на кластерах необходимо учитывать отсутствие общей для всех узлов оперативной памяти. Тем не менее в случае распараллеливания алгоритма Копперсмита всё достаточно просто: каждый узел хранит свою порцию строк матриц (x) и Ct(x), а один шаг
алгоритма выглядит следующим образом:
1) управляющий узел получает от каждого рабочего узла часть строк матрицы Ct(0), вычисляет матрицу перехода Tt и рассылает её копии каждому рабочему узлу;
2) каждый рабочий узел выполняют умножение выделенных ему строк матрицы Ct(x) (Gt(x)) на TtD/x (rtD);
3) управляющий узел вычисляет вектор степеней $t+1 и дожидается завершения вычислений на всех рабочих узлах.
Объём данных, передаваемых по сети на каждом шаге, минимален и составляет несколько килобайтов: 512m/(m/ + n/) байтов при пересылке матрицы C't(0)mx(m+n) и 512(m/ + n/)2 байтов при пересылке матрицы (rt)(m+n)x(m+n). Трудоёмкость приведения Ct(0) к ступенчатому виду мала по сравнению с остальными вычислениями, поэтому простой узлов кластера при этой операции является несущественным.
На отдельных узлах кластера возможно дальнейшее распараллеливание вычислений по отдельным ядрам, как это описано в предыдущем пункте. Тем не менее из-за особенностей работы оперативной памяти распределение вычислений по большому количеству ядер внутри одного узла не является эффективным.
Необходимо также отметить, что на каждом шаге алгоритма используется весь объём данных, полученных на предыдущем шаге, поэтому выполнение алгоритма на нескольких не связанных между собой вычислителях не представляется возможным.
4. Результаты экспериментов
Приведём результаты экспериментальных исследований времени работы алгоритма Копперсмита. Реализация алгоритма выполнена на языке C+—+, для компиляции программного кода использован 64-разрядный компилятор среды разработки Microsoft Visual Studio 2008. В качестве вычислителя использована ЭВМ с процессором Intel Xeon (2 процессора 2,53 ГГц по 4 ядра каждый) под управлением 64-разрядной ОС Windows 7. В таблице представлено время работы алгоритма при различных размерах задачи. Рассматривалась последовательность, получаемая в качестве результата работы первого этапа блочного алгоритма Копперсмита при решении СОЛУ из M уравнений; m и n — параметры алгоритма; L — длина последовательности; threads — число потоков.
Время работы алгоритма
M m n L threads Время, с
100 000 2 1 2443 1 71
100 000 2 1 2443 2 91
100 000 2 1 2443 4 90
100 000 2 1 2443 8 93
1 000 000 2 1 23537 1 6578
1 000 000 4 2 11818 1 12969
1 000 000 4 2 11818 2 8849
1 000 000 4 2 11818 4 5726
1 000 000 4 2 11818 8 4700
Из таблицы видно, что при небольших размерах задачи распараллеливание не приводит к ускорению работы алгоритма, а наоборот, замедляет его. Это связано с тем, что на каждом шаге алгоритма необходимо выполнить небольшой объём вычислений, и синхронизация потоков на каждом шаге отнимает время, сравнимое с временем вычислений. При увеличении размера входа видно, что параллельная реализация работает быстрее. Однако увеличение числа потоков в 2 раза не приводит к соответствующему снижению времени работы алгоритма из-за особенностей работы оперативной памяти.
5. Субквадратичные алгоритмы
Трудоёмкость алгоритма Копперсмита составляет O(m3L2) элементарных операций, где L — длина матричной последовательности, а m и n — размеры её элементов (m ^ n). В то же время предложено большое число алгоритмов построения векторных аннулирующих многочленов с асимптотически меньшей трудоёмкостью. Например, если поле P допускает быстрое преобразование Фурье (БПФ), то изложенный в работе [7] алгоритм решает эту задачу за O(m3L log2 L) операций. В работе [5] предложена субквадратичная версия алгоритма Копперсмита — рекурсивный алгоритм Копперсмита — Томэ, который также использует БПФ и имеет трудоёмкость O (m3L log L) операций.
Несмотря на то, что эти алгоритмы имеют в асимптотике меньшую трудоёмкость, при относительно небольших размерах входа алгоритм Копперсмита работает значительно быстрее, поскольку не имеет большой константы в формуле трудоёмкости. Кроме того, у алгоритма Копперсмита на порядок ниже требуемый объём оперативной памяти, что также является немаловажным.
Заключение
Рассмотренный в данной работе алгоритм Копперсмита используется для построения векторных аннулирующих многочленов для матричных последовательностей над
полем. В частности, этот алгоритм является составляющей частью блочного алгоритма Видемана — Копперсмита нахождения решений больших разреженных систем линейных уравнений над полем.
Проведены исследования, связанные с эффективной реализацией алгоритма Коп-персмита на современной вычислительной технике. Предложен способ представления входных и промежуточных данных, позволяющий эффективно распараллеливать наиболее трудоёмкие операции каждого шага алгоритма. Рассмотрен подход к оптимизации основной операции — умножения полиномиальных матриц большой степени на скалярную матрицу — за счёт выполнения предварительных вычислений.
Отдельное внимание уделено многопоточной реализации алгоритма, а также выполнению алгоритма на вычислителях кластерного типа. Приведены экспериментальные результаты.
ЛИТЕРАТУРА
1. Coppersmith D. Fast evaluation of logarithms in fields of characteristic two // IEEE Trans. Inform. Theory. 1984. V.IT-30(4). P. 587-594.
2. Montgomery P. L. A block Lanczos algorithm for finding dependencies over GF(2) // EUROCRYPT’95. LNCS. 1995. V.921. P. 106-120.
3. Coppersmith D. Solving linear equations over GF(2) via block Wiedemann algorithm // Math. Comp. 1994. No. 62(205). P. 333-350.
4. Wiedemann D. H. Solving sparse linear equations over finite fields // IEEE Trans. Inform. Theory. 1986. V. IT-32(1). P. 54-62.
5. Thome E. Subquadratic computation of vector generating polynomials and improvement of the block Wiedemann algorithm // J. Symbolic Comput. 2002. No. 33. P. 757-775.
6. Ахо А., Хопкрофт Д., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979. 536 с.
7. Beckerman B. and Labahn G. A uniform approach for the fast computation of matrix-type Pade approximants // SIAM J. Matrix Anal. Appl. 1994. No. 15(3). P. 804-823.
