Псевдоскелетные аппроксимации для блочных матриц, порожденных асимптотически гладкими ядрами тема диссертации и автореферата по ВАК РФ 01.01.07, кандидат физико-математических наук Горейнов, Сергей Анатольевич

• Горейнов, Сергей Анатольевич
• кандидат физико-математических наук
• 2001

Оглавление диссертации кандидат физико-математических наук Горейнов, Сергей Анатольевич

Введение.

Глава 1. Мозаично-скелетонный метод.

1.1 Аппроксимация для асиптотически гладких ядер

1.2 Оценки асимптотической гладкости.

1.3 Аппроксимация для гармонических ядер.

1.4 Аппроксимация для осцилляционных ядер.

Глава 2. Алгоритмы аппроксимации сложности 0(п2)

2.1 Аппроксимация блока по методу Ланцоша.

2.2 Описание алгоритма.

2.3 Численные эксперименты.

2.3.1 Уравнение электрического поля.

2.3.2 Задача обтекания.

Глава 3. Псевдоскелетные аппроксимации.

3.1 Возможность псевдоскелетной аппроксимации.

3.1.1 Экстремальные подматрицы.

3.1.2 Оценки псевдоскелетной аппроксимации.

3.2 Аппроксимации при помощи подматрицы максимального объема.

3.3 Принцип выбора псевдоскелетной компоненты.

3.4 Численные эксперименты.

Введение диссертации (часть автореферата) на тему «Псевдоскелетные аппроксимации для блочных матриц, порожденных асимптотически гладкими ядрами»

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

К(х,-у)и(у) сЯу = ^х), которое каким-то образом дискретизируется (скажем, коллокацией или по методу Галеркина), и возникает линейная система Апхп = Ъп, размерность п которой растет при увеличении точности дискретизации. Ключевой вопрос — как быстро растет объем памяти, необходимый для хранения дискретизированного оператора, и число арифметических действий, необходимое для решения задачи с заданной точностью, в зависимости от п. Предполагается, что задача поставлена корректно и ядро К(х,у) сингулярно.

Предположим, что для решения линейной системы используется какой-либо итерационный метод (возможно, с предобусловлива-нием). Тогда основные вычислительные затраты связаны с операцией матрично-векторного умножения на Аи, что приводит к следующей задаче матричной аппроксимации, составляющей содежание проблемы: найти „близкую" к Ап матрицу Ап, такую, что

• матрично-векторное умножение у = Апх выполняется быстро;

• объем памяти, необходимый для хранения Ап, мал;

• погрешность в решении, допускаемая при замене Ап на Ап, сравнима с погрешностью дискретизации.

Говоря „быстро" и „мал" мы имеем в виду о(п2) арифметических операций и о(п2) ячеек памяти. Сложнее указать строгий критерий близости Аи и Аи, обеспечивающий допустимую погрешность в решении. В дальнейшем мы будем понимать под этим условие || Ап-АпЦр ^ £ 11 Ап 11 р, где £ не зависит от п, потому что для тех конкрет 5 — ных уравнений, с которыми мы имели дело, этого (при надлежащем выборе е) в самом деле достаточно.

Нельзя здесь не упомянуть следующую гипотезу, принадлежащую Н.С. Бахвалову [2] й относящуюся к корректно поставленным краевым задачам математической физики:

