автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.11, диссертация на тему:Решение систем линейных алгебраических уравнений с матрицами специальной структуры на векторно-конвейерной ЭВМ

российская академия наук ИНСТИТУТ ПРОБЛЕМ КИБЕРНЕТИКИ

На правах рукописи УДК 681.322:512.64

ЛОКОТЬ Елена Алексеевна

РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ С МАТРИЦАМИ СПЕЦИАЛЬНОЙ СТРУКТУРЫ НА ВЕКТОРНО-КОНВЕЙЕРНОЙ ЭВМ

Специальность 05.13.11 — математическое и программное обеспечение вычислительных машин и систем

автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Москва 1992

Работа выполнена в Институте проблем кибернетики Российской академии паук

Научный руководитель: доктор технических наук, профессор Ю. Г. ДАДАЕВ

Официальные оппоненты: доктор физико-математических наук Е. Е. ТЫРТЫШ НИ КОВ, кандидат физико-математических наук Н. Я. МАРЬЯШКИН

Ведущая организация: НИВЦ МГУ

Защита состоится « » ¿р1992 г. в ¿Г часов

на заседании специализированного совета К003.78.01 Института проблем кибернетики по адресу: 117333, Москва, ул. Вавилова, 37.

С диссертацией можно ознакомиться в библиотеке Института. Автореферат разослан «

» ^¡о^^/с^л-_ 1992 г.

Ученый секретарь А. 3. ИШМУХАМЕТОВ

специализированного совета

К003.78.01

Общая характеристика работы.

Актуальность. Одним из этапов решения таких задач, как численное решение дифференциальных уравнения, задач интерполяции и других, является решение систем линейных алгебраических уравнений Ах = г с лотрицали специальной структуры, то есть матрицами, ненулевые элементы которых расположены на небольшом (сгносительно порядка матрицы) числе диагоналей.

К таким матрицам относятся: трехдиагональные, пятидиагональные, блочно-трехдиагональные, ленточные и другие матрицы.

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

- каждое функциональное устройство может работать независимо от других и параллельно с ними;

- при выполнении различных операций использован принцип конвейера ;

- введены векторные регистры.

Реализация конвейера позволяет максимально загрузить функциональные устройства выполнением полезной работы, но для достижения этой цели необходимо организовать данные специальным образом, объединяя в отдельные потоки независимые однотипные операции.

Существование векторных регистров потребовало представления алгоритмов в виде последовательностей операций над векторами большой длины.

Для реализации на векторной архитектуре алгоритм необходимо представить в виде последовательности групп операций. • Все операции одной группы должны быть независимы и обладать возможностью быть выполненными одновременно на имеющихся в системе функциональных устройствах. Если алгоритм представим в описанном выше виде, то можно надеяться, что при его реализации на векторно-конвейерной ЭВМ будет достигнуто высокое быстродействие.

Можно выделить несколько способов построения эффективных векторных алгоритмов, а именно:

- модификация известных алгоритмов с целью их эффективной реализации на векторной архитектуре;

- создание новых алгоритмов, ориентированных на рассматриваемую

вычислительную архитектуру;

- комбинирование различных алгоритмов.

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

Это один из примеров того, что эффективность векторной реализации того или иного алгоритма во многом зависит от структуры матрицы.

Цель работы состоит в разработке и реализации пакетов ЛЕНТА и ИТЕРА, при помощи программ которых мокно эффективно решать системы линейных уравнений с матрицами специальной структуры на ВК ЭВМ прямыми и итерационными методами.

Научная новизна. При разработке пакета ЛЕНТА наряду с эффективной векторизацией известных векторных алгоритмов был эффективно реализован алгоритм конвейерной прогонки и разработан новый эффективный векторный алгоритм решения трехдиагональных систем линейных уравнений, алгоритм "деление на блоки".

Структура пакетов такова, что возможен выбор оптимальной программы,, решающей данную конкретную задачу.

Программы пакета ИТЕРА эффективно реализуют основные итерационные методы и ориентированы на возможное расширение класса итерационных методов. Программы пакетов взаимосвязаны, а именно, программы пакета ЛЕНТА используются при реализации итерационных методов.

Теоретическая и практическая ценность. Разработанные пакеты программ позволяют эффективно решать системы линейных алгебраических уравнений с :

- трехдиагональными матрицами;

- произвольными матрицами, имещими заполненную широкую ленту;

- пятидиагональными разреженными матрицами.

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

Программы пакета ИТЕРА можно использовать при расширении класса итерационных алгоритмов, основанных на расщеплении матрицы системы, при использовании попеременных и вложенных итераций.

Апробация работы. Основные результаты работы докладывались на семинарах отдела базового программного обеспечения Института Проблем кибернетики АН СССР, на международной информационной встрече по комплексной научной программе "Новые алгоритмы и архитектуры обработки гагформации" стран-членов СЭВ (г.Москва, 1989), на семинаре в Центре информатики и вычислительной техники Болгарской академии наук.

Публикации. Представленные на защиту результаты опубликованы в работах [1-4].

Структура диссертации. Диссертация состоит из введения, двух глав, заключения и списка литературы, содержащего 60 наименований.

Содержание работы.

Во введении обосновывается актуальность темы, определяется цель работы и формулируются ее основные положения. Приводится краткое содержание диссертационной работы.

В главе 1 представлены результаты разработки программ пакета ЛЕНТА, позволяющих эффективно решать прямыми методами системы линейных алгебраических уравнений с матрицами специальной структуры. Описаны алгоритмы решения трехдиагональных систем линейных уравнений. Рассмотрены возможности использования программ пакета для решения систем с произвольным числом ненулевых диагоналей прямыми методами.

