CC BY
101
УДК: 519.612
Р.С. Суровцев, С.П. Куксенко, Т.Р. Газизов
Многократное решение системы линейных алгебраических уравнений с помощью блочного LU-разложения для вычисления емкостной матрицы системы проводников и диэлектриков при изменении ее параметров
Приведены аналитические оценки ускорения многократного решения СЛАУ за счет использования блочного L ^/-разложения. Усовершенствован алгоритм вычисления ряда емкостных матриц системы проводников и диэлектриков. Выполнены численные оценки ускорения усовершенствованным алгоритмом в зависимости от различных параметров. Получено максимальное ускорение в 15 раз.
Ключевые слова: система линейных алгебраических уравнений, блочное L^-разложение, емкостная матрица, многократные вычисления.
В настоящее время имитационное моделирование различных процессов (механических, тепловых, электрических и пр.) получило широкое распространение в различных областях практической деятельности. Это моделирование часто основано на численных методах, сводящих задачу к решению системы линейных алгебраических уравнений (СЛАУ). Поэтому от умения быстро решать СЛАУ зависит эффективность моделирования в целом. Особая необходимость в этом возникает при моделировании в диапазоне параметров, например в задачах оптимизации по большому числу параметров. Тогда получим [1]
Akxk = bk, (1)
где Ak - плотная, неособенная и квадратная матрица СЛАУ порядка N; bk - вектор свободных членов, а xk - вектор неизвестных, k = 1, 2, ... .
Необходимо отметить, что в ряде случаев изменение некоторого параметра моделируемой структуры приводит к изменению лишь части элементов матрицы СЛАУ системы (1), а остальные элементы при этом остаются неизменными. Например, в ряде задач электромагнитной совместимости необходимо многократное вычисление матрицы электростатической индукции (далее емкостной матрицы C) системы проводников и диэлектриков методом моментов [2]. Необходимость многократного решения СЛАУ возникает при учете частотной зависимости относительной диэлектрической проницаемости диэлектрика (sr) [3]. Тогда для каждой частотной точки диапазона вычисляется емкостная матрица, для чего необходимо решить СЛАУ с различными векторами свободных членов, что можно представить в матричном виде как
SÄ=Vk, (2)
где Sk - квадратная и плотная матрица порядка N, являющаяся результатом применения метода моментов к анализируемой структуре; Vk - неизменная в ходе вычислений матрица размера NxNcond, состоящая из задаваемых потенциалов на подобластях, на которые разбиты границы структуры, а £k - искомая матрица размера N*NcoNd, дающая распределение плотности заряда на этих границах, Ncond - количество проводников, не считая опорного, k = 1, 2, ..., m, где m - число решаемых СЛАУ.
Порядок матрицы СЛАУ складывается из количества подобластей на границах проводник-диэлектрик (NC) и диэлектрик-диэлектрик (ND), а элементы матрицы вычисляются из параметров этих подобластей. При изменении sr изменяются лишь элементы с индексом больше NC, на главной диагонали матрицы СЛАУ, соответствующие подынтервалам диэлектрик-диэлектрик.
Традиционно для решения СЛАУ с плотной матрицей используются прямые методы, вычислительные затраты которых пропорциональны N3, что неприемлемо при многократных вычислениях с изменяющейся матрицей, поскольку затраты становятся пропорциональными mN3. Между тем специфика частичного изменения матрицы СЛАУ - это ресурс для уменьшения времени многократных вычислений. Поэтому актуальна разработка алгоритмов решения СЛАУ, учитывающих частичное
изменение матрицы СЛАУ при многократных вычислениях. Значительно ускорить вычисления позволяет блочная версия ¿^-разложения с последующим решением СЛАУ, поскольку при многократных вычислениях нет необходимости каждый раз выполнять полное LU-разложение матрицы СЛАУ, а нужно пересчитывать только блоки, соответствующие изменившимся элементам исходной матрицы.
С помощью блочного ¿^-разложения усовершенствован алгоритм вычисления ряда емкостных матриц, учитывающий изменение элементов матрицы СЛАУ только на главной диагонали [4]. Для данного алгоритма выполнены аналитические оценки затрат, показавшие значительный ресурс ускорения вычислений до 134 раз [5]. Алгоритм апробирован на практических задачах [6, 7]. Однако при моделировании в диапазоне параметров структуры возможно изменение не только ег, но и её размеров. Например, при оптимизации связанной линии передачи варьируется толщина диэлектрика (кс) между проводниками (рис. 1), где диапазон изменения кс снизу ограничен величиной к, а сверху - суммой к и t.
-Ж-
ГС
Рис. 1. Вид поперечного сечения связанной линии
При изменении параметра кС изменяются элементы в нижней и правой частях матрицы СЛАУ с индексами, соответствующими подынтервалам диэлектрик-диэлектрик (индексы строк и столбцов больше N0. Оценки применимости и эффективности блочного ¿^-разложения при таком изменении матрицы СЛАУ авторам неизвестны. Между тем они бы позволили расширить сферу применения алгоритмов на основе блочного ¿^-разложения к другим задачам с частичным изменением матрицы СЛАУ при многократных вычислениях, например при адаптивном итерационном выборе оптимальной сегментации структуры проводников и диэлектриков [8].
Цель работы - выполнить оценку эффективности применения блочного ¿^-разложения для решения СЛАУ при многократном вычислении емкостной матрицы в диапазоне параметров.
Для достижения поставленной цели необходимо решить следующие задачи: рассмотреть алгоритм блочного ¿^-разложения и выполнить возможные аналитические оценки эффективности его применения; выполнить совершенствование исходного алгоритма; выполнить оценку ускорения усовершенствованным алгоритмом на основе вычислительного эксперимента.
Алгоритм многократного решения СЛАУ с помощью блочного ¿^-разложения и аналитические оценки его эффективности. Рассмотрим алгоритм блочного ¿^-разложения на примере разложения матрицы 8, представленной в виде
S =
A C
B D
где блоки размера: А - N^N4, В - N^N0, С - N^N4, Б - N^N0.
Тогда ¿^-разложение матрицы 8 будет иметь следующий блочный вид (в случае, если блок А -неособенная матрица) [9]:
S =
A C
B D
^ L =
I
CA-1
U =
A
0
B
D - CA -1B
(3)
где I - единичная матрица. Из (3) видно, что необходимо несколько раз выполнять обращение блока А, что увеличивает затраты времени при программной реализации ¿^-разложения. Также увеличиваются затраты памяти на хранение отдельно матриц Ь и и. Поэтому для минимизации затрат, перепишем (3) в виде (матрица 8'), удобном для хранения и дальнейшего использования:
S' =
-S'11 S'12 " A-1 A-1B
_S'21 S'22 _ С D - CA -1B
(4)
Тогда алгоритм блочного ¿^-разложения в удобном для программной реализации виде можно представить последовательностью действий: 8'ц = А-1, 8'21 = С, 8\2 = 8'цВ, 8'22 = Б - 8'218'12. Видно, что при необходимости многократного решения СЛАУ изменение элементов блока А приводит к необходимости полного пересчета верхнетреугольной части матрицы, поэтому применение блочного ¿^-разложения будет неэффективно. Однако при изменении элементов любого из других блоков (В, С или Б) необходимо пересчитать лишь блоки матрицы 8', зависящие от изменившихся блоков
w
s
t
s
r
h
С
исходной матрицы. Таким образом, в самом общем случае (изменение блоков В, С, Б) при да-кратных вычислениях будет однократно выполняться вычислительно затратное первое разложение матрицы 8' (включающее обращение блока А) и последующие т - 1 вычислений только блоков ^12, 8^ и 8*22.
Для аналитической оценки ускорения многократного решения СЛАУ рассмотрим (аналогично работам [5, 10]) отношение (в) общего времени решения т СЛАУ последовательным алгоритмом ЬИ-разложения ко времени решения блочным алгоритмом:
Р Т + (т—Щ ' (5)
где Т1 - время первого решения, в которое входит нахождение блока 8'11 размером и
последующее решение СЛАУ с нахождением матрицы неизвестных из (2); Ts - время вычисления блоков 8'12, 8'22 с последующим нахождением матрицы неизвестных. Из (5) получим оценку максимально возможного ускорения:
Ртах =НтТ^Ц— = ^ . (6)
т^ю Т1 +(т — 1)Т$ Ts
Из выражений (5)-(6) следует, что чем больше да, тем меньше ускорение зависит от времени первого решения. Также видно, что величина ускорения обратно пропорциональна времени вычисления блоков 8'12, 822 и определяется размерами этих блоков. Таким образом, при больших т и N и малом N0 можно получить значительное ускорение многократного решения СЛАУ.
Совершенствование алгоритма вычисления ряда емкостных матриц с помощью блочного ЬЦ-разложения. Для ясности сначала рассмотрим исходный алгоритм (алгоритм 1).
Алгоритм 1. Исходный алгоритм вычисления ряда емкостных матриц.
1. Вычислить элементы матрицы воздействия Ук (размера
2. Для к от 1 до т.
3. Вычислить элементы матрицы 8к (размера ^Н).
4. Выполнить ЬИ-разложение матрицы 8к.
5. Найти матрицу решения £к из уравнения 8к£к = Ук.
6. Вычислить элементы матрицы Ск, основываясь на элементах Ък.
7. Увеличить к.
Как видно из алгоритма 1, на каждом шаге к выполняются заполнение матрицы СЛАУ, решение СЛАУ с помощью ¿^-разложения и вычисление емкостной матрицы Ск.
Как было отмечено выше, порядок матрицы N складывается из количества подобластей на границах проводник-диэлектрик - N0 и диэлектрик-диэлектрик - N0. Но, как в задаче оптимизации связанной линии, может изменяться высота не всей границы раздела воздушной среды и диэлектрической подложки, а только ее часть, находящаяся между проводниками (рис. 2). Поэтому, в такой задаче N1) складывается из количества подобластей с постоянными (^О0™5') и изменяющимися (N0™) параметрами. Таким образом, общий порядок неизменяющихся элементов составляет =N0 + так что N = ^„й + N0™. Далее для уменьшения количества вводимых индексов и N0™ приняты равными N и N0 соответственно. Для общего случая перед вычислениями необходимо изменить нумерацию границ структуры, чтобы подобласти, соответствующие изменяющимся строкам и столбцам матрицы СЛАУ, нумеровались последними. Тогда алгоритм вычисления ряда емкостных матриц с учетом (4) можно представить в следующем виде.
Алгоритм 2. Усовершенствованный алгоритм вычисления ряда емкостных матриц.
1. Вычислить элементы матрицы 81 (размера ^Ж).
2. Изменить нумерацию границ структуры так, чтобы подобласти, соответствующие изменяющимся строкам и столбцам матрицы 81, нумеровались последними.
3. Вычислить А1 = Аг1 (размеры блока NA^NA).
4. Вычислить элементы матрицы воздействия Ук.
5. Х0 = А1У0 (размер блоков Х0 и У0 - ^^тж).
6. Для к от 1 до да.
7. Вк = А:Вк.
8. Бк = Бк - СкВк.
9. X! к = У1 - СкХ0 (размер блоков Хи и У! - ^^тж).
10. = Dt-1Xlt.
11. Sot = Xo - BkZik.
12. Вычислить элементы емкостной матрицы C k.
13. Вычислить элементы изменяющихся блоков Bk+1, Ck+1, Dk+1.
14. Увеличить k.
Из алгоритма 2 видно, что для реализации блочного LU-разложения необходима дополнительная матрица X порядка NxNcond, в которой хранятся результаты промежуточных вычислений.
Вычислительный эксперимент. В пакете Microsoft Visual Studio выполнена программная реализация алгоритмов 1 и 2 для вычислительной оценки ускорения. Использована рабочая станция (без параллельных вычислений, т.е. работало одно ядро процессора) со следующими параметрами: платформа AMD FX(tm)-8350 Eight-Core Processor; частота процессора 4,01 ГГц; объем ОЗУ 32 Гб; число ядер - 8; операционная система Windows 7x64. Количество СЛАУ m выбрано, исходя из того, что в задаче оптимизации связанной линии число изменений матрицы СЛАУ ограничено только значением t, которое однозначно определяется типовыми значениями толщины проводящего слоя печатных плат и может варьироваться от 5 до 105 мкм. Поэтому для простоты число изменений hC принято равным m = 10, 20, ..., 100. Так как в связанной линии (см. рис. 1) два проводника (не считая опорного), то Ncond = 2.
Сначала выполнена оценка зависимости отношения (ß) времени вычислений исходным алгоритмом 1 (ТИ) ко времени вычислений усовершенствованным алгоритмом 2 (Ту) от порядка блока матрицы СЛАУ с неизменяющимися элементами (NA). Данная оценка позволит определить величину минимального порога NA/N, после которого усовершенствованный алгоритм будет работать быстрее, чем исходный. Вычисление выполнено для крайних точек диапазона изменения m. Полученные ускорения (ß = ТИ/ТУ) для матрицы размером N = 1000 сведены в табл. 1.
Таблица 1
Оценка ускорения решения СЛАУ за счет использования усовершенствованного алгоритма __ вычисления т емкостных матриц при У=1000___
Na/N 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9
m = 10 0,17 0,23 0,31 0,41 0,55 0,68 0,83 0,98 1,17
m = 100 0,17 0,23 0,31 0,42 0,57 0,77 1,05 1,48 2,57
Из табл. 1 видно, что ускорение возрастает при росте но для разных т характер роста отличается. При NA/N < 0,4 ускорение не зависит от т, однако при дальнейшем увеличении ^/У рост ускорения для т = 100 выражен сильнее чем для т = 10. Так, значение ^/у при котором временные затраты на вычисление исходным и усовершенствованным алгоритмами равны, для т = 10 близко к 0,8, а для т = 100 близко к 0,7. Максимальное ускорение для т = 10 составляет 1.17 раза, а для т = 100 - 2,57 раза. Подводя итог, можно сделать следующие выводы. Величина порога ^/у после которого усовершенствованный алгоритм работает быстрее, чем исходный, зависит от т и при малых т выше, чем при больших. При малых N и т использование усовершенствованного алгоритма неэффективно, что согласуется с аналитическими оценками.
Следующей выполнена оценка зависимости времени вычисления исходным и усовершенствованным алгоритмами от т для N4^ = 0,9; 0,99; 0,999. В качестве примера в табл. 2 сведено нормированное по минимальному значению время вычисления ряда емкостных матриц при N = 1000.
Таблица 2
Нормированное по минимальному значению время вычисления т емкостных матриц при У=1000
Na/N m
10 20 30 40 50 60 70 80 90 100
Исходный алгоритм
- 1,00 2,02 3,01 4,03 5,02 6,01 7,03 8,03 9,03 10,03
Усовершенствованный алгоритм
0,9 1,0 1,36 1,75 2,13 2,52 2,90 3,27 3,69 4,03 4,46
0,99 1,0 1,04 1,09 1,14 1,19 1,24 1,28 1,33 1,38 1,43
0,999 1,0 1,0 1,0 1,01 1,02 1,02 1,03 1,03 1,05 1,04
Из табл. 2 видно, что время вычислений исходным алгоритмом имеет линейную зависимость от т, которая обусловлена тем, что исходный алгоритм основан на последовательном ЬИ-разложении. Для усовершенствованного алгоритма зависимость от т нелинейна и при увеличении ^^ становится слабее. Например, время вычислений при увеличении т от 10 до 100 для NA/N = 0,9 возрастает в 4,46 раза, а для = 0,99 - в 1,43 раза. Для = 0,999 зависимости времени от т практически не наблюдается и для всех да нормированное время близко к 1. Такую слабую зависимость времени вычислений от т при NA/N = 0,999 можно объяснить высокими вычислительными затратами на однократное обращение блока А размером NA = 999 и малыми затратами на т вычислений оставшихся блоков (состоящих из одного столбца и одной строки) и решений СЛАУ. Таким образом, за счет малых затрат на две последние операции увеличение да не приводит к значительному росту времени вычислений усовершенствованным алгоритмом.
Для оценок пределов ускорения усовершенствованным алгоритмом и зависимости времени вычислений от N в табл. 3 сведено ускорение (Р) вычислений т емкостных матриц для N = 1000, 2000 и 3000 при NA/N = 0,9; 0,99; 0,999.
Таблица 3
Оценка ускорения решения СЛАУ за счет использования усовершенствованного алгоритма
вычисления т емкостных матриц при N = 1000, 2000, 3000
N МА т
10 20 30 40 50 60 70 80 90 100
1000 900 1,17 1,68 1,96 2,14 2,26 2,36 2,43 2,50 2,54 2,57
990 1,37 2,58 3,65 4,70 5,60 6,50 7,28 8,03 8,75 9,35
999 1,39 2,74 3,98 5,36 6,65 7,99 9,26 10,56 11,81 13,00
2000 1800 1,31 1,89 2,22 2,44 2,59 2,69 2,77 2,83 2,90 2,94
1980 1,51 2,89 4,16 5,32 6,37 7,34 8,25 9,06 9,87 10,57
1998 1,54 3,07 4,59 6,10 7,58 9,06 10,35 11,92 13,38 14,75
3000 2700 1,35 1,95 2,29 2,52 2,67 2,78 2,88 2,95 2,99 3,04
2970 1,55 2,97 4,25 5,44 6,53 7,54 8,48 9,34 10,12 10,87
2997 1,59 3,14 4,70 6,30 7,76 9,26 10,78 12,26 13,69 15,17
Как видно из табл. 3, подтверждаются оценки, сделанные выше. Ускорение вычислений возрастает по мере роста т и NA. Пределы роста Р разные и зависят от размера блоков с изменяющимися элементами. Например, при увеличении NA от 900 до 990 ускорение возрастает в 3,64 раза. Примечателен пример с NA = 999, так как в данном случае изменяются только одна строка и один столбец исходной матрицы. Максимальное ускорение при таком изменении матрицы СЛАУ составляет 13 раз. Также из табл. 3 видно, что при увеличении порядка исходной матрицы ускорение возрастает. При увеличении N от 1000 до 3000 максимальное ускорение возрастает от 13 до 15,17 раза. Таким образом, полученные оценки зависимости (вЫ) говорят о слабом росте ускорения при увеличении N.
Заключение. Выполнено усовершенствование алгоритма вычисления ряда емкостных матриц с помощью блочного Ш-разложения для случая, когда в структуре проводников и диэлектриков изменяются любые параметры. Проведены аналитические и вычислительные оценки эффективности усовершенствованного алгоритма для многократных вычислений. Так, показано, что при малых NA и т использование усовершенствованного алгоритма не эффективно, что согласуется с аналитическими оценками. Также выявлено, что при NA/N = 0,999 для всех т нормированное по минимальному значению время вычислений близко к 1. Оценка влияния размера матрицы СЛАУ на ускорение многократных вычислений показала его рост, хотя и незначительный. Примечательно, что применение блочного ¿^-разложения позволяет получить ускорение, начиная с довольно малых значений т > 10 и не с самых высоких значений NA/N > 0,9, что перекрывает весьма широкий круг задач на практике.
Дополнительно необходимо отметить, что в дальнейшем целесообразны более детальные сравнительные оценки усовершенствованного и исходного алгоритмов, которые позволят предварительно, до проведения вычислительного эксперимента, определить необходимость применения одного из алгоритмов. Примечательно, что для таких оценок могут быть получены аналитические выражения эффективности усовершенствованного алгоритма. Между тем детальные вычислительные и
аналитические оценки ускорения усовершенствованным алгоритмом позволят расширить сферу применения алгоритмов на основе блочного ¿^-разложения к другим задачам с частичным изменением матрицы СЛАУ при многократных вычислениях.
Вычисление емкости связанной линии методом моментов выполнено за счет гранта Российского научного фонда (проект № 14-19-01232) в ТУСУРе. Алгоритм для решения СЛАУ разработан при поддержке грантов РФФИ 14-07-31267 и 14-29-09254. Оценки ускорения выполнены в рамках государственного задания № 8.1802.2014/K Минобрнауки России.
Литература
1. Calrago C. Incremental incomplete LU factorization with application to time-dependent PDEs / C. Calrago, J.P. Chehab, Y. Saad // Numer. Lin. Algebra with Appl. - 2010. - Vol. 17(5). - P. 811-837.
2. Gazizov T.R. Analytic expressions for Mom calculation of capacitance matrix of two dimensional system of conductors and dielectrics having arbitrary oriented boundaries // Proc. of the 2001 IEEE EMC Symposium. - Montreal, Canada, 2001. - Vol. 1. - P. 151-155.
3. Djordjevich A.R. Wideband frequency-domain characterization of FR-4 and time-domain causality / A.R. Djordjevich, R.M. Biljic, V.D. Likar-Smiljanic, T.K. Sarkar // IEEE Trans. Electromag. Compat. -2001. - Vol. 43. - P. 662-666.
4. Куксенок С.П. Совершенствование алгоритма вычисления методом моментов ёмкостных матриц системы проводников и диэлектриков в диапазоне значений диэлектрической проницаемости диэлектриков / С.П. Куксенко, Т.Р. Газизов // Электромагнитные волны и электронные системы. - 2012. - №10. - C. 13-21.
5. Суровцев Р.С. Аналитическая оценка вычислительных затрат на решение СЛАУ при многократном вычислении емкостной матрицы в диапазоне изменения диэлектрической проницаемости диэлектриков / Р.С. Суровцев, С.П. Куксенко, Т.Р. Газизов // Записки научных семинаров ПОМИ. -2014. - №428. - С. 196-207.
6. Суровцев Р.С. Вычисление матрицы емкостей произвольной системы проводников и диэлектриков методом моментов: зависимость ускорения за счет блочного LU-разложения от порядка матрицы СЛАУ / Р.С. Суровцев, С.П. Куксенко // Изв. вузов. Физика. - 2012. - Т. 55, № 9/3. - С. 126-130.
7. Суровцев Р.С. Использование блочного ¿^-разложения для ускорения вычисления временного отклика связанных линий передачи с учётом частотной зависимости диэлектрической проницаемости подложки / Р.С. Суровцев, В.К. Салов, С.П. Куксенко // Инфокоммуникационные технологии. - 2013. - Т. 11, №3. - С. 64-69.
8. Аширбакиев Р.И. Адаптивный итерационный выбор оптимальной сегментации границ проводников и диэлектриков в задачах электростатики / Р.И. Аширбакиев, В.К. Салов // Доклады Том -ского государственного университета систем управления и радиоэлектроники. - 2013. - № 3(29), ч. 1. - С. 159-161.
9. Highman N.J. Accuracy and Stability of Numerical Algorithms. - Philadelphia: SIAM, 1961. -680 p.
10. Ахунов Р.Р. Многократное решение СЛАУ с частично изменяющейся матрицей итерационным методом / Р.Р. Ахунов, С.П. Куксенко, В.К. Салов, Т.Р. Газизов // Записки научных семинаров ПОМИ. - 2013. - № 419. - С. 16-25.
Суровцев Роман Сергеевич
Аспирант каф. телевидения и управления (ТУ) ТУСУРа
Тел.: 8 (382-2) 41-34-39
Эл. почта: [email protected]
Куксенко Сергей Петрович
Канд. техн. наук, доцент каф. ТУ Тел.: 8 (382-2) 41-34-39 Эл. почта: [email protected]
Газизов Тальгат Рашитович
Д-р техн. наук, ст. науч. сотрудник, зав. каф. ТУ Тел.: 8 (382-2) 41-34-39 Эл. почта: [email protected]
Surovtsev R.S., Kuksenko S.P., Gazizov T.R.
Multiple solution of linear algebraic systems by means of block LU-decomposition for computation of capacitance matrix of system consisting of conductors and dielectrics, and having its parameters varied
Analytical evaluations for speed up of the multiple solution of linear algebraic systems via block LUdecomposition have been obtained. A computational algorithm for some sequence of capacitance matrixes of the system consisting of conductors and dielectrics has been improved. Numerical evaluations of the speed up obtained via the improved algorithm, depending on different parameters, have been carried out. The maximum acceleration equals 15 times.
Keywords: linear algebraic system, block LU-decomposition, capacitance matrix, multiple solution.
