Повышение эффективности решения системы линейных алгебраических уравнений итерационным методом Текст научной статьи по специальности «Математика»

 1926,40 - 500 -

5-Ю-2 49,38 0,14 12 +

10"2 12,74 0,28 3 +

GMRES(20) 510-3 8,73 0,39 2 +

10-3 4,78 0,44 1 4-

5-10“5 5,00 0,68 1

10-4 8,84 4,08 1 н-

5-10-s 17,79 12,25 1 •+-

5-10-2 18,73 0,14 3

GMRES(30) 102 6,70 0,28 1 +

5-Ю-3 6,76 0,39 1 -ь

Из табл. 1 видно, что из всех рассмотренных методов выгоднее использовать итерационные методы с предобуеловливаиием при оптималь-

ном значении параметра т. Так, требуемое решение методом В1СС81аЬ при т = 5- КГ4 было получено в 9 раз быстрее (наибольший выигрыш), чем при использовании метода Гаусса. Использование метода ЕИСС^аЬ позволяет получить результат в 1,6 и 1,4 раза быстрее (при оптимальных значениях т), чем при использовании методов ЕНСС и СМК.Е5(т), соответственно.

Как видно из таблицы, для каждого итерационного метода существует оптимальное значение параметра т. Исследования по выяснению зависимости оптимального значения т от матрицы (структуры объекта), а также пояснения наличия этого факта можно найти в [9, с. 26-28]. В работе [10, с. 54-57] показано, что можно получить ускорение решения СЛАУ итерационными методами за счет снижения точности вычисления (числа итераций) до достаточной для получения заданных характеристик (диаграммы направленности антенны).

Литература

1, Харрингтон Р.Ф. Применение матричных методов к задачам теории поля IIТИЭЭР. 1967. №2.

2, Lee J„ Zhang J., Lu С. Incomplete LU preconditioning for scale dense complex linear systems from electromagnetic wave scattering problems. Computational Physics, 2003. Vol. 185.

3, Saad Y. Iterative methods for sparse linear systems. PWS Publishing Company. 1996.

4, Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. Новосибирск, 2000.

5, Ильин В.П, Методы неполной факторизации для решения линейных систем. М., 1995,

6, Тарап К. Sarkar at al. Survey of Numerical Methods for Solution of Large Systems of Linear Equations for Electromagnetic Field Problems, IEEE Trans, on Antennas and Propagat, Vol. AP-29. Nov, 1981, № 6.

7, http://www.cerfacs.fr/algor/softs/

8, Guennouni A. EL., Jbilou K„ Sadok H, A block version of BICGSTAB for linear systems with multiple right-hand sides. Electronic Transactions on Numerical Analysis. 2003, Vol. 16.

9, Газизов T.P., Куксенко С.П. Оптимизация допуска обнуления при решении СЛАУ итерационными методами с предобуеловливаиием в задачах вычислительной электродинамики II Электромагнитные волны и электронные системы, 2004, № 8,

10, Куксенко С.П., Газизов Т.Р. Ускорение решения СЛАУ в задачах вычислительной электродинамики: Сборник научных трудов VII Всероссийской научно-практической конф, «Проблемы информационной безопасности государства, общества и личности (16-18 февраля 2005)». Томск, 2005.

УДК 512.643:512.644

И.С. Костарев, С.П. Куксенко, Т.Р. Газизов

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

Томский государственный университет систем управления и радиоэлектроники

Введение

Необходимость в решении системы линейных алгебраических уравнений (СЛАУ) возни-

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

или рассеяния электромагнитной волны сложными объектами, являющихся одними из основных задач электромагнитной теории, может быть получено с помощью интегральных уравнений, сводящихся методом моментов (МОМ) к СЛАУ [1, с. 5-19] вида А х х = b, (1)

где А - плотная квадратная матрица размером N х .¡У; b - ненулевой вектор размером N; ж- вектор неизвестных (вектор решений).

11ри компьютерном моделировании основные вычислительные затраты состоят из суммы затрат на формирование матрицы и затрат непосредственно на решение СЛАУ. Следовательно, выбор наиболее эффективного способа решения СЛАУ позволит снизить общие временные затраты.