В п.1.1 обосновывается необходимость пакета ЛЕНТА, описывается его структура.

Матрицы специальной структуры разделены на три класса (через 1 обозначим общее число ненулевых диагоналей, через »га - ширину ленты):

1. *а « 1 < 2п - ленточные матрицы| возникают при использовании разностных схем высокого порядка точности для решения обыкновенных

дифференциальных уравнений и дифференциальных уравнений в частных производных.

2. иа = 1 « п/2 ; например, трехдиагональные матрицы - возникают при использовании разностных схем второго порядка точности для решения обыкновенных дифференциальных уравнений; пятидиагональные -возникают при использовании разностных схем с более высоким порядком точности для решения обыкновенных дифференциальных уравнений.

3. иа » 1 ; их будем называть разреженными 1-диагональными матрицами; возникают при разностной аппроксимации дифференциальных уравнений в частных производных.

Разделение на такие классы продиктовано особенностями решения задачи Ах = г на Евкторно-конвейерной вычислительной архитектуре. Каждый класс порождает свою проблему выбора векторизуемого алгоритма и проблему хранения необходимых массивов в памяти.

Структура пакета ЛЕНТА такова, что для решения системы Ах = г с матрицей рассматриваемой специальной структуры можно либо найти в пакете оптимальную программу для решения системы с данной матрицей, либо использовать имеющиеся в пакете программы в качестве ядер программ пользователя.

Пакет ЛЕНТА включает в себя программы по следующим разделам:

1. Вычисления по схеме Горнера и схеме цепной дроби.

2. Решение трехдиагональных систем линейных уравнений с диагонально-доминантными матрицами.

3. Решение трехдиагональных систем линейных уравнений с произвольными невыровденными матрицами.

4. Умножение ленточных матриц на векторы и матрицы.

5. Перепаковка ленточных матриц.

6.' Отыскание треугольных разложений симметрических ленточных матриц.

7. Отыскание треугольных разложений ленточных матриц.

8. Решение систем линейных уравнений с треугольными ленточными матрицами.

9. Решение систем линейных уравнений с произвольными ленточными матрицами.

Программы пакета ЛЕНТА испо-г'^зуют ленточное хранение элементов ненулевых диагоналей матрицы А. Для удобства пользователя в раздел 5

включены программы, позволяющие переходить от одной формы хранения матрицы к другой.

Программы разделов 6-9 предназначены для решения систем с матрицами первого класса (ленточными). Степень параллелизма алгоритмов этого раздела wa и задача их реализации довольно проста.

В работе рассмотрены вопросы реализации на ВК ЭВМ алгоритмов решения систем линейных уравнений с матрицами второго и третьего классов. Основные проблемы, возникающие при реализации таких алгоритмов заключаются в следующем :

-степень векторизации основной части алгоритма исключения по Гауссу для матриц второго класса не превосходит величины max (wua.wla)i - при решении систем линейных уравнений с матрицами третьего класса требуется дополнительная память объемом n(wa-i), когда хранение самой матрицы занимает nl ячеек оперативной памяти.

Программы разделов 1-3, предназначенные для решения трехдиаго-нальных систем линейных уравнений, могут быть использованы и для решения систем с матрицами второго и третьего классов.

Программы раздела 4 позволяют решать задачи, требующие дополнительно матрично-векторного умножения.

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

В п.1.2.1 описаны основные алгоритмы, предназначенные для решения систем линейных уравнений Ах=г с трехдиагональными матрицами, допускающими LU-разложение. Представлено описание основных приемов реализации этих алгоритмов на ВК ЭВМ с целью получения эффективных вычислений.

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

В пакете ЛЕНТА реализовано четыре известных ранее алгоритма и два новых.

тизс - программа, реализующая векторизованный алгоритм прогонки:

1. Алгоритм ьи - разложения :

1;. = с. Ь. ., 1 = 2.......

I 11-1 '

а, = 1 / а,;

Д = 1 / (а,- 1 = 2.......

= с^.,« 1 = 2.....п;

ъ? = й.ъ., 1 = 1.....11-1.

VII

2. Алгоритм прямой прогонки : у, = V

= г.- 1 = 2.......

у" = у. Д., 1 » 1,...,п-1.

3. Алгоритм обратной прогонки : Хп - УпЙп;

Х1 - »Г - Х1ЧЬГ- 1 = П-1.....1.

Основная часть вычислений алгоритма может быть осуществлена только на скалярных функциональных устройствах. Время работы векторных функциональных устройств для обменов с памятью и вычисления 1;., 1 ^, ь", у*, 1 = 1,...,п занимает лишь девятую часть от общего времени вычисления:

75п - 67 + ЗТ(КЗ) (т+1), где Т(ВБ) + тах(5,УЬ) - время блокировки памяти при считывании, п = 64т + в.

• В результате работы модуля трбс получается ьи-разложение матрицы а, и при решении системы ах=в, трехдиагональной системы с несколькими (к) правыми частями, при достаточно больших к, доля скалярных вычислений от общего времени 35п + 5пк будет невелика.

тисы - программа, реализующая метод векторизованной встречной прогонки.

Основная часть вычислений, как и для алгоритма прогонки, может быть осуществлена только на скалярных функциональных устройствах. Но так как процедура прогонки применяется одновременно к двум системам порядка п/2 (сверху вниз и снизу вверх), степень параллелизма алгоритма встречной прогонки в 2 раза больше, то при его реализации на ВК

!

ЭВМ использовался конвейерный принцип работы скалярных функциональных устройств. Векторные вычисления программы trcn занимают шестую часть общего времени вычислений :

49п + [1 ОТ (RS ) + 21 (WS) + T(S)](1+1) - 44, где T(ws) + max(5,VL) - время блокировки канала для записи в память,

Т(S) - время зацепления с командой считывания из памяти,

641 + в = (п-2)/2.

trpl - программа, впервые реализующая алгоритм конвейерной прогонки, предложенный в работе М.А.Фрумкина "Параллельная и конвейерная сложность некоторых алгоритмов линейной алгебры" (препринт ИСК АН СССР, 1903). Алгоритм основан на конвейеризации основных циклов прогонки. Для этого используются вычисления по схеме цепной дроби и по схемам Горнера, основанные на применении формул интерполяции для вычисления значений линейной и дробно-линейных функций, подробно описанных в пп. 1.3.2, 1-3.3.

Для использования формул интерполяции узлы интерполяции вычисляются перед каждым этапом алгоритма по формуле

z* = z, z2 * z + 1, z3 = z 4- 2, 1 » 1,...,к, i i i

где z - параметр для текущих вычислений по схемам Горнера и цепной

дроби.

Время выполнения скалярных вычислений - ЮЗ(к-1) тактов. При больших п это составляет незначительную часть от общего времени вычислений по схеме конвейерной прогонки :

т + 11 з (1с - 1) + Юк, где Юк тактов необходимо на формирование узлов интерполяции, Т » 20п + 1110Т(RS) + 3T(S) + 3T(WS)] + 13k + 7T(RS) + T(S) + T(+), T(+) - время зацепления с операцией сложения.

