XXI Зарубин B.C. Математическое моделирование в технике (1081441), страница 57
Текст из файла (страница 57)
Операция умножения матрицы на вектор составлена из серии скалярных умножений, и такую операцию следует отнести к более высокому уровню иерархии. Наконец, умножение матриц состоит из серии умножений матрицы на вектор. Перечисленные операции можно осуществить различными способами, и в зависимости от структуры векторов и матриц среди этих способов может быть найден наиболее рациональный с точки зрения производительности ЭВМ способ. Пример 7.1. Пусть необходимо вычислить вектор я = Ах, равный произведению матрицы А размера т х и и вектора х размера и. Казалось бы, удобнее каждую из строк этой матрицы с номером 1 = Т, пз скалярно умножить на вектор х и полученные скалярные произведения записать в качестве координат вектора я.
Для этого понадобится произвести тп умножений и тп сложений (если каждой из координат я; вектора я сначала присвоено нулевое значение). Но можно поступить иначе, если использовать операцию, состоящую из умножения вектора на число и сложения результата с другим вектором. Для этого при начальных нулевых значениях ьп и каждом фиксированном номере ~ = 1, и столбца матрицы А следует во внутреннем цикле по 1 = 1, т координату х вектора х умножить на элемент А;.
этой матрицы и результат сложить со значением г;, присвоив я; полученное значение суммы. Несложно проверить, что при этом количество умножений и сложений останется прежним. Различие состоит в том, что в первом случае доступ к элементам матрицы идет по строкам при фиксированных номерах г, а во втором — по столбцам при фиксированных номерах ~. Объем вычислительной работы, выполняемой ЭВМ-программой, измеряют количеством операций сложения или умно- 7.л. Вычислительные операции линейной алгебры 379 жения чисел с тиавающей точкой (или запятой), или галопов (сокращение английских слов Йоа$1пя ро1п1 орегаФюп). Таким образом, каждая из рассмотренных выше процедур умножения матрицы на вектор требует выполнения И = 2тп флоп.
Пример 7.2. Выясним, нельзя ли уменьшить число Ф операций в случае умножения симметрической матрицы А порядка н на вектор ю размера гь Ясно, что если не использовать симметрию матрицы, то Лс = 2п~. Можно попытаться уменьшить Ф за счет уменьшения объема памяти, требуемой для хранения симметрической матрицы. Ее можно хранить в одном векторе а размера п(п+ 1)/2 как, например, нижнюю треугольную матрицу [ПЦ, причем А17 — — а;+~„~узд~ 1р ь' = 1, чь, у = 1, ь'. Тогда для каждого фиксированного номера 1 = 1., и строки следует выполнить два внутренних цикла по номерам у' столбцов: сначала при ~ = г+1, и вычислить я; = я; + а +~„уз)б цх~, а затем при у' = 1,1 определить я; = я;+ а;+~„уд1 1)хл.
Но такой способ требует также 2пз флоп. Причиной неудачи является то, что алгоритм обращается к элементам симметрической матрицы не в порядке их расположения в массиве, в котором они хранятся, т.е. не в порядке возрастания номера координат вектора а. Расположим эти элементы в векторе а в порядке следования их по диагоналям, начиная с главной: при 1 >~' и й >О имеем А;~.ь; =а;+„ь ьр, 1уз, г'=1, и. Если теперь для каждого значения и = 1, и — 1, 1 = 1, и — к, вычислить сначала я; = я;+а1+сх;+у„а затем я;+ь = г;+я+ а1+са;, где ~ = п1с — к(к — 1)/2, то для нахождения всех координат искомого вектора л = Аж потребуется 2п(п — 1) флоп.
Таким образом, экономия по числу операций составит 2п флоп. Умножение матрицы А размера т х п на матрицу В размера и х г должно в общем случае потребовать 2тпг флоп. Если учитывать конкретную структуру матриц, то это число можно существенно сократить. Пример 7.3. Пусть необходимо перемножить две верхние треугольные матрицы А и В порядка и. Несложно убедиться, 380 7. АЛГОРИТМИЗАЦИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ что в результате получим также верхнюю треугольную матри цу С = АВ того же порядка.
Ее элементами будут укороченные скалярные произведения строк матрицы А на столбцы матрицы В с учетом того, что о;ь1аб = 0 при к < 1 или у < к, г, у', к = 1, и, Отсюда следует, что элементы матрицы С равны Таким образом, вычисление элемента с; при 1 < у' требует 2(у — г+ 1) флоп. В итоге, пренебрегая для больших значений и слагаемыми низших степеней, получаем* д1= ~ 2 2(у — ю'+1)=2~> ~~у 1=2~~) а=1 1=1 ~=1 1 1=1 из 2и из =') (и — г+1)'+ ~ (и — г+1) = — +из+ — = —. 3 3 3 ь=1 1=1 Итак, для умножения верхних треугольных (и аналогично нижних треугольных) матриц порядка и требуется примерно в 6 раз меньше арифметических операций по сравнению с умножением произвольных матриц того же порядка.
В вычислениях с матрицами больших размеров наряду с экономией арифметических операций актуальной является и экономия памяти ЭВМ, используемой для хранения исходных матриц и записи результатов вычислений. Вопрос использования памяти особенно важен, поскольку память современных ЭВМ организована по иерархическому принципу. Процессору, выполняющему арифметические операции, требуется наименьшее время для обращения к так называемой изида-иалтлтии (буферной памяти между регистрами центрального процессора и основной частью оперативной памяти). Несколько большее *Смл Голуб Дон., Вон Лоун Ч.
7.2. Вычислительные операции линейной алгебры 381 время требуется для записи и считывания данных в основной части оперативной памяти. Доступ к внешней памяти, имеющей большую емкость, требует еще больших затрат времени. Поэтому при разработке алгоритмов матричных вычислений следует стремиться к уменьшению обмена данными с внешней памятью, а данные, попавшие в кэш-память, целесообразно использовать максимальным образом. Пример Т.4. Пусть требуется вычислить произведение С = АВ квадратных матриц А и В порядка и, причем в кэшпамяти можно разместить лишь М чисел с плавающей точкой (М = Зп (( пз). Построим алгоритм вычислений так, чтобы каждый столбец матрицы В попадал в кэш-память только один раз.
Для этого при каждом фиксированном значении 1 = 1, и следует создать в кэш-памяти нулевой столбец С матрицы С и загрузить в кэш-память столбец В матрицы В, а затем для фиксированного значения й = 1, и загрузить столбец Аь матрицы А. Во внутреннем цикле по г =1, п проводим вычисление элемента сау = с; + а;ьбьуч Эти циклы повторяем для каждого значения й =1, и, так что в итоге получаем в кэш-пзмяти столбец С искомой матрицы С, который следует переписать в основную оперативную память, после чего перейти к следующему значению номера 2'. Подсчитывая количество чисел с плавающей точкой, проходящих при выполнении описанного алгоритма в обоих направлениях по каналу связи между кэш-пзмятью и основной оперативной памятью, получаем п п Зп4 д1 — "~ '« ~2п+ ~~» Ьп~ — 2пз+пз -2пз+ М 1=1 1=1 Отсюда видна важность увеличения кэш-памяти.
Предположим, что число М таково, что в кэш-памяти можно разместить по а столбцов матриц В и С и один столбец матрицы А, т.е. М - (2са+ 1)п. Тогда описанный выше алгоритм 382 7. АЛГОРИТМИЗАЦИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ можно использовать для вычисления матрицы С по блокам, каждый из которых состоит из а столбцов. Прн этом количество чисел, проходящих в обоих направлениях по каналу связи между кэш-памятью и основной оперативной памятью, будет равно 1Ч1 = ~~» 2иу'+ [и/а)* »» йп = 2п~+ип[п/а[*, где [и/а)* — наименьшее натуральное число, равное числу и/а или большее его.
Если — Е 1Ч, то 1»11 = 2п + — =2п + —. и 2»»к 2 2»»4 4» а М вЂ” и' При п » 1 имеем 1»«1 - -—, т.е. увеличение М в К» 2а+ 1 1»эз а' 3 позволит получить зкономию затрат времени в а раз. Покажем, что в данном случае кэш-память даст больший эффект, если при перемножении квадратных матриц А и В порядка п использовать их блочное представление, причем каждую из них представить состоящей из (и/8)2 блоков в виде квадратных подматриц порядка б (полагаем, что — Е Ы). Тогда, согласно правилу умножения блочных матриц [П1], блок С12, 1,,7 = 1, п/»э, матрицы С равен »»/Р С12 = ~~, А1КВК2, К=1 где А1к и Вкл — блоки матриц А и В соответственно.
Если для фиксированных значений 1, 1 = 1, и/~3 создать в кэш-памяти нулевой блок С1,1 и во внутреннем цикле по К = 1, п/11, вызывая в кэш-память блоки А1к и Вку, провести матричные вычисления См = С12 + А1кВкл, то в итоге получим блок С12 искомой матрицы С, который затем следует переписать в основную оперативную память. При реализации такого алгоритма по каналу связи между кзш-памятью и основной оперативной памятью в обоих направлениях пройдет количество чисел »»1»»»»Я2»»Я1 к = 2 2»4(4 «.
2 2кК ) = «2 1=12=1 К=1 Р 7.3. Алгоритмы векторно-конвейерных вычислений 383 При сравнении Д7г с Д71 следует учесть, что кэш-память объемом М = (2о+1)и позволит выбрать 11 из условия М = Зонг. Следовательно, при и >) 1 получим з 2иг+ и 1 а 3 и~/ — и 2а М вЂ” и ~/ЗММ и +2— У Так, при и = 21а и М = 2м по затратам времени на пересь|лку чисел из основной оперативной памяти в кэш-память и обратно последний алгоритм экономнее предыдущего в 4,9 раза.
7.3. Алгоритмы векторно-конвейерных вычислений *См:. Савельев А.Я. Скалярное умножение векторов — вычислительная операция, занимающая в силу своей простоты низший уровень иерархии вычислительных операций линейной алгебры (см. 7.2). Благодаря специальной организации вычислений можно добиться существенного повышения производительности ЭВМ при выполнении этой операции. Основной особенностью такой организации является так называемая конвейеризация выполнения операций, входящих в скалярное умножение векторов.