Дж. Деммель - Вычислительная линейная алгебра, страница 19

Иными словами, ни при какой организации матричного умножения (с выполнением прежнего числа 2пз арифметических операций) нельзя получить значение 0, большее, чем 0(т/М). С другой стороны, наш краткий анализ игнорирует ряд практических вопросов: 1. Реальная программа должна уметь работать и с прямоугольными матрицами, где при оптимальном разбиении блоки не обязательно являются квадратными.

2. Наилучшие размеры подматриц сильно зависят от структуры кзш-памятн и регистров компьютера. 3. Могут существовать специальные аппаратные возможности для выполнения умножения и сложения в одном цикле. Может также быть возможным одновременное выполнение нескольких независимых операций умножения- сложения. Подробное обсуждение этих вопросов для конкретной высокопроизводительной машины, а именно рабочей станции 1ВМ НЯ6000/590, можно найти в [1] (см. также РАКАЬЬЕЬ НОМЕРКЕ или Ы1р://уттувт.гэ6000лЬш.соуп/геэоцтсе/ /сес)тпо1ойу/еээЬЫтп1).

На рис. 2.5 показаны скорости трех уровней ВЬАБ'а для этой машины. По горизонтальной оси откладывается размерность матрицы, а по вертикальной — скорость в мегафлопах. Пиковая производительность машины составляет 266 мегафлопов. Верхняя кривая (с максимальным значением НЗ2: уровни 1, 2 и 3 пакета ВЫЗ х с о 200 Ф и 160 8 а о 100 ст 00 100 200 600 400 600 600 Размерность векторов и матриц Рис. 2.5. Производительность пакета ВЬАЯ иа компьютере 1ВМ ВЯ 0000/590.

2.6. Блочные алгоритмы как средство повышения производительности 81 тз 250 мегафлопов) соответствует умножению квадратных матриц. Средняя кривая (с максимальным значением гз 100 мегафлопов) относится к операции умножения квадратной матрицы на вектор. Нижняя кривая (с максимальным значением 75 мегафлопов) отвечает операции вахру. Заметим, что для больших матриц скорость увеличивается. Такой рост наблюдается всегда, поэтому мы постараемся строить алгоритмы с как можно ббльшими матрицами во внутренних матричных умножениях.

Оба приведенных выше алгоритма для умножения матриц выполняют 2пз арифметических операций. Оказывается, что существуют другие алгоритмы матричного умножения, в которых число операций гораздо меньше. Первым из таких алгоритмов был метод Штрассена [3); он же допускает наиболее простое объяснение. В этом методе сомножители рассматриваются как блочные 2 х 2-матрицы, умножение которых сводится к семи умножениям и 18 сложениям блоков половинного порядка. Эта идея используется рекурсивно, что дает аСНМПтатИЧЕСКуЮ СЛОжНОСтЬ П1'гг 7 = Пг В1 ВМЕСТО П'.