Описанные выше три алгоритма, реализованные в пакете ЛЕНТА, являются модификациями метода прогонки, причем последний из них значительно выигрыввет по времени реализации на векторно-конвейерной ЭВМ при больших порядках матрицы. Первые эффективны для небольших п.

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

В пакете ЛЕНТА представлены три модуля, реализупцие алгоритм нечетно-четной редукции с обратным ходом (схеме без деления), алгоритм

Уонга и новый алгоритм "деление на блоки".

Общее 8тих алгоритмов состоит в редуцировании системы, то есть сведении ее к трехдиагонвльной системе меньшего порядка.

тйжюв - программа, реализупцая алгоритм нечетно-четной редукции с обратным ходом.

Идея алгоритма состоит в том, что путем элементарных преобразований из системы исключаются все нечетные переменные. В новой системе вновь исключаются все нечетные переменные и так до тех пор, пока не останется одно уравнение.

Модифицированное уравнение трехдиагональной системы можно записать в виде : .

" <СЯ1-4С«1/ а21-.>Х21-1 + <*«Г Ь11С«1/в«1-Г Ь21°«и/в2и.)Х2Г

- <ьгЛи1/а2и1)х2и1 = Г2Г (С21/а21-1)Г21-Г <ь«1/а«и1>р.и»' где 1 = 1,2,..., [т/2 ], т = п, [п/2] , ,..., 1 •

Когда все четные переменные найдены, нечетные определяются по формулам :

. = (г„ - с„х„ - Ъ„. ,х„. )/а„. ,, 1 = 1.....Гп/21.

21-1 21 21 21-2 21-1 21 21-1 1 1

Каждый шаг алгоритма нечетно-четной редукции векторизуется (длина обрабатываемого вектора [т/2]). Алгоритм реализован на ВК ЭВМ не только для тривиального случая, когда порядок системы п = 2к- 1, но и для произвольного п. 1

При реализации прямого хода, в случае отсутствия у последнего четного уравнения системы последующего нечетного уравнения, формально добавляется уравнение

I О » х. + х. + 0 » х. , = О, к-1 к к*1

что не влияет на результат вычисления и достигается путем уточнения последних компонентов векторов данных при их считывании из памяти.

При реализации обратного хода на каадом шаге удваивается количество неизвестных, которые могли бы быть вычислены на предыдущем шаге, и проверяется, сколько из них не выходит за границы массива длины п, отведенного для записи результата.

При реализации алгоритма скалярные функциональные устройства не использовались.

Время работы программы тшшое равно

23(п-1 ) + (ЗТ(В) + 15Т(Ш5) + ЗТ(«В) + 8] Цов^ +1-1), где п/2 = 641 + в.

Алгоритм нечетно-четной редукции представляет собой параллельный процесс обработки векторов, длина которых понижается в два раза на каждом шаге алгоритма. Известно, что обработка векторов небольшой длины (VI < 32) не эффективна. Поэтому, чтобы к концу нечетно-четного исключения не произошло торможения процесса вычислений, разумно в некоторый момент перейти на вычисления по схеме векторизованной прогонки. Тестирование показало, что такой переход наиболее эффективен при длине вектора к* = 28.

тиоизс - программа, реализующая такую модификацию нечетно-четной

редукции в пакете ЛЕНТА. Время работы программы -• •

(23 - 2 23п/(2 + 1)] + к"(72 - + 0г,

где т - время работы программы тшиюЕ,

О = ЗТ (Б) + 15Т(КЗ) + ЭТ№) + 8,

С2 = Т(>) - 67, Т(>) - время зацепления с операцией сдвига.

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

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

Программа тгшсм реализует известный алгоритм Уонга. После разбиения основной системы на подсистемы элементарные преобразования, приводящие матрицу к диагональному виду, применяются ко всем к блокам одновременно. Отыскание компонентов решения далее осуществляется векторным делением. Редуцированная система имеет порядок к.

Алгоритм хорошо векторизуется и длина обрабатываемых векторов равна числу блоков - к. Решение редуцированной системы реализуется на скалярных функциональных устройствах.

