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

Блочные алгоритмы как средство повышения производительности 73 матрицы Р [244]. На практике можно взять и две диагональные матрицы Рсе н Р„1 с тем, чтобы решать систему [Р„мАР„1)В = Рсемб, л = Р,ыай Методы итерационного уточнения и уравновешивания реализованы соответственно ЬАРАСК-подпрограммами вйегг в и вйеейп. В свою очередь, к ним обращаются программы-драйверы типа вйевчх. 2.6.

Блочные алгоритмы как средство повышения производительности В конце раздела 2.3 было указано, что, в зависимости от используемого компьютера и решаемой задачи, изменение порядка трех вложенных циклов в алгоритме 2.2, реализующем гауссово исключение, может изменить скорость работы алгоритма на несколько порядков.

В этом разделе мы исследуем, почему это возможно, и опишем некоторые линейно-алгебраические программы, учитывающие такую возможность и написанные с особой тщательностью. Речь идет о так называемых блочных алгоритмах, самые внутренние циклы которых оперируют с квадратными или прямоугольными блоками матрицы, а не ее строками или столбцами. Соответствующие программы имеются в общедоступных библиотеках, таких, как ЬАРАСК [фортранная версия по адресу ХЕТЫВ/1арас)с)' и ЯсаЬАРАСК (см. ХЕТЫВ/эса)арас)с). ЬАРАСК и его версии на других языках предназначены для персональных компьютеров, рабочих станций, векторных компьютеров и параллельных машин с разделяемой памятью.

В их число входят Япп ЯРАНС-сепсег 2000 [238], ЯО1 Рожег С)за11епйе [223], ПЕС А1р)за Яегчег 8400 [103] и Сгау С90/390 [253, 254]. Для параллельных машин с распределенной памятью, таких, как 1ВМ ЯР-2 [256], 1п1е1 Рагайоп [257], серия Сгау ТЗ [255] и сети рабочих станций [9], разработана библиотека ЯсаЬАРАСК. Упомянутые библиотеки и подробные руководства для пользования ими доступны в )ХЕТЫВ'е [10, 34].

Более подробное обсуждение алгоритмов для высокопроизводительных [в особенности, параллельных) компьютеров можно найти в Интернете по адресу РАНАЬЬЕЬ НОМЕРАСЕ. Первоначальным стимулом разработки библиотеки ЬАРАСК была крайне неудовлетворительная производительность ее предшественниц, библиотек ЫХРАСК и Е1ЯРАСК (также доступных в ХЕТЫВ'е) на некоторых мощных компьютерах. Рассмотрим, например, приводимую ниже таблицу, показывающую скорость в мегафлопах Ы)ХРАСК-программы врота (реализующей метод Холесского) на машине Сгау УМР, суперкомпьютере конца 80-х годов. Метод Холесского — это вариант гауссова исключения, предназначенный для симметричных положительно определенных матриц.

Он подробно обсуждается в разд. 2.7; пока же нам достаточно знать, что он весьма схож с алгоритмом 2.2. В таблице указана также скорость нескольких других стандартных линейно- алгебраических операций. Сгау УМР— это параллельный компьютер, в котором одновременно могут работать до 8 процессоров. Поэтому мы приводим два столбца данных: для однопроцессорного и восьмипроцессорного режимов. 1 Имеется также С-версия 1.АРАСК'а, называемая СЬАРАСК 1см.

ХЕТЫВ/с1араск). ЬАРАСК++ (см. ХЕТЫВ/с++/1араск++) и ЬАРАСК90 1см. ХЕТЫВ/1араск90) являются соответственно С-ЬЧ вЂ” и Ротегаа90-интерфейсами для ЬАРАСК'а. Глава 2. Решение линейных уравнений 1 Проц. 8 Проц. Мвксимвльнвя скорость Умножение двух матриц (н = 500) Умножение матрицы на вектор Сп = 500) Решение уравнений ТХ = В (и = 500) Решение системы Тх = 5 (н = 500) Ь11с1РАСК (СЬо!евйу, и = 500) 1,АРАСК (СЬо!евку, и = 500) 1 АРАСК (СЬо1ея1су, н = 1000) 2640 2425 2285 2398 584 72 1414 2115 ЗЗО 312 311 309 272 72 290 301 Максимальная скорость компьютера, приведенная в верхней строке, является верхней границей для последующих чисел.

Данные по основным линейно- алгебраическим операциям в следующих четырех строках получены с помощью программ, специально разработанных для достижения высокой производительности на машине Сгау УМР. Эти данные достаточно близки к максимальной скорости компьютера; исключение составляет программа для задачи Тх = Ь (регпение единственной треугольной системы линейных уравнений), использующая наличие 8 процессоров незффективно. Задача ТХ = В— это решение треугольной системы со многими правыми частями ( — квадратная матрица).