Пусть известно, что правая часть f(x) принадлежит некоторому компактному множеству, a N£ — минимальное число значений ft, достаточное для того, чтобы при любой правой части f по этим значениям можно было восстановить с точностью £ решение и(х). Тогда для решения корректной задачи математической физики достаточно использовать 0(Ne) значений и дополнительно произвести 0(N£logaN£log^ арифметических операций.

Первые шаги в получении быстрых приближенных методов решения интегральных уравнений, согласованные с этой гипотезой, были сделаны в 60-е годы [11, 13]. Нас, однако, будет интересовать только часть этой задачи — именно, быстрое приближенное умножение на интегральный оператор.

За последние 10-15 лет здесь появилось довольно много подходов (см., например, [20]), как правило, для какого-либо конкретного ядра (log |х — -у 1/|х — е1К'х-у1/|х — у |) и с помощью какого-либо специального метода — мультипольных разложений [52,35,53,54,47,21], кластеризации граничных элементов [36], интерполяции на регулярную (иерархическую) сетку [14,27], использования „многоуровневого анализа" и локальных волн (multiresolution analysis, wavelets) [25,37,49]. С алгебраической точки зрения, большинство этих подходов явно или неявно базируются на одном и том же наблюдении: для весьма широкого класса ядер К(х,-у) и для всех достаточно больших п можно указать такое блочное разбиение матриц Ап, что подавляющая часть блоков будет иметь хорошее малоранговое приближение. Целью настоящей работы является исследование и использование этой „алгебраической" идеи для получения эффективных алгоритмов приближенного матрично-векторного умжножения для достаточно широкого класса порождающих ядер. Именно эта идея напрямую исследуется в [62] в случае, когда ядро является асимптотически гладкой функцией. Так, вслед за Брандтом [26], называют функцию К(х,у), если d;k(tc,u)| ^ срЦх-цЦв-p, . vp^o. (o.i)

Здесь DJ обозначает оператор дифференцирования по у порядка р; ср — положительная константа, зависящая от р; g — некоторая константа; норма || • || — евклидова.

В этом случае основным инструментом для получения оценок на е-ранг блока служит формула Тейлора: в общем члене переменные х и у разделены, поэтому ранг суммы таких членов легко оценивается и не зависит от n, а остаточный член с учетом (0.1) определяет погрешность аппроксимации. Отметим, что теорема об аппроксимации всей матрицы следует из такого рода оценок только в том случае, когда делаются какие-либо предположения о характере зависимости ср от р — не слишком обременительные, но приводящие к приемлемому результату. В настоящей работе используется определение (1.1), принадлежащее Е. Б. Тыртышникову [62].

Главный результат работы [62] выглядит следующим образом. Пусть S С Ет — ограниченная область, в которой заданы наборы узлов (xii)}^ и {у ())}]!=!, удовлетворяющие условию квазиравномерности, то есть существуют константы Ti, Т2, такие, что для всякого m-мерного куба О С S, содержащего по крайней мере две точки сетки, имеет место следующее неравенство: mes Q , ^, mes О. Ti —71 ^ И-(^) ^ т2-— п. mes S mes S

Символом ц(П) здесь обозначено число узлов, принадлежащих области С1. Пусть, далее, К(х,у) — асимптотически гладкая функция. Будем говорить, что К(х,у) порождает матрицы Ап, если для любого п справедлива формула Ап = [К(х(г),у(Л)]у=1. Тогда для всякой Ап, порожденной ядром К, существует мозаично-скелетонная аппроксимация Ап, такая, что

An — ÂU||F = О(пе), (0.2)

Лп занимает 0(nlog nlogm 1 /е) ячеек памяти и за столько же арифметических операций может быть умножена на вектор. Символ О относится здесь к пределу п —» оо, а £ — произвольное достаточно малое положительное число.

Несколько изменив технику доказательства, автору удалось получить такую же теорему при ослабленном требовании к сеткам; достаточно, чтобы

ТТ"| рс С~)

4П) ^ т-^-п, vacs. (0.3) mes S

Это означает, что в качестве S можно выбирать, например, многообразие меньшей размерности, чем m, с сохранением указанных оценок аппроксимации.