В программе ТШ)0Ю предполагается, что к < 64. Это упрощает реализацию алгоритма. Время работы программы ТМ)СХС -

2бп + б4к ♦ [8Т(КЗ) + 8Т(*5) + 7Т(Б) + Т(/) + Т(») + Т(+)]р + + 4Тда) - 5Т(*3) - 7Т(8) - Т(/) - Т(») - Т(+),

"деление на блоки", с учетом возможнос-

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

Программа тевьо реализует новый алгоритм который был разработан при создании пакета ЛЕНТА тей вычислительной архитектуры ВК ЭВМ.

Решение сводится к одновременному решению к трехдиагональных систем линейных уравнений и редуцированной трехдиагональной системы порядка 2к.

к

Если через а обозначить блочную матрицу а - • а., составлен-

1 = 1

ную из диагональных блоков разбиения матрицы а, а через в - матрицу в = а - а, то решение системы находится путем последовательного решения следующих систем :

А X =-г, (Е + В / а) I = В х, ¡1 = 1, х — х — х.

Так как матрица (Е + в / а) имеет вид (символами "*" обозначены места расположения ненулевых элементов)

11 1 0 * * -

0 11...

0 * * ... 0

0

*

0

1 »*...»

* >1

1

то система (Е - в/А)! = вх сводится к трехдиагональной системе линейных уравнений « I = « порядка 2к, где

1 ~ (íl^íp^ip*^•1гp^'••

,f<k-l>p♦l,fkp', 8

я — [я, . «я ,71 _ , ■ v? , • • • , w.. , , , w1 ) •

1' р+1* р гр+»" 2р (к-1>р+1 кр

Остальные компоненты векторов (щ равны о. Матрица М состоит из элементов, расположенных на пересечении 1р,(1 - 1)р + 1, 1 = 1,...,к, строк и столбцов матрицы (Е + в А"1).

Полностью алгоритм представим в следующем виде :

1. Решая систему А х = г, находим ьи-разложение матрицы А, векторы х и в х = 1».

2. Осуществляя прямую и обратную прогонки для решения систем Ас х = е( и А. х = ер, находим коэффициенты матрицы К.

3. Решая систему и 1 = находим вектор 1.

4. Осуществляя прямую и обратную прогонки для решения системы А х = г, находим вектор х.

5. Решение системы находится вычитанием х = х - х.

Учитывая, что для решения системы Аьх = ер не требуется прямая

прогонка, а результаты прямых прогонок систем А.х = е , А^х =

1 = 1,...,к отличаются лишь на множитель прямая прогонка для системы А х = г может не выполняться, а обратная прогонка -осуществляться по формулам

хр<1+«>- *р<1«1> ^ ир<1»1) + арс 1*1 >^р<1*1)• ^ ~

х. » г" . Г. . - х . Ь* , З = р-1......1,

1р*| 1Р+.1 м г •

где - вектор -результат прямой прогонки для системы А1х = в4:

Для повышения эффективности реализации алгоритма использовалась перегруппировка вычислений первого и второго шагов алгоритма, в результате которой стало возможым максимально использовать промежуточные результаты вычислений и тем самым сократить количество обменов с памятью.

1. ьи-разложение матрицы а, прямая прогонка для систем а х = г, А1Мх = е4, 1=0,...,к-1, и обратная прогонка для системы А1Мх = ер.

2. Обратная прогонка для систем Ах = г И А.х = е , 1=0,...,к-1.

3. Решете системы М Г = «г.

4. Обратная прогонка для системы А х = 1 и нахождение решения системы х = х - х.

Реализация первого шага осуществляется путем сквозной векторизации всех вычислений. Для получения ьи-разложения используется алгоритм прогонки с минимальным числом операций - 9п.

Узким местом алгоритма является решение системы Мх = * порядка 2(к - 1) < 126. Для решения системы мг = я можно выбрать один из рас смотрены! выше алгоритмов, однако, как показывает тестирование, для матриц порядка < 128 наиболее эффективными являются методы векторизованной и встречной прогонки. При реализации использовался алгоритм векторизованной прогонки. Для этих вычислений не требуется обменов с памятью. Элементы вектора правой части, результаты промежуточных вычислений и вектор-решение системы хранятся на векторных регистрах.

Общее время выполнения программы ТРЛЮ -

1бп + потде) + '5Т(»Б) + Т( + )]р + 71к, при этом время выполнения скалярных операций - б8к.

Программа тгакзь реализует одновременное решение к трехдиагональ-ных систем линейных уравнений путем сквозной векторизации алгоритма прогонки с минимальным числом операций одновременно к к системам.

Время работы программы : (9к + 1 )п + [бКйЗ) + ЗТ№)] (т + 1) тактов, где к = 64т + в.

В таблице 1 приведены основные характеристики предложенных в пакете программ.

Таблица 1.

Алгоритм Программа Условия для реализации алгоритмов Время реализации операций

последовательных векторных

Векторная прогонка ТЙЗС 67 (п-1) 8п

Встречная прогонка тнси 39(п-2) 10п

Конвейерная прогонка тюч, кр=п-1 ,к<6 задание 7к узлов интерполяции 113к 20п+1Зк+ +1(10Т(НЗ)+ +ЗТ(Б)+ЗТ№))

Нечетно- • четная редукция ТШФОЕ дополнительная память длины П - ■ ' ' 23(п-1)+ + (3(Т(Б)+Т(№5))+ +15Т(КЗ)+8)1о82п

Модифицированная нечетно-четная редукция ТНОЕЗО дополнительная память. длины п 67к* 23(п-к*-2) + + (3(Т(3)+Т(ИЗ))+ +15Т(НЗ)+8)* »Юв^п/к")