Решать СЛАУ можно прямыми (точными) и итерационными методами, каждые из которых характеризуются oiфеделенным быстродействием. При формальном подходе к решению системы (1) с N неизвестными (например, используя один из прямых методов - метод Гаусса) число арифметических операций ~А?3, что при большом числе неизвестных затрудняет решение системы и ограничивает крут рассматриваемых задач. Поэтому в последнее время широко используются для решения СЛАУ итерационные методы, так как для их сходимости требуется число итераций много меньше, чем N, причем основные вычислительные затраты ~N" приходятся на умножение матрицы на вектор в каждой итерации. Таким образом, проблема увеличения быстродействия требует определенного подхода, связанного уже непосредственно е увеличением скорости перемножения матрицы на вектор.

Сами алгоритмы, позволяющие быстро перемножить матрицу на вектор, были разработаны сравнительно недавно. О том, что данное направление является новым и далеко не полностью изученным, говорит тот фа кг, что основными источниками, которые удалось обнаружить, являются статьи [2, р. 177-185; 3-5; 6, р. 280-292; 1, р. 553-591; 8, р, 48-51; 9; 10, р. 177-185]. Причем в рассмотренных автором публикациях не представлены результаты вычислительных экспериментов, позволяющие реально оценить возможности того или иного алгоритма.

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

Метод сопряженных градиентов с использованием быстрого преобразования

Фурье

Метод сопряженных градиентов с использованием быстрого преобразования Фурье (СО-РРТ) был предложен в [2, р. 177-185]. Он является одним из первых алгоритмов, позволяющих увеличить скорость перемножения плотной матрицы на вектор. Первоначально этот алгоритм использовался для вычисления плотности тока в плоском проводнике от г-компоненты магнитного поля, причем, как утверждают авторы, применение его возможно для проводников любой толщины при условии, что вектор токов имеет только х- и у-компоненты. Как видно из самого названия алгоритма, он сочетает в себе два метода: - метод сопряженных градиентов и мегод быстрого преобразования Фурье (БПФ). Использование БПФ было направлено на непосредственное увеличение скорости перемножения матрицы на вектор. Основой для применения БПФ являлось то, что в методе СО матрица А являегся структурной, в частности Тепл и девой. Теплицева матрица - это одна из структурных матриц [3] вида (рис. 1).

' а0 ал ... ®2~N ai-N

а] ciQ <V,Y

а, ...

V ^.V - і «V - 2 а. с7р ^

Рис. 1. Теплицева матрица

Как утверждается авторами, благодаря такому виду матрицы и применению метода БПФ для умножения её на произвольный вектор можно уменьшить временные затраты на перемножение матрицы на вектор, с ~М2 до -/V \0g2N.

Метод СО-РРТ был рассмотрен и в работе [4], где он используется для решения проблем рассеяния электромагнитных волн от объемных объектов. Матрично-векторное перемножение представляет собой действие оператора Грина на ток, наводимый в объекте (рассеивателе), что можно записать как интеграл свертки