Данные в таблице относятся к большим векторам и матрицам (и = 500). Программа метода Холесского из Ь1)лгРАСК'а (см. шестую строку таблицы) работает намного медленнее предыдущих программ при том, что применяется к столь же большой матрице и выполняет математически сходные операции. Подобная плохая производительность побудила нас модифицировать зту и другие линейно-алгебраические программы с тем, чтобы заставить их работать столь же быстро, как более простые программы типа матричного умножения. Скорости этих модифицированных программ из ЬАРАСК'а приведены в двух последних строках таблицы.

Очевидно, что ЬАРАСК-программы гораздо ближе к максимальной производительности компьютера. Подчеркнем, что программы метода Холесского из обеих библиотек, ЬГгт'РАСК и ЬАРАСК, выполняют одни и те же операции с плавающей точкой; различен лишь порядок операций. Чтобы понять, как было достигнуто зто ускорение, мы должны разобраться, на что тратится машинное время при выполнении программы.

Это, в свою очередь, требует понимания того, как функционирует память компьютера. Оказывается, что память любого компьютера, от самой дешевой персональной машины до самого мощного суперкомпьютера, устроена иерархически и состоит из ряда различных видов памяти с очень быстрой, дорогостоящей и потому малой памятью на вершине иерархии и медленной, дешевой и очень большой памятью в ее подножии. Быстрая, мелея, дорогая Регистры Кэш-пвмять Оперативная память Диски Ленты Медленная, большая, дешевая 2.6. Блочные алгоритмы как средство повышения производительности 75 К примеру, регистры являются наискорейшим видом памяти; затем следуют кэш-память, основна; память, диски и, наконец, ленты как самый медленный, больпюй и дешевый тип памяти.

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

Это означает, что при непродуманной реализации алгоритма ббльшая часть времени будет им тратиться не на реальную работу, а на перемещение данных из основания иерархической памяти в регистры. Вот пример простого алгоритма, в котором, к сожалению, нельзя избежать потери большей части времени на обмен данными, а не на полезную арифметику. Предположим, что нужно сложить две п х и-матрицы настолько большие, что они могут поместиться лишь в большом и медленном уровне иерархической памяти. Для сложения матрицы частями должны быть пересланы в регистры, где и производятся сложения элементов, а суммы должны быть пересланы в обратном направлении. Таким образом, на каждое выполняемое сложение приходится ровно 3 пересылки между быстрой и медленной памятью (считывание двух слагаемых в быструю память и запись одной суммы в медленную память). Пусть 1агнь (секунд) — время выполнения операции с плавающей точкой и 1, (секунд) — время пересылки одного слова данных из одного уровня памяти в другой, причем 1юе» банна.

Тогда время работы всего алгоритма составит пз(д,н,ь + 31, ), что намного больше времени пзгагнь, требуемого для арифметики самой по себе. Это означает, что сложение матриц обречено выполняться со скоростью самого медленного уровня памяти, хранящего матрицы, а не с гораздо более высокой скоростью сложения чисел. В противоположность этому, как мы увидим позднее, для некоторых других операций, например матричного умножения, можно добиться скорости, характерной для наиболее быстрого уровня памяти, даже если данные первоначально хранятся в наиболее медленном уровне. Программа метода Холесского из 1 1ХРАСК'а работает столь медленно потому, что в ней не предусмотрена минимизация движения данных между уровнями памяти на таких машинах, как Сгау УМР'.

Напротив, для матричного умножения и трех других основных алгоритмов линейной алгебры, указанных в таблице, такая минимизация была проделана. 2.6.1. Основные подпрограммы линейной алгебры Разрабатывать для всякого нового типа компьютеров специальную версию каждой программы, вроде программы метода Холесского — это явно не эффективное решение. Требуется более систематический подход к этой проблеме. Учитывая, что операции типа матричного умножения столь распространены 1 Хотя в ней была предусмотрена минимизация другого параметра обмена между основной памятью и диском, а именно страничных отказов (раде увила) 76 Глава 2. Решение линейных уравнений в вычислениях, производители компьютеров каталогизировали их как «Основные подпрограммы линейной алгебрыь (Вал«с Ьепеаг А1деога Яибгоийпег), или ВЬАЯ ]169, 89, 87], и оптимизировали их для своих машин. Другими словами, для высокопроизводительных компьютеров (как и ряда других машин) существуют библиотеки подпрограмм, включающие в себя умножение матриц, матрично-векторное умножение и другие подобные операции, со стандартным интерфейсом на фортране или С, однако тела этих подпрограмм оптимизированы для каждого типа компьютеров.

Наша цель — использовать этот оптимизированный пакет ВЬАЯ, модифицируя алгоритмы типа алгоритма Холесского так, чтобы они ббльшую часть работы выполняли с помощью ВЬАЯ- подпрограмм. В этом разделе мы дадим общее описание пакета ВЬАЯ. В разд. 2.6.2 будет рассмотрена задача оптимизации матричного умножения. Наконец, в разд. 2.6.3 мы покажем, как модифицировать гауссово исключение с тем, чтобы ббльшая часть его работы производилась посредством матричного умножения. Рассмотрим пакет ВЬАБ более подробно.