MULTIPLE COMPUTATIONS OF CAPACITANCE MATRIX OF CONDUCTOR AND DIELECTRIC SYSTEM WITH VARYING PARAMETERS BY BLOCK LU-DECOMPOSITION

Full Text

Введение В настоящее время имитационное моделирование различных процессов (механических, тепловых, электрических и пр.) получило широкое распространение в различных областях практической деятельности. Это моделирование часто основано на численных методах, сводящих задачу к решению системы линейных алгебраических уравнений (СЛАУ). Поэтому от умения быстро решать СЛАУ зависит эффективность моделирования в целом. Особая необходимость в этом возникает при моделировании в диапазоне параметров, например, в задачах оптимизации по большому числу параметров. Тогда получим [1] , (1) где Ak - плотная квадратная матрица СЛАУ порядка N, bk - вектор свободных членов, а xk - вектор неизвестных, k = 1, 2, … Необходимо отметить, что в ряде случаев изменение некоторого параметра моделируемой структуры приводит к изменению лишь части элементов матрицы СЛАУ системы (1), а остальные элементы при этом остаются неизменными. Например, в ряде задач электромагнитной совместимости необходимо многократное вычисление матрицы электростатической индукции (далее емкостной матрицы C) системы проводников и диэлектриков методом моментов [2]. Необходимость многократного решения СЛАУ возникает при учете частотной зависимости относительной диэлектрической проницаемости диэлектрика [3]. Тогда, для каждой частотной точки диапазона вычисляется емкостная матрица, для чего необходимо решить СЛАУ с различными векторами свободных членов, что можно представить в матричном виде как , (2) где Sk - квадратная и плотная матрица порядка N, являющаяся результатом применения метода моментов к анализируемой структуре, Vk - неизменная в входе вычислений матрица размера N × NCOND, состоящая из задаваемых потенциалов на подобластях, на которые разбиты границы структуры, а - искомая матрица размера N × NCOND, дающая распределение плотности заряда на этих границах, NCOND - количество проводников, не считая опорного, k = 1, 2, … m. Порядок матрицы СЛАУ складывается из количества подобластей на границах проводник-диэлектрик NC и диэлектрик-диэлектрик ND, а элементы матрицы вычисляются из параметров этих подобластей. При изменении изменяются лишь элементы с индексом больше NC, на главной диагонали матрицы СЛАУ, соответствующие подынтервалам диэлектрик-диэлектрик. Традиционно для решения СЛАУ с плотной матрицей используются прямые методы, вычислительные затраты которых пропорциональны N 3, что неприемлемо при многократных вычислениях с изменяющейся матрицей, поскольку затраты становятся пропорциональными mN 3. Между тем, специфика частичного изменения матрицы СЛАУ - это ресурс для уменьшения времени многократных вычислений. Поэтому актуальна разработка алгоритмов решения СЛАУ, использующих частичное изменение матрицы СЛАУ при многократных вычислениях. Значительно ускорить вычисления позволяет блочная версия LU-разложения с последующим решением СЛАУ, поскольку при многократных вычислениях нет необходимости каждый раз выполнять полное LU-разложение матрицы СЛАУ, а нужно пересчитывать только блоки, соответствующие изменившимся элементам исходной матрицы. С помощью блочного LU-разложения усовершенствован алгоритм вычисления ряда емкостных матриц, учитывающий изменение элементов матрицы СЛАУ только на главной диагонали [4]. Для данного алгоритма выполнены аналитические оценки затрат, показавшие значительный ресурс ускорения вычислений до 134 раз [5]. Алгоритм апробирован на практических задачах [6-7]. Однако, при моделировании в диапазоне параметров структуры возможно изменение не только , но и ее размеров. Например, при оптимизации связанной линии передачи варьируется толщина диэлектрика hC между проводниками (см. рис. 1), где диапазон изменения hC снизу ограничен величиной h, а сверху - суммой h и t. Рис. 1 Вид поперечного сечения связанной линии При изменении параметра hС изменяются элементы в нижнем и правом блоках матрицы СЛАУ с индексами, соответствующими подынтервалам диэлектрик-диэлектрик (индексы строк и столбцов больше NC). Показана возможность применения блочного LU-разложения при таком изменении матрицы СЛАУ и предложена «вариативная перенумерация» подобластей моделируемой структуры для более эффективного решения [8]. Однако более детальные оценки применимости и эффективности блочного LU-разложения при таком изменении матрицы СЛАУ ещё не сделаны. Между тем, они бы позволили расширить сферу применения алгоритмов на основе блочного LU-разложения к другим задачам с частичным изменением матрицы СЛАУ при многократных вычислениях, например, при адаптивном итерационном выборе оптимальной сегментации структуры проводников и диэлектриков [9]. Цель работы - выполнить оценку эффективности применения блочного LU-разложения для решения СЛАУ при многократном вычислении емкостной матрицы структуры проводников и диэлектриков в диапазоне параметров. Для достижения поставленной цели представляется необходимым: выявить возможную для использования специфику матриц, проанализировать роль перенумерации, рассмотреть алгоритм и реализацию блочного LU-разложения и выполнить возможные аналитические оценки эффективности его применения; усовершенствовать исходный алгоритм; оценить ускорение усовершенствованным алгоритмом на основе вычислительного эксперимента. Использование специфики изменения матриц и перенумерации Рассмотрим специфику матриц, а точнее, специфику их изменения от матрицы к матрице. Отметим, что для применения блочного LU-разложения важна не степень изменений матрицы (весьма существенная, например, для решения итерационным методом), а расположение изменяющихся элементов. Кроме того, для повышения эффективности этого применения важна сама возможность изменения их расположения ближе к правой и нижней границам матрицы. Рассмотрим эти вопросы по отдельности. Для определенности и ясности целесообразно сделать это на рассматриваемой задаче вычисления ёмкостной матрицы методом моментов, а после обобщить на произвольную задачу. Как уже отмечалось, в исходном алгоритме вычисления первая часть индексов элементов матрицы СЛАУ соответствует границам проводников структуры, а вторая - диэлектриков. Таким образом, могут быть заданы такие изменения её параметров, которые затронут только диэлектрики или часть их, не затрагивая проводников. Кроме того, протяженность границ диэлектриков может быть много меньше, чем проводников. Таким образом, может существовать широкий ряд структур, естественным образом (то есть без участия пользователя) дающих изменения значений небольшого числа элементов матрицы СЛАУ, расположенных только в правом и нижнем блоках матрицы. В используемых авторами алгоритмах вычисления емкостной матрицы есть ещё одна специфика: в частных случаях ориентации границ проводников и диэлектриков ортогонально одной из декартовых осей, значения соответствующих этим ориентациям элементов матрицы СЛАУ вычисляются по своим аналитическим выражениям, и поэтому эти элементы сгруппированы в своих строках и столбцах. (Детальная иллюстрация этого приведена в работе [8].) Это приводит к еще большей структурированности матрицы СЛАУ и возможности изменений её элементов, расположенных (опять же, без участия пользователя) только в правых строках и нижних столбцах. В результате, существует большое число реальных задач, когда матрицы последовательно решаемых СЛАУ полностью готовы для решения с помощью блочного LU-разложения, причем это решение может быть эффективным в силу малого (относительно порядка матрицы) числа изменяющихся строк и столбцов. Между тем, в самом общем случае, изменяющиеся элементы матрицы СЛАУ могут быть расположены произвольно. Поэтому, для повышения эффективности применения блочного LU-разложения важна возможность изменения их расположения: как можно ближе к правой и нижней границам матрицы. Примечательно, что такая возможность вполне может существовать, и даже быть легко реализуемой. Так, в алгоритме рассматриваемой задачи значения элементов матрицы СЛАУ выражают собой связь между подынтервалами, расположенными на границах проводников и диэлектриков структуры, и во многом определяются значениями декартовых координат этих подынтервалов. При этом сама нумерация подынтервалов не влияет на значения элементов матрицы. Таким образом, если подынтервалы пронумерованы и элементы матрицы вычислены, то при изменении нумерации подынтервалов не требуется повторное вычисление элементов матрицы. Более того, не требуется даже перестановка элементов, а достаточно лишь изменить ссылки (указатели) на них, что требует минимальных вычислительных затрат. Отметим, что номера подынтервалов, изменивших свои координаты, можно легко найти по номерам границ (число которых гораздо меньше, например, на рис. 1 их 11: на проводниках - 8, на диэлектрике - 3), изменивших свои координаты - поэтому и здесь вычислительные затраты очень малы. В итоге оказывается, что если выполнена первая нумерация подынтервалов и вычислены элементы матрицы СЛАУ для любой заданной структуры и заданы изменения ее параметра, то можно и целесообразно выполнить перенумерацию её подынтервалов, перемещающую все изменяющиеся элементы матрицы в правый и нижний блоки матрицы. Важно отметить, что введение этой процедуры автоматически оптимизирует любую структуру для многократного решения СЛАУ с помощью блочного LU-разложения. Поэтому она представляется авторам важным звеном, позволяющим получить ускорение решения на практике: нет матриц с малыми размерами правого и нижнего блоков - нет значительного ускорения. Алгоритм многократного решения СЛАУ с помощью блочного LU-разложения и его реализация Рассмотрим алгоритм блочного LU-разложения на примере разложения матрицы S, представленной в виде , где блоки размера: A - NA×NA, B - NA×ND, C - ND×NA, D - ND×ND. Тогда LU-разложение матрицы S будет иметь следующий блочный вид (в случае, если блок A - неособенная матрица) [10]: ,, (3) где I - единичная матрица. Из (3) видно, что необходимо несколько раз выполнять обращение блока A, что увеличивает затраты времени при программной реализации LU-разложения. Также увеличиваются затраты памяти на хранение отдельно матриц L и U. Поэтому для минимизации затрат, перепишем (3) в виде (матрица S') удобном для хранения и дальнейшего использования: . (4) Тогда алгоритм блочного LU-разложения в удобном для программной реализации виде можно представить последовательностью действий: S'11 = A-1, S'21 = C, S'12 = S'11B, S'22 = D-S'21S'12. Видно, что при необходимости многократного решения СЛАУ изменение элементов блока A приводит к необходимости полного пересчета верхнетреугольной части матрицы, поэтому применение блочного LU-разложения будет неэффективно. Однако, при изменении элементов любого из других блоков (B, C или D) необходимо пересчитать лишь блоки матрицы S', зависящие от изменившихся блоков исходной матрицы. Таким образом, в самом общем случае (изменение блоков B, C, D), при m-кратных вычислениях будет однократно выполняться вычислительно затратное первое разложение матрицы S' (включающее обращение блока A) и последующие m-1 вычислений только блоков S'12 и S'22. Совершенствование алгоритма вычисления ряда емкостных матриц с помощью блочного LU-разложения При вычислении емкостной матрицы используется последовательное LU-разложение, поскольку требуется решать СЛАУ с несколькими векторами свободных членов, а данный подход широко используется для решения подобных задач и хорошо освещен в литературе. (Помимо этого подхода, подобные задачи могут быть решены с помощью блочных версий итерационных методов, однако в данной работе они не рассматриваются, в частности, поскольку для итерационных методов трудно выполнить точные аналитические оценки арифметических затрат до начала решения из-за неопределенности требуемого числа итераций.) Для ясности сначала рассмотрим исходный алгоритм 1. Алгоритм 1 Исходный алгоритм вычисления ряда емкостных матриц 1 Вычислить элементы матрицы воздействия Vk (размера N×NCOND) 2 Для k от 1 до m 3 Вычислить элементы матрицы Sk (размера N×N) 4 Выполнить LU-разложение матрицы Sk 5 Найти матрицу решения Sk из уравнения SkSk=Vk 6 Вычислить элементы матрицы Ck, основываясь на элементах Sk 7 Увеличить k Как видно из алгоритма 1, на каждом шаге k выполняются заполнение матрицы СЛАУ, решение СЛАУ с помощью LU-разложения и вычисление емкостной матрицы Сk. Как было отмечено выше, порядок матрицы N складывается из количества подобластей на границах проводник-диэлектрик - NC и диэлектрик-диэлектрик - ND. Но, как в задаче оптимизации связанной линии, может изменяться высота не всей границы раздела воздушной среды и диэлектрической подложки, а только ее часть, находящаяся между проводниками (см. рис. 1). Поэтому в такой задаче ND складывается из количества подобластей с постоянными (NDconst) и изменяющимися (NDvar) параметрами. Алгоритм 2 Усовершенствованный алгоритм вычисления ряда емкостных матриц 1 Вычислить элементы матрицы S1 (размера N×N) 2 Перенумерация подобластей структуры 3 Вычислить A1 = A1-1 (размеры блока NA×NA) 4 Вычислить элементы матрицы воздействия Vk 5 X0 = A1V0 (размер блоков X0 и V0 - NA×NCOND) 6 Для k от 1 до m 7 Bk = A1Bk 8 Dk = Dk - CkBk 9 X1 k = V1 - CkX0 (размер блоков X1k и V1 - ND×NCOND) 10 S1k = Dk-1X1k 11 S0k = X0 - BkS1k 12 Вычислить элементы емкостной матрицы Ck 13 Вычислить элементы изменяющихся блоков Bk+1, Ck+1, Dk+1 14 Увеличить k Таким образом, общий порядок неизменяющихся элементов составляет Nconst = NC + NDconst, так что N = Nconst + NDvar. Далее для уменьшения количества вводимых индексов Nconst и NDvar приняты равными NA и ND соответственно. Для общего случая перед повторными вычислениями необходимо изменить нумерацию границ структуры, чтобы подобласти, соответствующие изменяющимся строкам и столбцам матрицы СЛАУ, нумеровались последними. Способ вычисления ряда емкостных матриц с учетом (4) можно представить в виде алгоритма 2 - из которого видно, что для реализации блочного LU-разложения необходима дополнительная матрица X порядка N × NCOND, в которой хранятся результаты промежуточных вычислений. Оценки арифметических затрат Перед оценкой ускорения усовершенствованным алгоритмом необходимо отметить, что в алгоритмах 1 и 2 содержатся вспомогательные операции, которые необходимы для формирования элементов матрицы СЛАУ из исходной структуры проводников и диэлектриков, ее модификации при изменении параметров структуры, а также для вычисления элементов емкостной матрицы из полученных векторов решения. Эти затраты много меньше затрат на решение СЛАУ, поэтому далее они не учитывались. Из алгоритма 1 видно, что его основная часть состоит из LU-разложения матрицы Sk и решения уравнения . Оценка арифметической сложности (QLU) алгоритма 1 нахождения ряда емкостных матриц выполнена в [5]: Как видно из алгоритма 2, затраты на его реализацию состоят из однократного обращения блока A и m - 1 обращений блока D, а также нескольких матричных вычитаний и умножений. Оценки арифметической сложности на реализацию этих операций также выполнены в [5]. Полученные на их основе оценки арифметической сложности операций алгоритма 2 приведены в таблице 1. Таблица 1. Оценки арифметической сложности решения СЛАУ в алгоритме 2 № шага алгоритма Арифметическое действие Арифметическая сложность (Q) Однократно выполняемые операции 3 А1 = А1-1 5 X0 = A1V0 m-кратно выполняемые операции 7 Bk = A1Bk 8 Dk = Dk - CkBk 9 X1k = V1 - CkX0 10 S1k = Dk -1X1k 11 S0k =X0 - Bk S1k Далее, используя данные таблицы 2, получим арифметическую сложность первого решения СЛАУ и m - 1 последующих вычислений . Тогда для решения СЛАУ при реализации алгоритма 2 (Qf + mQm-1) требуется операций. Итоговым аналитическим выражением для ускорения решения СЛАУ является отношение () количества операций алгоритма 1 к количеству операций алгоритма 2: .(5) Аналитические оценки ускорения вычислений усовершенствованным алгоритмом () m емкостных матриц для N = 1000, 2000 и 3000 при NA/N = 0,9; 0,99; 0,999, выполненные по (5), сведены в таблицу 2. Число СЛАУ m выбрано, исходя из того, что в задаче оптимизации связанной линии число изменений матрицы СЛАУ ограничено только значением t, которое однозначно определяется типовыми значениями толщины проводящего слоя печатных плат и может варьироваться от 5 до 105 мкм. Поэтому, для простоты, число изменений hC принято равным m = 10; 20 … 100. Так как в связанной линии (см. рис. 1) имеются два проводника, не считая опорного, то NCOND = 2. Таблица 2. Аналитические оценки ускорения решения СЛАУ N NA m 10 20 30 40 50 60 70 80 90 100 1000 900 1,77 2,39 2,70 2,89 3,02 3,11 3,17 3,23 3,27 3,30 990 2,40 4,49 6,31 7,92 9,35 10,63 11,78 12,83 13,78 14,64 999 2,50 4,96 7,39 9,78 12,13 14,46 16,75 19,00 21,23 23,42 2000 1800 1,77 2,39 2,70 2,89 3,01 3,10 3,17 3,22 3,27 3,30 1980 2,40 4,48 6,30 7,91 9,33 10,61 11,76 12,81 13,75 14,62 1998 2,49 4,95 7,37 9,76 12,11 14,43 16,71 18,96 21,18 23,37 3000 2700 1,77 2,38 2,70 2,88 3,01 3,10 3,17 3,22 3,27 3,30 2970 2,40 4,47 6,29 7,90 9,33 10,61 11,76 12,80 13,75 14,61 2997 2,49 4,95 7,37 9,75 12,10 14,42 16,70 18,95 21,17 23,35 На основе данных таблицы 2 можно сделать несколько выводов. Ускорение практически не зависит от N. Наименьшая эффективность усовершенствованного алгоритма наблюдается при малых NA, например, при NA = 900 ускорение изменяется от 1,77 до 3,3 раза. Однако эффективность значительно возрастает при стремлении NA к N, например, при увеличении NA от 900 до 999 для m = 100 ускорение достигает 23 раз. Пример с NA = 999 примечателен, так как в данном случае изменяются только одна строка и один столбец исходной матрицы. Также видно, что ускорение сильно зависит от m. Так, по мере роста m от 10 до 100 ускорение для NA = 900 возрастает в 1,86 раза, а для NA = 999 в 9,37 раза. Таким образом, ключевыми параметрами, влияющими на время вычислений, являются порядок блока A и число решений СЛАУ. Отметим, что выполнение аналитических оценок ускорения позволяет быструю и априорную оценку эффективности и целесообразности использования того или иного метода (здесь, блочного LU-разложения) для многократного решения СЛАУ. В этой связи, принципиально важно оценить, при достижении каких значений параметров СЛАУ появляется ускорение решения. Как видно из таблицы 2, для параметра m это значение меньше 10. Поэтому выполнены аналогичные оценки для m < 10. Результаты оказались также почти независящими от N, но примечательны зависимостью от других параметров. Для наглядности они приведены графически (рис. 2). Видно, что ускорение начинается с m = 4 для NA = 999; 990; 900 и сохраняется даже при уменьшении NA до 700, хотя уже не увеличивается с ростом m. Рис. 2. Зависимости ускорения (аналитическая оценка) решения СЛАУ от числа решаемых СЛАУ для разных NA при N = 1000 Вычислительный эксперимент Выполнена программная реализация алгоритмов 1 и 2. Для вычислений использована рабочая станция (без параллельных вычислений, т.е. работало одно ядро процессора) со следующими параметрами: платформа AMD FX(tm)-8350 Eight-Core Processor; частота процессора 4,01 ГГц; объем ОЗУ 32 Гб; число ядер - 8; операционная система Windows 7×64. При вычислительной оценке параметры матрицы и число решений СЛАУ приняты такими же, как при аналитической. Результаты сведены в таблицу 3. Таблица 3. Вычислительные оценки ускорения решения СЛАУ N NA m 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 подтверждают ряд оценок, сделанных выше. Ускорение вычислений возрастает по мере роста m и NA, однако пределы этого роста ниже, чем полученные аналитически, в среднем в 1,5 раза. При этом характер зависимости ускорения от NA и m схож при аналитической и вычислительной оценках. Минимальное ускорение составляет 1,17 раза, максимальное -15,17. Интерес представляет зависимость ускорения от N: из табл. 3 видно, что ускорение возрастает при увеличении N, причем этот рост значителен (до 16%). Подводя итог, можно сказать, что данные, полученные аналитически и в результате вычислительного эксперимента, схожи по характеру зависимости от NA и m, но отличаются и по характеру зависимости от N, и по значениям. Заключение В работе рассмотрен комплекс вопросов по ускорению многократного решения СЛАУ с помощью блочного LU-разложения. В частности, выявлена возможная для использования специфика матриц и проанализирована роль перенумерации как нового этапа, позволяющего, при малых вычислительных затратах, значительное повышение эффективности решения и расширение круга анализируемых структур. Выполнено усовершенствование алгоритма вычисления ряда емкостных матриц с помощью блочного LU-разложения для случая, когда в структуре проводников и диэлектриков изменяются любые параметры. Получены формулы для аналитической оценки ускорения, позволяющие априорную оценку эффективности использования блочного LU-разложения при заданных параметрах решения. Выполнены аналитические и вычислительные оценки эффективности усовершенствованного алгоритма для многократных вычислений, показавшие максимальное ускорение 23 раза. Выявлено, что применение блочного LU-разложения начинает давать ускорение, начиная с довольно малых значений m и не с самых высоких значений NA/N, что перекрывает весьма широкий круг задач на практике. Дополнительно необходимо отметить, что в дальнейшем целесообразны более детальные сравнительные оценки усовершенствованного и исходного алгоритмов. Они позволят расширить сферу применения алгоритмов на основе блочного LU-разложения к другим задачам с частичным изменением матрицы СЛАУ при многократных вычислениях.

About the authors

Roman Sergeevich Surovtsev

Tomsk State University of Control Systems and Radioelectronics

Email: surovsevrs@gmail.com

Sergei Petrovich Kuksenko

Tomsk State University of Control Systems and Radioelectronics

Email: ksergp@mail.ru

Talgat Rashitovich Gazizov

Tomsk State University of Control Systems and Radioelectronics

Email: talgat.tu@tusur.ru

References

Supplementary files