| сіг' ■ ‘¿(г ~ г') ■ і (г') . где g(r) - функция Грина V

и/(/-)- наводимый ток.

Данный интеграл представляет собой свертку, и согласно [11] применение БПФ при вычислении свертки может свести вычислительные затраты до ~/V log2jV. Также этот интеграл можно продискретизировать и, используя МОМ, свести его к СЛАУ. Причем полученная матрица А будет иметь структуру типа Теплиц, что. в свою очередь, позволит увеличить скорость перемножения матрицы на вектор за счет использования БПФ.

Наравне с Теплицевыми матрицами используют и другие структурные матрицы, такие как матрицы Ваидермонда (V), Ханкеля (Н) [3]. Вид

этих матриц представлен на рис. 2.

Н

( а а \

- N і S а - N - 2 и ■ N : 3 “о

а . д. й.д,. 3 a-Ni 4 а,

а . N тЗ °-,v. 4 а.

aN.~2

\ ао а! Ор!_ J <V. J

Рис. 2, э. Матрица Ханкеля

' 1 1 1 ... 1 ^

а0 а, а., 4>V-J

2 2 -> 2

а0 0\ аг aN-l

N-1 N 1 ал («Л'.2)Л! ( V.) J

Рис. 2, б. Матрица Ваидермонда

Необходимо отметить, что использование БПФ для умножения плотной матрицы на вектор применяется в сочетании с другими итерационными методами, например, с методами BiCG, BiCGSTAB(L) [5].

Быстрый метод многополюсника

Впервые быстрый метод многополюсника (Fast Multipol Method - FM.M) был предложен в 1987 году [6, р. 280-292] и рассматривался для оценки потенциала в системах с большим количеством частиц, взаимодействие которых имеет гравитационный характер. Идея алгоритма заключается в том, что частицы, находящиеся в определенной области пространства, объединяют в кластеры и вычисляют их (кнастеров) взаимодействие с другими такими же кластерами, тем самым уменьшая общее количество операций.

В настоящее время данный алгоритм используется в разных областях и направлениях науки [7. р. 553-591; 8, р. 48-51]. В частности, он также применяется и для увеличения скорости умножения плотной матрицы на вектор в итерационных методах [4; 9]. Детальное описание метода Р.ММ можно найти в вышеупомянутых ссылках. Далее этот метод рассматривается в приложении к перемножению плотной матрицы па вектор [4].

При умножении плотной матрицы на вектор число операций достигает ІУ2. Такой алгоритм легко представить в виде одноуровневой системы (рис. 3).

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

Ну j)

N

a"'°v

J =

М.

(2)

// ■•У ,

і/У // Vх \\ \Ч ) ■■

* У -- ч : -

Рис. 3, Алгоритм стандартного перемножения матрицы на вектор

Эта сумма может быть представлена следующим образом:

. І

Ум.

(.У, -х,)

О'"’] -^.у)

(.Ум ~х\)2 ■■■ (Ум _ХА')2

В тоже время сумму (2) можно записать как:

игхг -

2 • У і &• хі = а • >’/ + Ь ~ 2 • У.; • 8>

_/=1,..., М,

где коэффициенты а, b и g получаются за C)(N) операций умножения (всего за 3N операций); суммы v(y¡) получаются за 3М операций. Следовательно, общее число операций будет 0(3 М + ЗЛО = 0(М + N). Таким образом, идея РММ алгоритма заключается в изменении порядка суммирования при условии, что элементы матрицы получены из специфической функции, которая может быть разложена на множители или на переменные.

Необходимо отметить, что алгоритм. РММ работает как с двумерными объектами, так и с трехмерными. Отличие заключается в том, что при работе с двумерными объектами FMM уменьшает сложность матрично-векторного умножения и требуемый объем памяти до 0(N), а с трехмерными объектами до 0(N4fi).

Многоуровневый быстрый алгоритм многополюсника

Многоуровневый быстрый алгоритм многополюсника. (Multilevel Fast Multipole Algoithm -MLFMA) [10, p. 177-185] был разработан на основе FMM и применяется для решения задач, связанных с рассеянием электромагнитных волн от объемных объектов. Для реализации MLFMA необходимо поместить некий объект в куб, который .целится на восемь меньших кубов. Каждый полученный куб, в свою очередь, рекурсивно делится на еще меньшие кубы до тех пор, пока длина ребра одного куба не будет равна примерно половине длины волны. Таким образом, получают структуру, состоящую из множества уровней, вложенных один в другой (которые состоят из кубов от самых больших и до самых маленьких). Далее из всех кубов выбирают только «непустые» кубы, которые в дальнейшем кодируются е использованием древовидной структуры данных [12, р. 923-947]. Таким образом получается, что вычислительные затраты .зависят только от непустых кубов.

Использование метода MLFMA для матрич-но-векторного умножения приводит к тому, что данное действие происходит в две стадии [12]. Первая стадия заключается в расширении многополюсника для каждого непустого куба на всех уровнях. Такое расширение называется внешним расширением многополюсника и используется для того, чтобы вычислить область вне куба. Вторая стадия заключается в построении локального расширения многополюсника от отдаленных кубов на всех уровнях. Вычисляя последовательно внешнее и локальное расширение многополюсника, можно уменьшить временные

затраты на перемножение матрицы на вектор для трехмерных объектов до N 10g2iV.

Особое внимание хотелось бы уделить обзорной, но емкой работе [4]. В ней представлено несколько быстрых методов решения интегральных уравнений (в том числе и рассмотренных выше), применяемых к поверхностным и объемным рассеивателям. Эти методы уменьшают время вычисления и объемы памяти, в частности и за счет уменьшения затрат, связанных с умножением плотной матрицы на вектор. Чтобы не повторяться, рассмотрим только часть из представленных в статье методов.

Матричный декомпозиционный алгоритм

и многоуровневый матричный декомпозиционный алгоритм

Матричный декомпозиционный и многоуровневый матричный декомпозиционный алгоритмы (The Matrix Decomposition Algorithm, Multilevel Matrix Decomposition Algorithm -MDA, MLMDA) позволяют увеличить скорость решения СЛАУ итерационными методами при решении задач рассеяния электромагнитных волн структурными объектами. .В отличие от FMM, MDA и MLMDA анализируют матрицы МОМ с использованием линейных алгебраических методов. MDA и MLMDA используют ограниченное число степеней свободы, которые характеризуют поле рассматриваемой области, находящейся на «большом расстоянии» от исходной области.

В MDA недиагональные блоки матрицы метода моментов соединены в большие блоки и анализируются с использованием многоуровневого алгоритма, который напоминает FFT, как показано на рис. 4. Требуемые объемы памяти компьютера, а так же вычислительная сложность MDA - 0(Nxl1), для MLMDA приближается к <9(i¥Iog22/V).

Z*

Рис. 4, а. Матрица МОМ

Рис. 4, б. Матрица МОМ после использования MDA

Алгоритм наискорейшего спуска

Алгоритм наискорейшего спуска (Fast Steepest Descent Path Algorithm - FASDPA) представляет собой новый двухуровневый алгоритм, совмещающий в себе два алгоритма FMM и MDA, Мало чем отличаясь от них, FASDPA начинается с пространственно]'« разложения рассеивателя на большое количество подрассеивателей, где взаимодействия между близлежащими подрас-сеивателями рассчитываются прямыми методами. Взаимодействия между отдаленными под-рассеивателями рассчитываются, как это происходит в MDA. Вычислительная сложность FASDPA алгоритма, в общем случае, достигает 0(М4ГЗ).

Алгоритм перемножения матрицы на вектор с использованием БПФ

И.з всех рассмотренных, выше алгоритмов для реализации был выбран алгоритм перемножения матрицы на вектор с использованием БПФ. Рассмотрим его подробнее. Пусть дана матрица типа Теплица. Тогда для вычисления произведения матрицы А на вектор ж за время 0(A/log2.A/) необходимо провести следующие операции:

1. Методом БПФ вычисляем преобразование Фурье от вектора а - [щ, аь ..., сік-\]

(где вектор а - вектор-столбец матрицы А) и вектора х = [а'о, х\, Xv-il размерами 2N по формулам:

f = FN(x) и g = Fn(&),

где F,.v - Фурье преобразование.

2. Вычисляем произведение f на g, обозначив его z

Z = [Z(b Z], ..., Z2N- lj. Zi

3. Вычисляем обратное преобразование Фурье от вектора %

У = Р*к(ъ).

Полученный результат, вектор у, и является искомым произведением матрицы на вектор. Таким образом, видно, что алгоритм состоит из вычисления двух БПФ размера 2Д;, 2N произведений базовых типов данных (их сложность пренебрежимо мала) и одного обратного БПФ размера 2А.

Так же на основе представленного выше алгоритма был реализован метод Шенхаге-Штрассена для перемножения длинных чисел. Метод отличается но своему построению от метода перемножения матрицы на вектор с использованием БПФ тем. что сначала числа, которые необходимо перемножить, представляются в виде полинома и записываются в вектор, начиная с элемента, обладающего наименьшим порядком. Далее полученный вектор дополняется нулями до тех пор, пока число элементов не будет соответствовать N = 2". После этого все пункты представленного выше алгоритма повторяются. Метод Шенхаге-Штрассена использовался для увеличения непосредственного процесса перемножения элементов при умножении матрицы на вектор [13, с. 1 09-1 11].

Сравнительный анализ

Алгоритмы перемножения матрицы на вектор с использованием БПФ были реализованы в среде МаАсас! и на языке С-І-+. Сравнительные результаты времени перемножения матрицы на вектор традиционным методом и методами с использованием БПФ (реализованные на языке С++) показаны в табл. 1 (анализ проводился на компьютере со следующими параметрами: процессор - 600 МГц, оперативная память -128 Мбайт).

Таблица 1 Результаты вычислений при N 512

Метод перемножения матрицы на вектор Время, ыс

Традиционный 25

Метод БПФ (для N = 2'', ¥> 0) 10

Метод БПФ (для N четного) 12

Метод Шенхаге-Штрассена 24.73

Из табл. 1 видно, что выгоднее оказывается метод непосредственного перемножения матрицы на вектор с использованием БПФ для N = Т (V > 0). Причем выигрыш данного метода по сравнению с традиционным методом перемножения матрицы на вектор составляет 2.5 раза.

Вычислительное экспериментирование

Было проведено сравнение по быстродействию метода перемножения матрицы на вектор с использованием БПФ (для .¥=2’' (т>0)) с традиционным методом перемножения дом различных размеров матриц [14, с. 112-114]. Результаты эксперимента представлены в табл. 2 и на рис. 5.

Зависимости времени перемножения матрицы на вектор от порядка матрицы представлены на рис. 5.

Из табл. 1, а так же рис. 5 видно, что использование алгоритма БПФ позволяет ускорить процесс перемножения матрицы на вектор по сравнению с традиционным методом перемножения. Причем следует отметить, что с ростом порядка матрицы N выигрыш растет, и п р и

Таблица 2

Зависимость времени перемножения матрицы на вектор от размера, матриц

Размеры, используемой матрицы A Время перемножения традиционным методом, И, мс Время, перемножения с использованием БПФ, 12, мс. Отношение временных результатов методов 11/12

32x32 0,03 0,058 0,517

64x64 0.12 0.16 0.75

128x128 0.92 1.2 0.78

256x256 3.7 3.2 1.16

512x512 25 10 2.5

1024x1024 110 25 4.4

2048x2048 1600 58 27

■Т. мс-

MV

MVFFT

Рис, 5. Зависимость времени перемножения матрицы на вектор от размера матрицы. МУ - традиционный метод матрично-векторного перемножения;

МУРРТ - метод перемножения'матрицы на вектор с использованием БПФ

Литература

1, Харрингтон Р.Ф, Применение матричных методов к задачам теории, поля IIТИЭЭР. 196?, № 2.

2, Rinke J, Wijngaarden. Fast determination of 2D current patterns in flat conductors from measurement of their magnetic field - Physica 295, 1998,

3, Olshevsky V,, Shokrollahi A. A Superfast Algorithm for Confluent Rational Tangential Interpolation Problem via Matrix-vector Multiplication for Confluent Cauchy-like Matrices.

4, Cui T.J., Chew W.C. Fast Algorithm for Electromagnetic Scattering by Buried Conducting Plates of Large Size—IEEE Trans, on Antennas and Propagat. Vol. 47. No. 6. June, 1999.

5, Fang S., Torres-VerdHn С., Guo-zhong G. A New Approximation for 3D Electromagnetic Scattering in the Presence of Anisotropic Conductive Media -3DEMII Workshop. February. 2003,

8. Greengard L, Rokhlin V. A Fast Algorithm for Particle Simulations - Journal of Computational Physics 135,1997.

7. Labreutche C, A Convergence theorem for the multipole method for 2 dimesional scattering problems, Matemating of Computation. Vol. 67. No 222. April, 1998.

8. Colfman R., Rokhlin V. The Fast Multipole Method for Electromagnetic Scatering Calculations - IEEE Trans, on Antennas and Propagat № 4. Januare. 1993,

9. Saad Y. Iterative methods for sparse linear systems. PWS Publishing Company, 1996.

10. Cal-Cheng Lu, Chew W.C., Multilevel Fast Multipole Algorithm for Electromagnetic Scattering by Large Complex Objects, IEEE Trans, on Antennas and Propagat, Vol. 45. №. 10. October, 1997,1998,

11. Гольденберг Л.М. Цифровая обработка сигнала. М., 1987.

12. Anderson C.R, “An implementation of the fast multipole method without multipoles,” SIAM J. Sci. Stat, Comput. Vol. 13. №. 4. Apr. 1992.

13. Костарев И.С. Решение СЛАУ итерационными методами II Материалы Всероссийской научно-технической конференции студентов, аспирантов и моподых специалистов «Научная сессия ТУСУР-2005» (26-28 апреля 2005 г.). Томск, 2005.

14. Костарев И.С., Куксенко С.П. Увеличение скорости решения СЛАУ с помощью быстрого преобразования Фурье // Материалы Всероссийской научно-технической конференции студентов, аспирантов и молодых специалистов «Научная сессия ТУСУР-2005» (26-28 апреля 2005 г.). Томск, 2005.

УДК 535.3:621.391.63

МЛ. Фролова

ПЛАНАРНЫЕ ЧЕТЫРЕХСЛОЙНЫЕ СТРУКТУРЫ» СОДЕРЖАЩИЕ ВОЛНОВОД С ПА-РАБОЛИЧЕСКИМ ПРОФИЛЕМ ПОКАЗАТЕЛЯ ПРЕЛОМЛЕНИЯ.

Томский государственный университет систем управления и радиоэлектроники

В настоящее время важную роль в устройствах интегральной оптики играют многослойные структуры, использующие волноводное распространение света [1, р. 1397-1413; 2, с. 1555-1559;

3, с. 305-310]. Многие технологические методы получения волноводов - диффузия примесей, ионный обмен в кристаллах ниобата и танталата лития, а также жидкофазная эпитаксия для кристаллов силленитов - приводят к формированию градиентных волноводов, характеризующихся плавным изменением показателя преломления по их толщине [2-3; 4, с. 103-108; 5, с. 80-82]. Профиль показателя преломления, то есть закон изменения Дп(г), является одной из важнейших, характеристик оптических волноводов, поскольку он определяет их модовый состав и распределение светового поля, которые должны приниматься во внимание при проектировании интегрально-оптических схем и их элементов. Важное преимущество использования технологий получения градиентных и многослойных волноводов в оптических устройствах состоит в том, что имеется возможность управления параметрами волноводных и изолирующих слоев. Так, например, технология получения кристаллов со структурой силленитов и их слоев путем изменения состава позволяет варьировать фундаментальные оптические и хироптические

характеристики: спектры отражения, дисперсию показателя преломления, оптическое вращение [6, с. 356-380].

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

Некоторые технологические методы получения волноводов, такие как нанесение эпитаксиальных слоев В^ТЮзо методами кристаллизации из раствора в расплаве и кристаллизацией из ограниченного объема, описанные в [6], приводят к тому, что новый кристаллический слой с показателем преломления п() формируется не сразу. На границе с подложкой возникают переходные слои с промежуточными параметрами решетки. Эти переходные слои имеют плавный профиль показателя преломления, который может быть аппроксимирован асимметричной параболической функцией (рис. 1, о) или комбинацией прямого участка и параболы (рис. 1, б), в зависимости от результатов операции шлифования поверхности. Параболическое распределение показателя преломления является наиболее подходящей функцией для этих целей, поскольку для нее имеется аналитическое решение волнового уравнения, кроме того, ее параметры легко изменять, вплоть до ступенчатой