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

Назначение

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

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

Метод Ланцоша для заданной симметричной матрицы А порядка N и начального вектора  v ∈ E^N,  v ≠ 0, строит последовательность векторов  v₁, v₂,..., v_j и последовательность трехдиагональных матриц T₁, T₂,..., T_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 -1 v1 }   ,     j = 1,... , N   , 

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

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

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

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

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

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

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

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

    int aes4r_c (integer *n, integer *kl, integer *kr, real *ff,
                 real *eps, integer *m, real *vj, integer *k, integer *l, real *r, 
                 real *r1, real *r2, real *ev, real *er, real *v, real *y, integer *ip)

   Параметры 

   Версии 

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

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

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

int main(void)
{
    /* Initialized data */
    static float vj[6] = { 1.f,1.f,1.f,1.f,1.f,1.f };

    /* Local variables */
    extern int aes4r_c(int *, int *, int *, U_fp, float *, int *, float *,
                       int *, int *, float *, float *, float *, float *,
                       float *, float *, float *, int *);
    static int i__, m;
    static float r__[97], v[30] /* was [6][5] */, y[12] /* was [6][2] */;
    extern int f0_c();
    static float r1[6], r2[6], er[2];
    static int ip;
    static float ev[2];

#define y_ref(a_1,a_2) y[(a_2)*6 + a_1 - 7]

    m = 5;
    ip = 0;
    aes4r_c(&c__6, &c__0, &c__2, (U_fp)f0_c, &c_b4, &m, vj, &c__4, &c__97, r__,
            r1, r2, ev, er, v, y, &ip);

    printf("\n %5i %5i \n", ip, m);
    if (ip < 0) {
        goto l55;
    }
    printf("\n %17.7e %17.7e \n", ev[0], ev[1]);
    for (i__ = 1; i__ <= 6; ++i__) {
        printf("\n %17.7e %17.7e \n", y_ref(i__,1), y_ref(i__,2));
    }
    printf("\n %17.7e %17.7e \n", er[0], er[1]);
l55:
    return 0;
} /* main */


int f0_c(float *x, float *z__, int *n, float *s)
{
    /* Initialized data */
    static float rab[6] = { 0.f,1e-5f,2e-5f,3e-5f,.1f,1.f };

    /* System generated locals */
    int i__1;

    /* Local variables */
    static int i__;

    /* Parameter adjustments */
    --z__;
    --x;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        z__[i__] = (rab[i__ - 1] - *s) * x[i__] - z__[i__];
/* l10: */
    }
    return 0;
} /* f0_c */


Результаты:

         ip  =  0   ,     m  =  3
        
         coбcтвeнныe значения:

         ev(1)  =  1.   ,     ev(2)  =  0.1
         
         coбcтвeнныe векторы:  

                   |   5.e - 7           2.e - 4  |
                   | - 6.e - 7           5.e - 5  |
                   | - 2.e - 6         - 6.e - 5  |
      y_ref  =   | - 3.e - 6         - 2.e - 4  |
                   |   1.e - 6         - 1.          |
                   |   1.                   3.e - 6  |
                     
         нeвязkи:     er(1)  =  4.e - 6   ,     er(2)  =  2.e - 5