Алгоритм 2.8. Метод Штрассена длл матричного умножения С = $2таззеп/А,В,п/ /' Алгоритм вычисляет произведение С = А з В п х п-матриц А и В Предполагаетсл, что и есть степень двойки */ 1/п=1 вычислить С = А * В /' умножение чисел */ е1зе Представить матриць1 в блочном виде А = А11 А12 В = В11 В12 где блоки АВ и В, имеют порядок п/2 Р1 —— Бстаззеп(А12 — А г В21+ Вгг п/2) Рг = Бьтаззеп(А11 + Агг, Вм + Вгг, и/2) Рг = %тавзен(А11 — Агп В11 + В12 п/2) Р4 = $1таззеп(А11 -'т Агг, Вгг, п/2) РВ = Бгтаззеп(А11, Вш — Вг2, и/2) Р — — БЬ.аззеп(Агг, Вгп — В11, п/2) Р, = Бстаззеп(А21 + Агг, В11, и/2) С11 = Р1 + Рг — РВ + РВ Сш =РВ+РВ С21 = РВ + Рт С22 = Р2 — Рг + Рь — Р7 возвратить матрицу С = ~ С Сы Сгг 21 22 епд 1/ Проверка по индукции того, что этот алгоритм правильно перемножает матрицы, является несложным, хотя и довольно скучным упражнением (см.

вопрос 2.21). Покажем, что его сложность 0(п"гг "). Обозначим через Т(п) число сложений, вычитаний и умножений, выполняемых в алгоритме. Поскольку алгоритм предусматривает 7 обращений к самому себе с матрицами порядка 82 Глава 2. Решение линейных уравнений и/2 и 18 сложений п/2 х п/2-матриц, можно составить рекуррентное соотношение Т(п) = 7Т(п/2) + 18(п/2)з. Замена переменного т = 1окэ и дает новую рекурсию Т(т) = 7Т(т — 1)+18(2 г)з, где Т(т) = Т(2 ). Можйо проверить, что решение этого линейного разностного уравнения относительно Т имеет вид Т(т) = 0(7 ) = 0(п"з' ) Значение алгоритма Штрассена состоит не только в этой меньш й асимптотической сложности, но и в сведении задачи к решению меньших п дзадач.

! Продолжая редукцию, в конечном счете получим подзадачи, помещающиеся в быстрой памяти; как только это достигнуто, можно пользоваться стайдартным алгоритмом матричного умножения. Для некоторых компьютероя и матриц относительно высокого порядка этим путем было достигнуто ускорение умножения [22]. Недостатками данного подхода являются большие затраты памяти и несколько худшая численная устойчивость, хотя и достаточная для многих задач [77]. Существует ряд еще более быстрых алгоритмов матричного умножения; в настоящее время рекордный по скорости алгоритм, предложенный Виноградом и Копперсмитом, имеет сложность 0(пз зтв). Однако подобные алгоритмы становятся зкономичней алгоритма Штрассена лишь для столь больших значений и, какие не встречаются в практических задачах.

Обзор этой тематики см. в [195]. 2.6.3. Реорганизация гауссова исключения с целью использования уровня 3 пакета ВЬАБ Мы модифицируем гауссово исключение с тем, чтобы вначале использовать уровень 2 пакета ВЬАБ, а затем уровень 3. Для простоты, предположим, что выбор главного элемента не нужен. В действительности, алгоритм 2.4 уже является алгоритмом уровня 2, поскольку ббльшая часть работы производится в его второй строке А(1+ 1: п,1+ 1: и) = А(1+ 1: и, г' + 1: и) — А(1 + 1: и, 1) ь А(ю', з + 1: и), описывающей модификацию ранга 1 подматрицы А(1 + 1: п,1+ 1: п). Прочая арифметика в алгоритме, те.

А(з+1: п,1) = А(1+1: п,1)/А(з,1), выполняется на самом деле путем умножения вектора А(1+ 1: п,1) на число 1/А(1,1) по той причине, что умножение значительно быстрее деления; это также ВЬАЯ-операция, но уровня 1. Нам придется немного изменить алгоритм 2.4, чтобы воспользоваться им внутри версии уровня 3.

Алгоритм 2.9. Реализация ЬВ-разложения б. з выбора главных элементов с помощью уровня и пакета ВЬАЯ. Матрица А имеет размер т х и, где т > п; на ее место записываются т х п-матрица Ь и п х и-матрица У. Для последующих ссылок наиболее важные строки помечены номерами. /ос 1 = 1 со ппп(т — 1,п) (1) А(1+1:т,ю') =А(1+1:т,1)/А(1,1) '/'< (2) А(1+1: т,и+1: и) = А(1+1: т,г+1: п)— А(1+ 1: т, г) А(Ь з'+ 1: и) епд /о Левая часть рис. 2.б иллюстрирует алгоритм 2.9 в применении к квадратной матрице.

К началу ю'-го шага алгоритма уже вычислены столбцы матрицы 2.6. Блочные алгоритмы как средство повышения производительности 83 Швг! реализации ьц-рззлозсення с помощью уровня 2 пакета ВЕАВ Шзг! реализации ьС1-рвзлоясення о помощью уровня 3 пакета ВСАЗ Рнс. 2.6. Реализации 1Л1-разложения с помощью уровней 2 н 3 пакета В1 АБ. с — 1 Ь и — Ь вЂ” с+1 с — 1 А= Ь и — Ь вЂ” с+1 < Аы Агг Асз 'с Агс Агг Агз Азс Азг Азз с'ы Е'~г О'сз О Агг Агз О Азг Азз Т 1 О где все матрицы разбиты на блоки одинаковым образом.

Ситуацию иллюстри- руют правая часть рис. 2.6. Применяя алгоритм 2.9 к подматрице находим Агг Тгг с. ТггПгг Ь и строки матрицы 11 с номерами от 1 до с' — 1; на с-м шаге должны быть вычислены с-й столбец в Т и с-я строка в Ьг, а для оставшейся подматрицы в А должна быть проведена модификация ранга 1. Подматрицы на рисунке пронумерованы в соответствии со строками алгоритма ((1) или (2)), в которых они подвергаются пересчету. Модификация ранга 1 в строке (2) заключается в вычитании из подматрицы, помеченной номером (2), произведения закрашенных столбца и строки.

Указанные вычисления будут реорганизованы в алгоритме уровня 3 следующим образом: модификация подматрицы (2) будет отлозссена на Ь шагов, где Ь вЂ” малое число, называемое размером блока, после чего Ь модификаций ранга 1 будут произведены единовременно как одно матричное умножение. Чтобы понять, как это сделать, предположим, что первые с — 1 столбцов в Т и строк в У уже вычислены. Тогда 84 Глава 2.

Ретение линейных уравнений Это позволяет написать Агг Агз 122Угг Агз Вгг О 1 ~ 1122 Хгг'Агз Взг 1 ) ~ О Азз — Взг . (1 22'Агз) В32 г О .433 — Взг огз Вгг О П22 ~' 23 В результате получено модифицированное разложение, где в В вычислены до- полнительные Ь столбцов, а в У вЂ” дополнительные 6 строк: с Аы Агг Агз 1 ~ Ьы О 0 1 оы еггг Сгз Агг Агг Агз = Вм В» 0 . О Сгг Сгз Азг Азг Азз Взг Вгз 1 0 О Азз Все это определяет алгоритм, состоящий из трех шагов и иллюстрируемый правой частью рис. 2.б: (1) Посредством алгоритма 2.9 найти разложение ~ — ~ = ~ В сГ22. Агг Вгг Азг зг (2) Вычислить матрицу Угз = ВгггАгз.