Впрочем, в некоторых случаях показатель степени в члене logm 1 /е можно уменьшить до размерности многообразия S. В настоящей работе доказан подобный результат для случая, когда S является простым несамопересекающимся кусочно-целым контуром. Имеется в виду следующее: пусть функция y(v) : [0, 1] i—> (С, параметризующая контур, такова, что для каждого интервала какого-либо конечного разбиения отрезка [0, 1] ^Hy(v) и 3y(v) представляются рядами Тейлора с бесконечным радиусом сходимости, имеют отделенные от нуля производные и при v —> оо растут не быстрее экспоненты. Тогда для матриц, порожденных асимптотически гладкой по у функцией К(х,у), справедливы приведенные выше оценки с m = 1. Эти результаты составляют содержание раздела 1.1. В разделе 1.2 исследована асимптотическая гладкость некоторых конкретных функций; кроме того, доказано, что всякая гармоническая функция является асимптотически гладкой.

В разделах 1.3 и 1.4 для получения оценок типа (0.2) в случае неко 8 — торых конкретных ядер К(х,'у) используются разложения, отличные от формулы Тейлора. В разделе 1.3 гармонические по переменному у ядра анализируются при помощи производящей функции для ультрасферических многочленов; в разделе 1.4 мозаично-скелетонные аппроксимации для фундаментальных решений уравнения Гельмгольца в двумерном и трехмерном случае строятся при помощи теорем сложения для цилиндрических функций. В целях замкнутости изложения некоторые классические результаты, используемые в работе, приведены с доказательствами.

Глава 2 посвящена численной проверке работоспособности мозаич-но-скелетонного метода, изложенного в главе 1, на довольно сложной трехмерной задаче, связанной с рассеянием монохроматической электромагнитной волны на идеально проводящем экране. При этом, в качестве первого шага, позволяется использовать 0(п2) арифметических операций при построении приближающей матрицы. Аппрокси-мационный алгоритм, используемый для блоков матрицы, никак не зависит от порождающего ядра: фактически, это метод Ланцоша для симметризованного блока. В разделе 2.1 исследуется сходимость метода в терминах задачи матричной аппроксимации (использовать хорошо известную теорию Каниэля-Пэжа-Саада [48,55] в наших целях оказывается сложнее, чем получить прямые оценки). Используемый алгоритм детально описывается в разделе 2.2. Собственно численные результаты составляют содержание раздела 2.3.

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

А^Св^ С=А"\ (0.4) 9 — где А — подматрица Л порядка г, находящаяся на пересечении строк И и столбцов С. Разложение (0.4) называется иногда скелетным. Оказывается, что для матрицы А = А + Е также существуют матрицы С и И (при этом С, например, содержит г столбцов А, а не А!) такие, что

А-СС1*||2^(т,т)£, е = ЦЕЦ2, (0.5) где С — некоторая матрица порядка г, а f — слаборастущая функция гит. Матрицу ССК будем называть псевдоскелетной компонентой А. Получение оценок такого вида происходит в разделе 3.1.

Теперь, если до вычисления матричных элементов, составляющих заданный блок, знать номера строк и столбцов, дающих „хорошую" псевдоскелетную компоненту, можно вычислить только С и И, избавившись тем самым от асимптотики 0(тг2) по арифметическим затратам на построение всей приближающей матрицы.

Раздел 3.2 посвящен характеризации строк и столбцов, составляющих псевдоскелетные компоненты. Именно, если в качестве А выбрать подматрицу А, для которой модуль детерминанта максимален среди всех подматриц А порядка г, имеет место оценка типа (0.5). Наличие такой характеризации позволит, как кажется автору, применять псевдоскелетные аппроксимации для достаточно широкого класса ядер. Что это возможно в простых случаях — например, для контурного уравнения с логарифмическим ядром, — видно из численных результатов для этого уравнения, приведенных в разделе 3.4. Принцип выбора псевдоскелетной компоненты, реализованной в соответствующем алгоритме аппроксимации, описан в разделе 3.3.

Перейдем к формулировке основных результатов.

Теорема 1.1. Пусть матрица Ап порождается асимптотически гладким ядром К(х,у), е Мт, на сетках в Мт, подчиненных условию (0.3); тогда существует матрица Ап такая, что ||АП—АП||Р = 0(п7е), Ап занимает 0(пк^пк^т 1) ячеек памяти и за столько же арифметических операций может быть умножена на вектор.

10

Здесь и далее, у ^ 1 — константа, зависящая от свойств ядра К; явные формулы для у имеются в основном тексте работы.

Теорема 1.2. Пусть матрица Ап порождается асимптотически гладким ядром К(х,-у), х,у € К2, на сетках, заданных на кусочно-целой кривой Б, имеющей } точек разрыва. Пусть сетки подчинены условию тлея Г

УГсБ, mes Г mes S где Г — пересечение S с любым квадратом, содержащее по крайней мере две точки сетки. Тогда существует матрица Ап такая, что ~ ~ 1 ||АП — Лп||р = (9(J3/2nr£), Ап занимает OQnlognlogячеек памяти и за столько же арифметических операций может быть умножена на вектор.

Теорема 1.3. Пусть функция К(х,у), х,у G Rm, — гармоническая по переменному у. Тогда К(х,-у) — асимптотически гладкая функция.

Теорема 1.4. Пусть матрица Ап порождается гармоническим по у ядром К(х,у), х,у е Жт, на сетках в Мт, подчиненных условию (0.3); тогда существует матрица Ân такая, что ||An — Ân||F = 0(пте), Ап занимает O(nlognlogm11) ячеек памяти и за столько же арифметических операций может быть умножена на вектор.

Теорема 1.5. Пусть х,у G Mm, m = 2,3, к > 0, и

Но1)(к|х--у|), m = 2;

К(х,у) = гк|х-у l/\x-y\} ТТ1 = 3.