Метод Уонга ТМ)С1С кр=п, к<б4 72(к-1) 2бп-8к+(8Т(ИЭ )+ +8Т (УКЗ+7Т (Б)+ +Т(/)+Т(»)+ +Т(+))Р

Метод "деление на блоки" ТЕВЬО .кр=п, к<б4 б8к 1бп+р(12Т(РЗ)+ +ЗТ(ЯЗ)+Т(+))

В п.1.2.2 описывается реализация программ раздела 3 пакета ЛЕНТА, позволяющих решать трехдиагональные системы линейных уравнений с произвольными невырожденными матрицами. В пакете представлены программы тгшсое и тшюе, реализующие алгоритмы прогонки с выбором ведущего элемента и ортогональной прогонки, соответственно.

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

Алгоритм ортогональной прогонки представлен в пакете ЛЕНТА модулем тшхж. В программе реализована модификация алгоритма, которая не использует операции извлечения квадратного корня.

В п.1.3 описываются программы раздела 1 пакета ЛЕНТА, позволяющие проводить вычисления по схемам Горнера и цепной дроби, ап + Ьп(ап_1 + ... + Ьг(а4 + Ь^) ...) и

ап + Ьп/(ап1 + ... + (а1 + ^А) ...), соответственно (1; - параметр схемы).

Выполнение вычислений по схеме Горнера требуется, например, при решении треугольных систем линейных уравнений. Схема цепной дроби обеспечивает получение элементов матриц ь и и ьи-разложения матрицы А.

Если получено ьи-разложение матрицы А, то решение системы АХ = И с несколькими (к) правыми частями представляет собой вычисление по к схемам Горнера для осуществления прямой и обратной прогонок. Аналогичных вычислений требует и нахождение матрицы, обратной матрицы а.

Вычисление по одной схеме Горнера (схеме цепной дроби) реализовано в пакете ЛЕНТА модулем БИОВЗ. Выбор схемы определяется параметром р. Если р = 0, то вычисления проводятся по схеме Горнера, а если р = 1, то - по схеме цепной дроби. Основной цикл вычислений модуля шовэ является скалярным и реализует строго последовательные вычисления. Векторные регистры используются лишь для обмена с памятью.

Из всего многообразия вариантов распараллеливания вычислений по схеме Горнера и схеме цепной дроби для реализации в пакете ЛЕНТА выбран конвейерный алгоритм, который впервые реализуется в пакете ЛЕНТА модулями corner и drob, СООТВвТСТВвШЮ.

Алгоритм основан на применении формул интерполяции для вычисления значений линейной и дробно-линейной функций. Если требуется вычислить, например, выражение

ап+ W« + ••• + Vai + V> •••>•

то, выбрав произвольное целое число к, 1 < к < п, и 1 так, чтобы kl > n, k (1 - 1) < п, положим для 1 = 1,...,к

- au + bu<aw-. + •••<aiu-„M + Ь1и-,>..4> ■••>• Тогда у(t) = ylk(yUk.t)-•. (yt (t)...). Поскольку yu(t) - линейные функции, то достаточно иметь их значения уu(ti') и Уи(*") для различных значений t/и t" аргумента и воспользоваться интерполяционной формулой

у, (t) = (у, (t')(t - t") - у, (t")(t - t'))/(t" - t*).

"U 'Ii i i Ii i i' i i

Вычисление значений yu(t[) и yu(t") могут выполняться одновременно для всех 1=1,...,к, используя векторные операции с векторами длины к. На векторных функциональных устройствах можно провести и вычисление величин i/(t" - t'). Выбор к < 64 исключает нарезку векторов и, следовательно, обмен с памятью результатами промежуточных вычислений. Обмен с памятью необходим только для чтения данных и записи результата.

Аналогичные рассуждения применяются и для конвейеризации вычислений по схеме цепной дроби. Вычисления по интерполяционной формуле y(t) = ((t - - tpcyctp - y(t;))y(t-'; -

- (t - t^'xt^ - tj)(y(tp - y(t-'))y(t;)) /

<(t - t•)(17 - t-')(y(t") - y(t;>> -.- (t - - t;)(y(t-) - y(t|"))

требуют трех значений - yu(tl')> и yu~ для ТР9Х Раз"

личных значений аргумента.

Сравнение времен реализации, например, corner и drobs (1=0) показывает, что если n, к и 1 таковы, что 1 < Ск / С19k - т(s)), где С = Т(накл.), то применение программы corner менее эф1ективно. чем применение программы drobs (i=o). При 1 < Ck/(36n-T(s)) применение программы drob менее эффективно, чем применение программы drobs(1=1). Кроме того, в процессе работы программы drobg получаются все

промежуточные значения y^t) = а.+ Ь1У1_,(<;)« i=1.....в то время

как конвейерные алгоритмы вычисляют лишь к таких значений. Для получения же остальных значений yt(t), i=Jl, j=i,...,k, требуется дополнительное вычисление по формулам

у, . (t) = b,. у,. ,(t) + а,. - для схемы Горнера и

U»j lltj'U+j-i ll»j г г

yu+j(t) = bu*/yiUj-i(t) + alUj - схемы цепной дроби

(1 = 1.....k-1, i = 1.....1).

Эти вычисления могут быть выполнены векторными функциональными устройствами, но потребуют еще по крайней мере 2kl + T(s) и 3kl + T(S) дополнительного времени. Это еще больше увеличивает диапазон эффективности применения последовательных схем.

Вычисления по нескольким схемам Горнера (цепной дроби) реализованы в пакете ЛЕНТА модулем gorv. Выбор схемы, как и для программы DROBS, определяется параметром р. Вычисления проводятся по нескольким схемам одновременно, используя векторные функциональные устройства и вектора длины 1. Если l>64, то осуществляется нарезка. Время работы : Э1п + T(s)(s + 1)n + 1 + T(s) при p = о и 4In + T(S) (в + 1 )n + 1 + T(S) при p = 1.

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

Время работы : 22n + 4T(S) - для схемы Горнера и

40n + 4T(S) - для схемы цепной дроби.

В п.1.4 приводятся примеры использования программ пакета ЛЕНТА для решения систем линейных уравнений с произвольной матрицей специальной структуры. Например, для решения треугольных 1-дигональных разреженных систем можно использовать программы вычисления по схеме Горнера, проводя соответствующие преобразования правой части системы.

Рассматривается пример решения систем линейных уравнений с разреженными l-диагональными матрицами.

Пусть матрица А - 1-диагонально разреженная матрица, причем

-1; = 1, i<k<it. i^1 - = 1, i<e<iz.

Через а обозначим ленточную матрицу с шириной ленты 1* + 1* + 1,

1|-ые, 1=1.....к, 1^-ые, j=ip..,s и главная диагонали которой совпадают с соответствующими диагоналями матрицы А. Через в обозначим матрицу А - к.

Решение основано на сведении решения 1-диагонально разреженной системы линейных уравнений к решению системы с матрицей, имеющей заполненную ленту ( Ä ) и заполненной матрице меньшего порядка

min (n, 2п - Ц^1 + I®*1)). Ясно, что при + 1®+1 < п использование такого алгоритма приво-

дит к более сложным вычислениям, чем LU- разложение, и потребует еще большей дополнительной памяти. Если lj+1+l**'>n, то, используя рассуждения, аналогичные проведенным для метода "деление на блоки", находим, что х = х - х, где Ах = г, Äx = f, (Е + A_<B)i = Вх, а матрица (Е + ¿"'в) имеет ненулевые элементы только на пересечении n - lj*1 первых, п - последних строк и столбцов и только n-lj*4 первые и п-1®*4 последние компоненты вектора f ненулевые.

Решение системы будет состоять из следующих этапов.

1. Находим LU-разложение матрицы к. Решая систему Äx = г, находим вектор х. Находим вектор Вх.

2. Используя найденное LU-разложение, вычисляем элементы n-lj*4 первых и n - последних столбцов матрицы А"1. Находим матрицу Е + А~'в.

3. Решая систему порядка 2n - n - I**1), находим вектор Г.

4. Используя полученное LU-разложение матрицы Ä, находим вектор х как решение системы кх = t.

5. Решение системы находится путем вычитания х = х - х.

Шаги 1 и 4 алгоритма можно выполнить, используя программы пакета ЛЕНТА. Необходимые для дальнейших вычислений столбцы матрицы Ä"1 находятся путем использования программ для решения систем с несколькими правыми частями. Шаг 3 представляет собой решение заполненной системы порядка 2n - lj*'). Если этот порядок больше ю, то алгоритм

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

1) T(LU) - времени, необходимого для нахождения LU-разложения заполненной ленточной матрицы;

2) [(n - lJM) + (п - 1^') + 2] Т(Ly=r), где Т(Ly=r) - время, необходимое для решения нижнетреугольной ленточной системы;

3) I(n - lj*') + (n - ljj*1) + 2] T(Ux=y), где T(Ux=y) - время.

необходимое для решения верхнетреугольной ленточной системы;

4) т((Е + А"'в)г = и) - времени, необходимого для решения заполненной системы порядка 2п - (1км+ 1®*').

Если п_1г и t n_1i' то время л"11 нахождения

решения такой системы можно оценить следующим образом :

0(Т(ьи) + 2Т(Ly=r) + 2T(Ux=Y) + (n - 1) + n + const). Применение же стандартного алгоритма LU-разложения потребует времени о(пэ + 2п* + п). Если же l^'-lj 1 и -» 1, то эффективней будет применение ьи-разложения.

Как пример использования описанного выше метода представлен алгоритм решения систем с циклической трехдиагональной матрицей.

Приводится несколько примеров использования программ пакета ЛЕНТА для решения систем, имеющих матрицу с узкой заполненной лентой.

В главе 2 описывается реализация итерационных методов для решения систем линейных уравнений с матрицей специальной структуры.

В п.2.1 обосновывается необходимость пакета ИТЕРА1, кратко описывается его структура.

В п.2.2 описывается реализация наиболее используемых итераций для решения трехдиагональных, заполненных ленточных и пятидиагональ-ных разреженных систем линейных уравнений, которые включены в пакет ИТЕРА 1.

Метод Якоби представлен двумя программами :

- bmjac - для решения трехдиагональных и произвольных ленточных систем линейных уравнений;

- BHJAC5 - для решения пятидиагональных разреженных систем.

Запись в матричном виде итерации Якоби

xk = D"1 (г - Ахк_1) + хк_1 позволяет полностью векторизовать метод. Умножение матрицы на вектор в программе bmjac осуществляется вызовом соответствующей программы пакета лента. в программе bmjac5 умножение матрицы на вектор реализовано непосредственно, используя возможности зацепления с остальными операциями итерационного шага.

Время работы в тактах для выполнения одной итерации :

- программы BMJAC - 4n + TQ, где т - время работы соответствующей программы умножения матрицы на вектор;

- программы BMJAC5 - 1бп - 71t - 412.

Вычисление итераций Гаусса-Зейделя представлено в пакете ИТЕРА1 программами :

BMCSK5 - для решения пятидиагональных разреженных систем линейных уравнений;

trgasb - для решения трехдиагональных систем линейных уравнений; bmgask - для решения произвольных ленточных систем линейных уравнений.

Для эффективной реализации итераций Гаусса-Зейделя хк = (г - £ a xk - Е a xk"'). i = 1,....n

i I I J 1 Г" I j J

)<i ' ' j >i ' '

на векторно-конвейерной ЭВМ алгоритм был представлен в следующем виде:

xk,t - xk + rk/a - Еа " хк)/а .

i i i и j<t ч j j 11

Теперь хк*' momio представить как хк" = хк" + хк*' , где вычис-

1 ' V * i " V v I 1 1 1

ление вектора х = (х ) выполняется в векторном режиме, так как = xk-' + D~'rk.

Вычисление же вектора хк" = (хкм) для каждого из рассматриваемых типов матриц имеет свои особенности. В программе bugasf. эти вычисления можно провести на векторных функциональных устройствах. Для трехдиагональных матриц эти вычисления представлены скалярным циклом, в для пятидиагональных разреженных матриц - двумя скалярными циклами.

Время работы в тактах для вычисления одной итерации :

- программы BUGASE - 6n + 41zn + TQ, где TQ - время работы программы BMVBG пакета ЛЕНТА для умножения ленточной матрицы на вектор;

- программы TRGASE - 29п;

- программы BMGSE5 - АЬп - Ыг- 21 .

Метод Гаусса-Зейделя с последовательной верхней релаксацией в пакете ИТЕРА1 представлен тремя программами :

trsor - для решения трехдиагональных систем линейных уравнений; bhsor - для решения ленточных систем линейных уравнений; BMSOR5 - для решения пятидиагональных разрежешшх систем линейных уравнений.

Идея алгоритма состоит в том, что если найдено k-oe приближение, то сначала вычисляется хк*' по формуле приближений Гаусса-Зейделя, а затем

хкм = хк + и и"*1 - хк).

J J ^ Г

Если учесть преобразование, позволяющее частично векторизовать итерационный процесс Гаусса-Зейделя, то получим :

х"р* = хк + иКО-'СЬ - Ахк)^ + £*♦*], где хк" = -Е в)1(х^1-

Время работы программы, реализующей указанный алгоритм, совпадает с временем работы программы, реализующей алгоритм Гаусса-Зейделя.

Это достигается за счет вычисления вектора (ш/а ) в зацеплении с

11 к

вычислением вектора (1/аи), 1 = 1.....п и вычисления Д^ х" - х"

в неявном виде :

Д. = ш ((В"'(Ь - Ахк))1 + хкм).

Таким образом, не теряя во времени для вычисления одной итерации, можно ускорить сходимость итераций Гаусса-Зейделя.

Программа веса пакета ИТЕРА1 реализует метод сопряженных градиентов для решения систем линейных уравнений с трехдиагональными и ленточными матрицами, который представляет собой алгоритм минимизации квадратичной формы 1/2(хтах) - хть последовательно по системе независимых направлений рк и корректен только для симметрических положительно определенных матриц. При этом итерации сходятся к точному решению не более, чем за п шагов.

Для вычисления произведения Арк используются программы пакета ЛЕНТА. Для реализации остальных вычислений используются стандартные модели реализации вычисления триад и скалярных произведений.

Время выполнения одной итерации : ЗОп + Т0, где т - время работы соответствущей программы умножения матрицы на вектор.

Чтобы расширить класс матриц систем линейных уравнений, к которым применим метод сопряженных градиентов, используется метод сопряженных градиентов с предварительной симметризацией матрицы. Этот метод основан на том факте, что для любой невырожденной матрицы а произведение м = ата является симметрической положительно' определенной матрицей. Исходя из этого, метод сопряженных градиентов применяется к системе мх = ать. В пакете ИТЕРА1 этот метод реализован программой высисв. Вычисление ата и атъ путем вызова соответствующих программ пакета ЛЕНТА.

В случае плохо обусловленной матрицы А эффективным является использование метода сопряженных градиентов с предобусловливанием, пре-

дставленного в пакете ИТЕРА программой bmpcq, a именно, метода сопряженных градиентов, примененного к системе мах = иь, где м = (sTs)"', S - матрица-предобусловливатель.

При реализации метода использовались те же приемы, что и при реализации метода сопряженных градиентов. Реализация решения системы U rktl = гк+* осуществляется путем использования программы решения системы с матрицами типа и. Выбор такой программы определяется структурой матрицы 1(, требуемой точностью и эффективностью вычислений. Это могут быть программы как пакета ЛЕНТА, так и пакета ИТЕРА.

Проверка сходимости процесса осуществляется по схеме :

I (rkM,rkM)| < Е.

Если "да", то проверить: |гк,,| < е. Если "да", то х*= хкм, иначе продолжить вычисления.

Время работы программы bmpcg : t(cg) + t(sm), где т(СО) = ЭОп + TQ, TQ - время работы соответствующей программы умножения матрицы на вектор, t(SM) - время работы программы, реализующей решение системы lírkM = rktl.

В п.2.2.7 приведены сведения о сходимости и скорости сходимости представленных в пакете итерационных методов для матриц специальной структуры.

В п.2.2.8 приведены результаты сравнения программ пакета.

Применение пакета ИТЕРА1 ограничено не только выбором итерационного метода решения системы Ах= i""™« А. Если матрица А не является трехдиагональной, плотной леь.^ ...„/ " идиагональ-ной разреженной матрицей, то либо необходимо вносить изменения в представленные в пакете ИТЕРА1 программы, использовать программы пакета, проигрывая в эффективности (решать девятидиагональную разреженную систему как систему с заполненной лентой). Для использования какого-либо другого итерационного метода требуется написание новой программы.

В п.2.3 приводятся обоснование и структура пакете ИТЕРА.

В пакет включена программа, которая делает пакет более универсальным и гибким.

Используется тот факт, что любой иторпциошшй процесс, в основа которого лежит расщепление матрицы а, представим к виде

Мхи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*1 = 14 а(г - Ах1), и программа, позволяющая решить систему с матрицей к (Мх1+1= 11+|), задающей некоторое расщепление матрицы А = м - ь, то обобщенный итерационный процесс можно реализовать по схеме :

Первый блок вычислений может быть опущен, если положить х° = 0. Для выполнения второго блока вычислений могут быть использованы программы пакета ЛЕНТА, реализующие матрично-векторное умножение.

Программа iter реализует вычисления третьего блока и включена в пакет ИТЕРА. Вычисление вектора г1"' = rk + а(г - Ахк) представляет собой триаду, вычисление которой легко векторизуется. При этом одновременно подсчитывается норма вектора г-Ахк. Если она меньше е, то в память на место а заносится 0 (показатель получения решения), в на место к - номер итерации, на которой требуемое решение было получено. Вычисленный вектор fkt| дублируется в памяти два раза - для использования его на этапе решения системы Mxktl= fk,i и на этапе вычисления jk-tz

Время работы программы: 5n + (3T(RS) + T(WS))(m + 1) + Т(£), где п = 64т + в, Т(£) - время суммирования последних компонент частичных сумм скалярных произведений.

Перед выполнением четвертого блока вычислений необходимо проверить содержимое ячейки ОП, которая была отведена для записи числа а. В случае наличия в ней нулевого элемента вычисления следует прекратить и перенести результат выполнения последней итерации хк на место правой части. В противном случае необходимо решать систему М xkfi - fk".

Приведена модель распределения памяти, обеспечивающая слаженное выполнение всех четырех блоков вычислений.

В п.2.3.2 приводятся примеры использования программ пакета ЛЕНТА для решения системы Мхк*' = Г1"'.

В п.2.3.3 приводятся примеры использования программ пакета ИТЕРА для решения системы Mxk*' = fk*'. Приводится блок-схема использования программы ITER для реализации вложенных итераций.

В п.2.3.4 приведена блок-схема использования программы iter для реализации поперелекных итерационных методов, общая идея которых состоит в том, что два или более итерационных процесса могут чередоваться.

В п.2.3.6 анализируются возможности использования пакета ИТЕРА. Главное достоинство пакета ИТЕРА состоит в том, что наряду с программами (bmgase, Trgase, bmcg и др.), эффективно реализующими наиболее распространенные итерационные методы для решения систем линейных уравнений с трехдиагональными, ленточными и пятидиагоналышми разреженными матрицами, в пакет включена программа iter, при помощи которой и с привлечением программ пакета ЛЕНТА можно огущйотвить решение

любой системы линейных уравнений, независимо от типа матрицы, итерационными методами, основанными на расщеплении матрицы системы.

Программа ITER может быть так же использована и для реализации метода Ричардсона :

xktl = xk + аь,(г - Ахк).

Задав последовательность итерационных констант с^, к = 1..... можно

многократно использовать программу ITER без каких либо изменений.

Программы вмсо и вмрсо, реализующие метод сопряженных градиентов и метод сопряженных градиентов с предобусловливанием матрицы системы, могут также быть использованы для решения системы линейных уравнений с произвольной матрицей. Необходимо только внести незначительные изменения, а именно, оформить программу умножения матрицы данной системы на вектор, и уметь решать систему с матрицей предобусловливания. Для этого могут быть использованы программы пакета ЛЕНТА.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ.

1. Разработаны пакеты программ ЛЕНТА и ИТЕРА для решения систем линейных уравнений с матрицами специальной структуры. Программы написаны с учетом структуры матриц; время выполнения реализованных в пакете алгоритмов оптимально.

2. Впервые реализован алгоритм конвейерной прогонки; разработан новый эффективный алгоритм решения трехдиагональных систем линейных уравнений.

3. Структура пакета ЛЕНТА предоставляет возможность выбрать для решения конкретной задачи оптимальную программу. Пакет прямых методов активно используется для решения систем итерационными методами.

4. Представленные в пакетах программы могут быть- использованы для решения более широкого класса задач.

5. Модули ITER и BMCG пакета итерационных методов могут быть без изменения использованы для эффективного решения систем с произвольными матрицами.

По теме диссертации опубликованы следующие работы :

1. Е.А.Локоть. Параллельные алгоритмы обработки матриц специальной структуры. Препринт ИСК Кибернетика, М., 1990.

2. Е.А.Локоть, М.А.Фрумкин, М.В.Храпченко. Пакет программ "ЛЕНТА" для обработки ленточных матриц на векторно-конвейерной ЭВМ. Вопросы кибернетики. Базовое программное обеспечение Супер ЭВМ, М., 1990, с.5-27.

3. Е.А.Локоть. Векторизация алгоритмов решения систем линейных уравнений с ленточными матрицами методами Якоби, Гаусса-Зейделя, сопряженных градиентов. Вопросы кибернетики. Средства моделирования и разработка прикладных программ для Супер ЭВМ, М.,1992, с.49-58..

4.. Е.А.Локоть. Параллельный алгоритм решения трехдиагональных систем линейных уравнений. Вопросы кибернетики (в печати).

¿Л* -