|
Текст подпрограммы и версий ( Фортран ) ags0r.zip ags0d.zip |
Тексты тестовых примеров ( Фортран ) tags0r.zip tags0d.zip |
|
Текст подпрограммы и версий ( Си ) ags0r_c.zip ags0d_c.zip |
Тексты тестовых примеров ( Си ) tags0r_c.zip tags0d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) ags0r_p.zip ags0e_p.zip |
Тексты тестовых примеров ( Паскаль ) tags0r_p.zip tags0e_p.zip |
Вычисление методом Ланцоша с выборочной ортогонализацией нескольких крайних собственных значений и соответствующих собственных векторов обобщенной проблемы собственных значений Ax = λ Bx с вещественными симметричными матрицами A и B и положительно определенной матрицей B.
Метод Ланцоша для обобщенной проблемы собственных значений
(1) A x = λ B x ,
где A и B - заданные вещественные симметричные матрицы порядка N, B > 0, и заданного начального вектора v ∈ EN, v ≠ 0, состоит в построении последовательности векторов Ланцоша v1, v2, ..., vj и последовательности симметричных трехдиагональных матриц T1, T2, ..., Tj по следующим формулам:
v1 = v / (vT Bv)1/2 , β1 = 0 ,
ui = Avi - βi Bvi - 1 ,
(2) αi = uiT vi ,
T1 = α1 , i = 1 ,
| | 0 |
| | . |
Ti = | Ti-1 | 0 | , i > 1 ,
| _______| βi |
| 0 ... 0 βi αi |
wi + 1 = ui - αi Bvi ,
решение системы Bui +1 = wi + 1 ,
βi + 1 = (uTi + 1 wi + 1)1/2 ,
vi + 1 = ui + 1 / βi + 1 i = 1, ..., j .
Метод Ланцоша (2) может быть рассмотрен как реализация для задачи (1) метода Ритца, где в качестве подпространства аппроксимации используется подпространство Крылова K j = span { v1, B- 1 Av1, ..., (B- 1 A)j - 1 v1 }, в котором задано скалярное произведение (x, y)B = xT By. При этом, если (θi j, si j), 1 ≤ i ≤ j, - собственная пара симметричной трехдиагональной матрицы Tj, получаемой на j - ом шаге метода Ланцоша (2), а Vj - матрица, столбцами которой являются векторы Ланцоша v1, v2, ..., vj, то (θi j, Vj si j) будет соответствующей парой Ритца для пучка (A, B).
Полезной особенностью метода Ланцоша является то, что часто уже при j ≈ 2 √N крайние пары Ритца оказываются хорошими приближениями к искомым крайним собственным парам пучка (A, B).
В подпрограмме AGS0R используется метод Ланцоша с выборочной ортогонализацией - устойчивая к влиянию ошибок округления модификамция простого метода Ланцоша (2), отличающаяся от него тем, что на некторых шагах процесса (2) перед вычислением vi + 1 проводится B - ортоганализация вектора ui + 1 относительно некоторых векторов Ритца.
Как видно из формул (2), метод Ланцоша использует исходную матрицу A только в операциях умножения матрицы на вектор, а матрицу B - в операциях умножения матрицы на вектор и решения системы с матрицей B, поэтому всю необходимую информацию о матрицах A и B удобно задавать в виде трех соответствующих подпрограмм, в которых пользователь может максимально учесть специфику матриц. Особенно удобен метод Ланцоша для разреженных или имеющих регулярную структуру матриц.
В подпрограмме AGS0R гарантированной точностью по невязке для вычисленных собственных пар является величина k √ε || B- 1 A ||, где k ~ 1, ε - относительная машинная точность, но часто возможно получение лучшей точности. Под невязкой для собственной пары (θ, y) здесь понимается величина || B- 1 Ay - θy ||B / || y ||B.
Выход из подпрограммы происходит, когда все требуемые собственные векторы будут вычислены с заданной точностью по невязке или если продолжение процесса (2) окажется невозможным или бесполезным.
Б.Парлетт. Симметричная проблема собственных значений. Численные методы. М.: Мир, 1983.
SUBROUTINE AGS0R (N, KL, KR, FA, FB, SB, EPS, M, VJ, K, L,
R, R1, R2, EV, ER, V, Y, IP)
Параметры
| N - | порядок исходной матрицы (тип: целый); |
| KL, KR - | задаваемые соответственно число наименьших и число наибольших собственных значений и соответствующих собственных векторов, которые требуется вычислить (тип: целый); |
| FA - |
имя подпрограммы вычисления Z: = AX - Z; первый
оператор подпрограммы должен иметь вид: SUBROUTINE FA (X, Z, N) Здесь: |
| X - | заданный вещественный вектор длины N; |
| Z - | вещественный вектор длины N, содержащий на входе в подпрограмму заданный вектор, а на выходе из подпрограммы - вычисленный вектор AX - Z; |
| N - | порядок матрицы A (тип: целый); |
| FB - |
имя подпрограммы умножения матрицы B на вектор;
первый оператор подпрограммы должен иметь вид: SUBROUTINE FB (X, Z, N) Здесь: |
| X - | заданный вещественный вектор длины N; |
| Z - | вещественный вектор длины N, содержащий на выходе из подпрограммы произведение B X; |
| N - | порядок матрицы B (тип: целый); |
| SB - |
имя подпрограммы решения системы с матрицей B;
первый оператор подпрограммы должен иметь вид: SUBROUTINE SB (X, Z, N) Здесь: |
| X - | вещественный вектор длины N, содержащий правую часть системы BZ = X; |
| Z - | вещественный вектор длины N, содержащий на выходе из подпрограммы SB вычисленное решение системы BZ = X; |
| N - | порядок матрицы B (тип: целый) |
| EPS - | заданная требуемая точность по невязке для вычисляемых собственных значений и соответствующих собственных векторов (тип: вещественный); |
| M - | целочисленная переменная, содержащая на входе в подпрограмму задаваемое максимальное число шагов алгоритма, а на выходе из подпрограммы - число выполненных шагов; |
| VJ - | вещественный вектор длины N, содержащий на входе в подпрограмму заданный начальный вектор; если VJT * B * VJ = 0, то в качестве начального будет взят вектор (1, ..., 1)T; |
| K - | заданное максимальное число вычисляемых собственных значений и соответствующих собственных векторов матриц Tj, K ≥ KL + KR + 2 (тип: целый); |
| L - | длина рабочего вектора R, L ≥ M * (K + 9) + 8 * K (тип: целый); |
| R - | вещественный вектор длины L, используемый как рабочий; |
| R1, R2 - | вещественные векторы длины N, используемые как рабочие; |
| EV - | вещественный вектор длины (KL + KR), в первых KL компонентах которого на выходе из подпрограммы содержатся в порядке неубывания вычисленные KL наименьших собственных значений, а в последних KR компонентах в порядке невозрастания - KR наибольших собственных значений; |
| ER - | вещественный вектор длины (KL + KR), содержащий невязки для вычисленных собственных пар; |
| V - | вещественный двумерный массив размерности N * M, в столбцах которого в ходе работы подпрограммы располагаются вычисленные векторы Ланцоша; |
| Y - | вещественный двумерный массив размерности N * K, в I - ом, 1 ≤ I ≤ KL + KR, столбце которого на выходе из подпрограммы содержится вычисленный B - нормированный собственный вектор, соответствующий выходному собственному значению EV ( I ); |
| IP - | простая целочисленная переменная, содержащая на входе в подпрограмму признак решаемой задачи, при этом: |
| IP= 1 - | если требуется выполнить ровно M шагов алгоритма; |
| IP= 0 - | если допустим выход из подпрограммы при j < M, j - число выполненных шагов; в этом случае вычисления будут прекращены после того, как на некотором шаге j будет получено заданное число приближений к собственным парам задачи (1) с заданной точностью по невязке; |
| на выходе из подпрограммы IP содержит признак выхода, причем: |
| IP= -1 - | если заданное значение параметра L не удовлетворяет неравенству L ≥ M (9 + K) + 8K; |
| IP= -2 - | если заданное значение параметра K не удовлетворяет неравенству K ≥ KL + KR + 2; |
| IP= -3 - | если обе описанные выше ошибки имели место; |
| IP= 0 - | если требуемые собственные значения и соответствующие собственные векторы вычислены с заданной точностью; |
| IP= 1 - | если выполнено заданное число M шагов алгоритма, а заданная точностью еще не достигнута; |
| IP= 2 - | если на некотором шаге j, j < M, алгоритма потребовалось вычислить более K собственных пар матрицы Tj; в этом случае для продолжения процесса не хватает памяти; |
| IP= 3 - | если заданная точность не может быть достигнута; |
| IP= 4 - | если получено слишком малое значение βj + 1; |
| IP= 5 - | если резко возрастает объем работ по проведению выборочной ортогонализации |
Версии
| AGS0D - | вычисление методом Ланцоша с выборочной ортогонализацией нескольких крайних собственных пар обобщенной проблемы собственных значений Ax = λ Bx, где A и B - вещественные симметричные матрицы, B > 0, в случае, когда информация о матрицах A и B задана с двойной точностью |
Вызываемые подпрограммы
|
AVZ6R - AVZ6D | упорядочивание вектора по возрастанию его компонент с запоминанием произведенных перестановок; |
|
AVZ2R - AVZ2D | нахождение максимальной по модулю компоненты и ее индекса из всего вектора или из заданного подмножества компонент этого вектора; |
|
AM05R - AM05D | вычисление произведения матрицы на вектор |
Замечания по использованию
| 1. |
Наименование подпрограмм, соответствующих формальным параметрам FA, FB, SB, должны быть указаны в операторе EXTERNAL в программе пользователя, вызывающей AGS0R; | |
| 2. |
Подпрограммы FA, FB, SB должны сохранять вектор X; | |
| 3. |
Если начальный вектор был выбран неудачно, то возможно, что вычисленные собственные пары не обязательно будут соответствовать первым KL и последним KR собственным парам задачи (1). Поэтому, вообще говоря, для проверки необходимо обратиться к AGS0R с разными начальными векторами. Особенно это необходимо, если на выходе из подпрограммы IP = 4 или мало выходное значение параметра M; | |
| 4. |
Неудачный выбор начального вектора может привести к выходу из подпрограммы на шаге с номером j < KL + KR по признаку IP = 4. В этом случае будут вычислены все j собственных пар, причем все j собственных значений будут расположены в массиве EV в порядке неубывания. Понятно, что в этом случае повторное обращение к подпрограмме совершенно необходимо; | |
| 5. |
Если на выходе из подпрограммы IP < 0, то никакой информации об искомых собственных парах получено не будет; если же IP > 0, то на выходе из подпрограммы будут получены KL первых и KR последних пар Ритца, вычисленных на шаге прерывания (номер этого шага - выходное значение M), т.е. будут получены приближения к искомым собственным парам, но точность их еще не является удовлетворительной. Эти приближения могут быть использованы для построения нового начального вектора и повторного обращения к AGS0R; | |
| 6. |
Если на выходе из подпрограммы IP = 2, то необходимо увеличить значение параметра K и соответственно L. Рекомендации для выбора значения K, которые были бы удовлетворительными для любых матриц, не могут быть даны, т.к. требуемое значение K зависит от спектра исходной матрицы и равно числу пар Ритца, невязки которых станут меньше √ε || B- 1 A || к тому моменту, когда невязка для заданного числа крайних пар Ритца станут меньше EPS; | |
| 7. |
В подпрограмме AGS0D параметры EPS, VJ, R, R1, R2, EV, ER, V, Y имеют тип DOUBLE PRECISION; | |
| 8. |
Подпрограммы AGS0R1(AGS0D1), AES0R0(AES0D0), AV04R1(AV04D1) используются в качестве рабочих; | |
| 9. |
Для подпрограммы AGS0D в подпрограммах FA, FB и SB параметры X и Z, должны иметь тип DOUBLE PRECISION; | |
| 10. |
Если для матрицы B легко выполнимо разложение Холецкого
B = LLT, то метод Ланцоша можно применять
непосредственно к неявно заданной матрице
C = L- 1 A (LT) - 1,
т.е. использовать не AGS0R, а подпрограмму, реализующую метод
Ланцоша для стандартной проблемы собственных значений
Пример использования SUBROUTINE FA (X, Z, N)
DOUBLE PRECISION X, Z, RAB
DIMENSION X(N) Z(N), RAB(3)
DATA RAB /1.D0, 9.D0, 14.D0/
DO 10 I = 1, N
Z(I) = RAB(I)*X(I) - Z(I)
10 CONTINUE
RETURN
END
SUBROUTINE SB (X, Z, N)
DOUBLE PRECISION X, Z, RAB
DIMENSION X(N), Z(N), RAB(3)
DATA RAB /2*1.D0, 5.D - 1/
DO 10 I = 1, N
Z(I) = X(I)*RAB(I)
10 CONTINUE
RETURN
END
SUBROUTINE FB (X, Z, N)
DOUBLE PRECISION X, Z, RAB
DIMENSION X(N), Z(N), RAB(3)
DATA RAB /2*1.D0, 2.D0/
DO 10 I = 1, N
Z(I) = X(I)*RAB(I)
10 CONTINUE
RETURN
END
DIMENSION VJ(3), EV(3), ER(3), R(82), R1(3), R2(3),
* V(3, 3), Y(3, 5)
DOUBLE PRECISION VJ, EV, ER, R, R1, R2, V, Y
EXTERNAL FA, FB, SB
DATA VJ /3*1.D0/
M = 3
IP = 0
CALL AGS0D (3, 3, 0, FA, FB, SB, 1.D - 7, M, VJ, 5, 82, R,
* R1, R2, EV, ER, V, Y, IP)
Результаты:
EV = (1., 7., 9.) ,
ER = (1.D - 15, 4.D - 15, 9.D - 16) ,
| - 1. 5.D - 17 - 9.D - 18 |
Y = | 2.D - 16 2.D - 15 1. |
| 2.D - 16 0.70710678118655 - 3.D - 16 |
IP = 0
M = 3
|