Пусть матрица Ап порождается ядром К(х,у) на сетках, подчиненных условию (0.3) и принадлежащих ограниченному множеству диаметра а; тогда существует матрица Ап такая, что ||АП — Ап||р = 0(тг£), А занимает ©(пк^тг^к, ос, е)) ячеек памяти и за столько же арифметических операций может быть умножена на вектор. Символом f обозначена функция n m fm(K, ex, е) = ка + logj)2, m = 3.

11 —

Теорема 2.4. Пусть все проекции фг стартового вектора Ланцоша qi на собственные векторы матрицы АТА ненулевые и все собственные значения At матрицы АТА различны. Тогда погрешность малоранговой аппроксимации Ак, порождаемой процессом Ланцоша для А на шаге 1с, оценивается по формуле М где — погрешность наилучшей к-ранговой аппроксимации А по фробениусовой норме; т^)Тг — чебышевское уклонение на множестве

Ат>Aji] U [Aj+i,An] для многочлена степени к, принимающего в точке

Aj значение 1.

Теорема 3.1. Пусть A, F G Cmxn, rank(A-F) ^ к, и ||F||2 ^ £ для некоторого £ > 0. Тогда в А существует подматрица А порядка к такая, что опирающаяся на нее псевдоскелетная компонента CGR для некоторой матрицы G допускает оценку

А — CGR||2 ^ const.£>/kmax(m,n).

Теорема 3.2. Пусть выполнены условия теоремы 3.1. Тогда в А существует подматрица А порядка к такая, что опирающаяся на нее псевдоскелетная компонента CA^R допускает оценку А - CAJ.R||2 ^ const.£kVmn.

