Подпрограмма:  AES0R

Назначение

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

Математическое описание

Метод Ланцоша для заданной симметричной матрицы А порядка N и начального вектора v ∈ Е^N ,  v ≠ 0 ,  строит последовательность векторов Ланцоша  v₁, v₂, ..., v_j и последовательность трехдиагональных матриц  Т₁, Т₂, ..., Т_j по следующим формулам процесса (1):

      v1 = v / (vTv)1/2  ,     β1 = 0  , 

      ui = Avi - βivi-1  ,     αi = uiTvi  , 

  T1 =  α1 ,                            i = 1 ,

          |               |  0 |
          |               |  .  |
  Ti =  |    Ti-1     |  0 |  ,          i > 1 ,
          | _______| βi |
          | 0 ... 0 βi   αi |

   wi+1 = ui - αivi ,    βi+1 = (wTi+1wi+1)1/2 ,    vi+1 = wi+1 / βi+1  ,
    i = 1, ..., j 

Метод Ланцоша может быть рассмотрен как реализация для спектральной задачи Az = λz метода проекций Ритца - Галеркина на последовательность подпространств Крылова

     K( j ) = span { v1 , Av1, ... , Aj-1v1 }  ,   j = 1, ..., N  , 

при этом числа Ритца (собственные значения сужения оператора Р_jА на К ^( j ), где Р_j - ортогональный проектор на К ^( j ) совпадают с собственными значениями Θ_i ^( j ) , i = 1, ..., j матрицы Т_j , получаемой на j - ом шаге метода Ланцоша, а векторы Ритца равны V_js_i^( j ), где s_i^( j ) - соответствующие Θ_i ^( j ) собственные векторы Т_j , а V_j - матрица, столбцами которой являются векторы Ланцоша  v₁, v₂, ..., v_j .

Полезной особенностью метода Ланцоша является то, что часто уже при j ~ 2√N  крайние пары Ритца (Θ_i ^( j ), V_js_i ^( j )) оказываются хорошими приближениями к искомым крайним собственным парам матрицы А.

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

В подпрограмме АЕS0R используется метод Ланцоша с выборочной ортогонализацией - устойчивая к ошибкам округления модификация простого метода Ланцоша, отличающаяся от него тем, что на некоторых шагах процесса (1) проводится ортогонализация w_j+1 относительно некоторых векторов Ритца.

В подпрограмме АЕS0R метод Ланцоша используется как итерационный: после выполнения заданного числа М шагов процесс (1) прерывается (заканчивается один крупный шаг) и, если вычислены еще не все требуемые векторы, то строится новый начальный вектор (линейная комбинация крайних векторов Ритца, точность которых еще не является удовлетворительной) и процесс (1) начинается заново (начинается следующий крупный шаг).

В подпрограмме АЕS0R предполагается, что спектральная норма || А ||₂ ~ 1 . Гарантированной точностью по невязке для вычисляемых собственных пар является величина k√ε || А ||₂ , где ε - относительная машинная точность, k ~ 1 , но часто возможно получение лучшей точности.

Выход из подпрограммы происходит, когда все требуемые собственные пары будут вычислены с точностью ЕРS по невязке, где ЕРS задаваемая величина, или по заданному максимально допустимому числу NI крупных шагов.

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

Использование

    SUBROUTINE  AES0R ( N, KL, KR, KL0, KR0, M, NI, VJ, F, EPS,
                                             L, R, EV, Y, ER, IP) 

   Параметры 

   Версии :   нет 

   Вызываемые подпрограммы 

   Замечания по использованию  

Пример использования

F0 - подпрограмма пользователя, вычисляющая W:  = (A - S*I)V - W
       для A = diag (0., 1.E-5, 2.E-5, 3.E-5, 0.1, 1.)

       SUBROUTINE  F0 (X, Z, N, S)
       DIMENSION  X(N), Z(N), RAB(6)
       DATA RAB /0., 1.E-5, 2.E-5, 3.E-5, 0.1, 1./
       DO 10  I = 1, N
       Z(I) = (RAB(I) - S)*X(I) - Z(I)
  10 CONTINUE
       RETURN
       END

       DIMENSION  Y(6, 2), VJ(6), EV(2), ER(2), RAB(211)
       DATA VJ /6*1./
       EXTERNAL F0
       KL0 = 0
       KR0 = 0
       IP = 0
       NI = 1
       CALL  AES0R (6, 0, 2, KL0, KR0, 5, NI, VJ, F0, 1.E-6, 211, RAB,
      *                          EV, Y, ER, IP)

Результаты:

       EV = (1., 0.1) ,   ER = (2.2E-11, 2.2E-9) ,   IP = 0 ,   NI = 1 ,
  
       Y1  =  (-1.1E-11,  1.1E-11,  1.1E-11,  -1.2E-11,  -2.0E-13,  1.)T
       Y2  =  (-1.1E-8,  1.1E-8,  1.1E-8,  -1.1E-8,  -1.,  1.2E-10) .