Дж. Деммель - Вычислительная линейная алгебра (1156793), страница 18
Текст из файла (страница 18)
В таблице 2.1 подсчитаны число обращений к памяти и число операций с плавающей точкой для трех родственных ВЬАБ-подпрограмм. Например, число обращений к памяти для операции вахру в верхней строке таблицы равно Зп + 1, поскольку требуется переслать из медленной памяти в регистры и значений км п значений у, и значение о, а затем записать и значений у; снова в медленную память. В последнем столбце дано значение (точнее, его старший член по и) отношения д числа флопов к числу обращений к памяти. Смысл коэффициента д в том, что он показывает приблизительное число флопов, приходящихся на одно обращение к памяти, иначе говоря, сколько полезной работы приходится на время, затрачиваемое на перемещение данных.
Эта характеристика показывает, насколько быстро алгоритм может работать в принципе. Предположим, например, что в алгоритме выполняется 7 операций с плавающей точкой с затратой 1„иь секунд на каждую и пъ обращений к памяти, каждое из которых требует 1, секунд. Тогда общее время работы алгоритма составит гп 1тет 1 этет 1 У 1«н«ь+гп'айшет = У'1«неь' 1+ ] =У'1«г!«ъ ' 11+ ] У 1енеь Ч 1енеь Мы предполагаем, что арифметика и обращения к памяти не производятся одновременно. Таким образом, чем больше д, тем ближе время работы алгоритма к наилучшему возможному значению 2" 1 „еь; так быстро алгоритм работал бы, если бы все данные находились в регистрах. Это означает, что алгоритмы с большими значениями д являются наилучшими строительными блоками для других алгоритмов. В таблице 2.1 отражена иерархия операций.
Операции типа эахру работают с векторами, производят 0(п ) флопов и имеют наихудшие значения д. Они составляют уровень 1 пакета ВЬАЯ, или ВЬАЯ1 [169]; к этому уровню относятся скалярное произведение векторов, умножение вектора на число и другие простые операции. Такие операции, как умножение матрицы на вектор, работают с матрицами и векторами, производят 0(п~) флопов и имеют несколько лучшие значения д. Они составляют уровень 2 пакета ВЬАБ, или ВЬАЯ2 ~89, 88]; он включает в себя решение треугольных систем уравнений и матричные 2.6.
Блочные алгоритмы как средство повышения производительности 77 Таблица 2.1. Подсчет числа операций с плавающей точкой (/) н числа обращений к памяти (т) для ВЬАБ-подпрограмм. модификации ранга 1 (т. е. операцию А + хут, где х и у — векторы-столбцы). Операции типа матричного умножения работают с парами матриц, производят 0(п ) флопов и имеют наилучшие значения д. Они составляют уровень 3 пакета ВЬАЯ, или ВЬАВЗ (87, 86]; к их числу относится решение треугольной системы уравнений со многими правыми частями. Директория ХЕТЬ1В/Ыаз содержит документацию и (неоптимизированные) реализации всех ВЬАЯ-подпрограмм.
Краткое описание всего пакета можно найти на странице ХЕТЬ1В/Ыаз/Ыазог.рз. Это описание дано также в [10, приложение С) (см., кроме того, ХЕТЬ1В/!арве!с/1пй/!арасй1п8.1пш1). Поскольку наивысшие значения д соответствуют уровню 3, мы попробуем модифицировать наши алгоритмы так, чтобы они опирались на операции типа матричного умножения, а не на матрично-векторное умножение или операцию вахру. (1,1ХРАСК-подпрограмма метода Холесского организована как последовательность обращений к процедуре вахру.) Подчеркнем, что модифицированные алгоритмы станут быстрее лишь при использовании оптимизированных ВЬЮ-подпрограмм. 2.6.2.
Как оптимизировать матричное умножение Исследуем подробно задачу реализации матричного умножения С = А. В+ С с целью минимизировать движение данных и тем самым повысить производительность. Мы увидим, что получаемая производительность сильно зависит от деталей реализации. Чтобы упростить обсуждение, используем следующую модель компьютера. Предполагается, что матрицы хранятся по столбцам, как это принято в фортране. (Приводимые ниже примеры легко модифицировать для случая, когда матрицы хранятся по строкам, как в С.) Считаем, что в иерархической памяти только два уровня, быстрый и медленный. При атом медленная память достаточно велика для того, чтобы в ней поместить три и х и-матрицы А, В и С. В то же время, быстрая память состоит лишь из М слов, где 2п ( М «пз; зто означает, что быстрая память может вместить в себя два столбца или две строки матрицы, но не матрицу целиком.
Мы предположим еще, что программист в состоянии управлять движением данных. (На 78 Глава 2. Решение линейных уравнений практике движение данных может управляться автоматическими устройствами типа контроллера кэш-памяти. Тем не менее основная схема оптимизации не меняется и в этом случае.) Простейший возможный алгоритм матричного умножения состоит из трех вложенных циклов, которые мы аннотировали, чтобы указать движения данных.
Алгоритм 2.6. Неблочная версия магаричного умножения (аннотированная с целью показать работу памяти) /ог1=1 еоп ( Счигпать сгпроку г' матрицгя А в быструю память ) /ог з = 1 со п ( Считать С; в быструю памлть ) ( Считать столбец у матрицы В в быструю память ) /ог й = 1 со п С;,=С,,+А;, В„ епб /ог ( Записать С; обратно в медленную память ) епб /ог епб /ог Самый внутренний цикл состоит в модификации элемента Сы скалярным произведением 1-й строки матрицы А и ~-го столбца матрицы В, что иллюстрируется следующим рисунком: Два внутренних цикла (по з и к) можно еще описать как модифицацию г-й строки матрнпы С посредством векторно-матричного произведения: г-я строка матрицы А умножается на матрицу В.
Это показывает, что мы не сможем получить производительность, более высокую, чем скорость операций из уровней ВЬАБ1 и ВЬАЯ2, поскольку именно такие операции используются во внутренних циклах алгоритма. Приведем подробный подсвет числа обращений к памяти: пз при и-кратном считывании матрицы В (по одному разу для каждого значения г), п~ при считывании матрицы А (считываем по очереди каждую строку и храним ее в быстрой памяти, пока она не становится далее ненужной) н 2пз обращений для матрицы С (каждый элемент СВ поочередно считывается в быструю память, хранится там, пока не закончится его модификация, а затем снова записывается в медленную память).
Всего получается пз + Зпз обращений, что дает д = 2пз/(пз+ Зпз) 2, т.е. значение, не лучшее, чем в ВЬАЯ2, и далекое от наилучшего возможного значения п/2 (см. табл. 2.1). Если М « и, так что в быстрой памяти нельзя разместить даже одну полную строку матрицы А, то 2.6. Блочные алгоритмы как средство повышения производительности 79 значение д уменьшается до 1, поскольку алгоритм превращается в многократное вычисление скалярных произведений, что соответствует уровню ВВАЯ1. Каждая перестановка трех циклов по г, у и и дает иной алгоритм, но все они имеют приблизительно одно и то же значение коэффициента у.
Мы предпочтем этим алгоритмам метод, использующий блочное разбиение: С рассматривается как блочная матрица размера Ю х Ф, где Сб — блоки порядка п/Л; таким же образом разбиваются на блоки матрицы А и В. Рисунок иллюстрирует это для случая )ч' = 4, указывая заодно модификацию алгоритма. Алгоритм 2.7. Блочная версил матричного умношсения /аннотированная с целью показать работу памяти) /ог з = 1 2о Л /ог у = 1 2о гч ( Считать Сб в быструю память ) /сей=1 зоЛ ( Считать Агь в быструю памятпь ) ( Считпать Вез в быструю память ) СБ = С'~ + А*"В"' еаза /ог ( Записать С'~ обратно в медленную память ) спи /ог епа /ог Число обращений к памяти подсчитывается следующим образом.
Каждый блок С'~ считывается и записывается обратно по одному разу, что дает 2пз обращений; матрица А в целом считывается Х раз (поскольку каждая п/Х х п/Ф-подматрица Ась считывается 1ч' раз), это дает Лпз обращений; столько же обращений требует матрица В (каждая подматрица В"' считывается Х раз).
Всего получается (2)ч'+ 2)пз ш 2Хпз обращений. Таким образом, чтобы минимизировать число обращений к памяти, нужно взять как можно меньшее Х. Однако Х подчиняется ограничению М > 3(п/И)з, означающему, что быстрая память должна одновременно вмещать в себя по одному блоку матриц А, В и С. Отсюда получаем Х п/3/М, а потому д ш (2пз)/(21ч'пз) ~/М/3, что значительно лучше, чем в предыдущем алгоритме. В частности, значение д растет вместе с М независимо от п; можно ожидать поэтому, что алгоритм эффективен для любого порядка п и его 80 Глава 2. Решение линейных уравнений эффективность растет с увеличением объема М быстрой памяти. Оба этих свойства весьма привлекательны. Можно показать, что, в действительности, алгоритм 2.7 асимптотически является оптимальным 115Ц.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