Символом А{ обозначена т-псевдообратная к матрице А.

Теорема 3.4. Пусть выполнены условия теоремы 3.1 и, кроме того, подматрица А невырождена и имеет максимальный по модулю детерминант среди всех подматриц А порядка к. Тогда опирающаяся на А

Л -J псевдоскелетная компонента СА R допускает оценку

A — cA-1R||c ^ (к+ 1)£.

Следствие теоремы 3.5. Пусть выполнены условия теоремы 3.1 и, кроме того, подматрица А порядка к невырождена и с!е1 (А)| где у™ах(А) — максимум модуля детерминанта подматриц Л порядка к, а^, 1 ^ -V > 0 — фиксированное число. Тогда опирающаяся на Л псевдоскелетная компонента С А-1 К допускает оценку

ЦА-СА-'ИЦс^^г.

Все перечисленные результаты являются новыми и, за исключением теорем 3.1 и 3.2, принадлежат автору; последние получены в соавторстве с Е. Е. Тыртышниковым и Н. Л. Замарашкиным [9].

Постановка задачи принадлежит Е. Е. Тыртышникову [20]. Результаты настоящей работы следует рассматривать как обобщение подхода, указанного в [62].

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

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

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

1. Н. И. Ахиезер. Классическая проблема моментов и некоторые вопросы анализа, связанные с нею. М.: Физматлит, 1961.

2. Н.С.Бахвалов. Об оптимальных методах решения задач. Appl. Mat. 13: 27-43 (1969).

3. В. В. Воеводин. Линейная алгебра. М.: Наука, 1980.

4. Д. Гайер. Лекции по .теории аппроксимации в комплексной плоскости. М.: Мир, 1986.

5. Ф. Р. Гантмахер. Теория матриц. М.: Наука, 1988.

6. Я. JI. Геронимус. Теория ортогональных многочленов. M.-J1.: ГИТТЛ, 1950.

7. Дж. Голуб, Ч. Ван Лоун. Матричные вычисления. М.: Мир, 1999.

8. С. А. Горейнов. Мозаично-скелетонные аппроксимации матриц, порожденных асимптотически гладкими и осцилляционными ядрами. Матричные методы и вычисления. — под ред. Е. Е. Тыртышнико-ва. с. 42-76. М.: ИВМ РАН, 1999.

9. С. А. Горейнов, Е. Е. Тыртышников, Н. Л. Замарашкин. Псевдоскелетные аппроксимации матриц. Доклады РАН 343(2): 151-152, 1995.

10. С.А.Горейнов, Е.Е.Тыртышников, Н.Л.Замарашкин. Псевдоскелетные аппроксимации при помощи подматриц наибольшего объема. Мат. заметки 62(4): 619-623 (1997).

11. К.В.Емельянов, А.М.Ильин. О числе арияметических действий, необходимом для приближенного решения интегрального уравнения Фредгольма II рода. ЖВМ и МФ 7(4): 905-910 (1967).

12. Г. И. Марчук, Ю. А. Кузнецов. Некоторые аспекты итерационных методов. Вычислительные методы линейной алгебры, с. 4-20. — под ред. Г. И. Марчука. Новосибирск, 1972.

13. В.И.Лебедев. Об итерационном КР-методе. ЖВМ и МФ 7(6): 1250-1269 (1967).87 —

14. Ю. М. Нечепуренко. Быстрые,численно устойчивые алгоритмы для широкого класса линейных дискретных преобразований. — М.: 1985. Препринт ОВМ АН СССР №98.

15. А. Ф. Никифоров, В. Б. Уваров. Основы теории специальных функций. М.: Наука, 1974.

16. Б. Парлетт. Симметричная проблема собственных значений. М.: Мир, 1983.

17. Г. Сегё. Ортогональные многочлены. М.: Физматлит, 1962.

18. Ю. Г. Смирнов. О разрешимости векторных интегродифференци-альных уравнений в задаче дифракции электромагнитного поля на экранах произвольной формы. ЖВМ и МФ. 34:1461-1475 (1994).

19. П. К. Суетин. Ряды по многочленам Фабера. М.: Наука, 1984.

20. Е. Е. Тыртышников. Методы быстрого умножения и решение уравнений. Матричные методы и вычисления. — под ред. Е. Е. Тыр-тышникова. с. 4-41. М.: ИВМ РАН, 1999.

21. Дж. Уилкинсон, К. Райнш. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра. М.: Машиностроение, 1976.

22. С. R. Anderson. An implementation of the fast multipole method without multipoles. SI AM J. Sci. Stat. Comput. 13(4): 923-947 (1992).

23. E. Anderson, Z. Bai et al. LAPACK Users' Guide. 2nd edition, SIAM, 1995.

24. M. Berry. Large-scale sparse singular value computations. The Intern. J. of Supercomp. Appl. 6(1): 13-49 (1992).

25. G. Beylkin, R. Coifman, V. Rokhlin. Fast Wavelet Transforms and Numerical Algorithms I, Yale University Preprint, 1990.

26. A. Brandt. Multilevel computations of integral transforms and particle interactions with oscillatory kernels. Computer Physics Communications 65: 24-38 (1991).88 —

27. A. Brandt, A. Lubrecht. Multilevel matrix multiplication and fast solution of integral equations. J. Comput. Phys. 90: 348-370 (1990).

28. S. Chandrasekaran, I. Ipsen. On Rank-Revealing QR Factorizations. Research Report YALEU/DCS/RR-880, December 1991.

29. C. Eckart, G. Young. A principal axis transformation for non-hermitian matrices. Bulletin of American Mathematical Society 45: 118-121, 1939.

30. G. Golub, F. Luk, M. Overton. A block Lanczos method for computing the singular values and corresponding singular vectors of a matrix. ACM Transactions on Math. Soft. 7(2):149-169 (1981).

31. S. A. Goreinov, E. E. Tyrtyshnikov, N. L. Zamarashkin. A Theory of Pseudoskeleton Approximations. Linear Algebra and Its Applications 261: 1-21, 1997.

32. S. A. Goreinov, E. E. Tyrtyshnikov. Maximal-Volume concept in approximation by low-rank matrices. Accepted for publication in Contemporary Mathematics.

33. S. A. Goreinov, E. E. Tyrtyshnikov, A. Yu. Yeremin. Matrix-free iterative solution strategies for large dense linear systems. Numerical Linear Algebra with Applications 4(4): 273-294 (1997).

34. R. L. Graham, D. E. Knuth, O. Patashnik. Concrete Mathematics. 2nd edition, Addison-Wesley, 1994.

35. L. Greengard, V. Rokhlin. The rapid evaluation of potential fields in three dimensions. Lecture Notes in Mathematics. 1360: 121-141 (1988).

36. M. R. Hackbusch, Z. P. Nowak. On the fast matrix multiplication in the boundary elements method by panel clustering. Numerische Mathematik 54(4): 463-491 (1989).

37. A. Harten, I. Yad-Shalom. Fast multiresolution algorithms for matrix-vector multiplication. SI AM J. Numer. Anal. 31(4): 1191-1218 (1994).

38. W. K. Hayman. Power series expansions for harmonic functions. The Bulletin of the London Mathematical Society 2: 152-158 (1970).89 —

39. W. Hoffman. Iterative algorithms for Gram-Schmidt orthogonalization. Computing. 41:335-348 (1989).

40. Y. P. Hong, C.-T. Pan. Rank-revealing QR factorizations and singular value decomposition. Mathematics of Computation 58(197): 213-232, 1992.

41. C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Proc. 2nd Symposium Large-Scale Digital Calculating Machines, pp. 164-206 (1951).

42. V. I. Lebedev. Zolotarev polynomials and extremum problem. Russian Journal on Numerical Analysis and Math. Modelling. 9(3):191—314 (1994).

43. B. Maskew. Program VSAERO theory document: A computer program for calculating nonlinear aerodynamic characteristics of arbitrary configurations. NASA Contractor Report CR-4023 (AMI-8416). 1987.

44. L. Mirsky. Symmetric gauge functions and unitarily invariant norms. The Quarterly Journal of Mathematics, Oxford II series 11: 50-59, 1960.

45. S. V. Myagchilov, E. E. Tyrtyshnikov. A fast matrix-vector multiplier in discrete vortex method. Russian Journal on Numerical Analysis and Math. Modelling. 7(4):325-342 (1992).

46. C. C. Paige. The computation of eigenvalues and eigenvectors of very large sparse matrices. Ph. D. thesis. The University of London. 1971.

47. T. von Petersdorff, C. Schwab, R. Schneider. Multiwavelets for second kind integral equations. Preprint of University of Maryland. 1991.

48. S. M. Rao, D. R. Wilton, A. W. Glisson. Electromagnetic scattering by surfaces of arbitrary shape. IEEE Transactions on Antennas and Propagation AP-30: 409-418, 1982.

49. S. Rjasanow, E. E. Tyrtyshnikov. The ATS algorithm. Manuscript, 1998.

50. V. Rokhlin. Rapid solution of integral equations of classical potential theory. J. Comput. Physics. 60: 187-207 (1985).

51. V. Rokhlin. A fast algorithm for particle simulations. J. Comput. Physics. 73: 325-348 (1987).

52. V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions. J. Comput. Phys. 86: 414-439 (1990).

53. I. Saad. Numerical methods for large eigenvalue problems. Manchester University Press, 1992.

54. I. Saad. Iterative solution of indefinite symmetric linear systems by method using orthogonal polynomials over two disjoint intervals. SIAM Journal on Numerical Analysis. 20(4):784-814 (1983).

55. I. Saad, M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sei. Stat. Comput. 7:856-869 (1986).

56. I. Saad. ILUT: A dual threshold incomplete LU factorization. Numer. Linear Algebra Appl. l(4):387-402 (1994).

57. E. Schmidt. Zur Theorie der linearen und nichtlinearen Integralgleichungen. I Teil. Entwicklung willkürlichen Punktionen nach System vorgeschriebener. Mathematische Annalen 63: 433-476, 1907.91 —

58. G. W. Stewart, J. Sun. Matrix Perturbation Theory. Academic Press, 1990.

59. E. E. Tyrtyshnikov. A brief introduction to numerical analysis. Birkhauser, 1997.

60. E. E. Tyrtyshnikov. Mosaic-skeleton approximations. Calcolo 33(1-2): 47-57 (1996).

61. E. E. Tyrtyshnikov. Incomplete cross approximation in the mosaic-skeleton method. Accepted for publication in Computing.

62. R. R. Underwood. An interative block Lanczos method for the solution of large sparse symmetric eigenproblems. Ph. D. thesis. Stanford University. 1975.

63. G. N. Watson. A treatise on the theory of Bessel functions. 2nd edition, Cambridge University Press, 1992.