Это равносильно решению треугольной системы линейных уравнений со многими правыми частями (Агз), что является одной операцией уровня 3. (3) Вычислить матрицу Азз = Азз — 132Угз, используя операцию матричного умножения. Более формальная запись алгоритма имеет следующий вид: Алгоритм 2.10. Ре лизация ЬУразложения без выбора главных элементов с помощью уровня Я пакета ВВАЯ. Матрица А имеет размер пхи; на ее место записываются В и У. Строки алгоритма нумеруются, как и в предыдущем его описании, и соответствуют частям правого рис.

Я.б. /ось=1 Соп — 1 эзер Ь (1) С помощью алгоритма з.9 найти разложение А(г: п,ь': 2+6 — 1) = '2 Угг Взг (2) А(г': 2+Ь вЂ” 1,2+Ь: и) = Вгг' А(г': г+Ь вЂ” 1,2'+Ь: и) /' вычислить матрицу Угз/ (3) А(2+6: п,2+ Ь; п) = А(2+ 6: п,г'+ Ь; и) — А(в+ 6: п,г: 2+Ь вЂ” 1) . А(г: 2+6 — 1,2+ 6: п) /' вычислить матрицу Азз '/ епЫ /ог 2.б. Блочные алгоритмы как средство пввьппвппя пропзводптельпостп 85 Чтобы получить максимальную производительность алгоритма, нужно еще выбрать размер блока Ь. С одной стороны, хотелось бы взять как можно большее Ь, поскольку, как мы видели, скорость перемножения матриц увеличивается с ростом их размеров.