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

Назначение

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

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

Метод Ланцоша для заданной симметричной матрицы А порядка 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 )) оказываются хорошими приближениями к искомым крайним собственным парам матрицы А.

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

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

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

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

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

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

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

    int aes0r_c (integer *n, integer *kl, integer *kr, integer *f, real *
                 eps, integer *l, real *r, real *ev, real *y, real *er, integer *ip)

   Параметры 

   Версии :   нет 

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

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

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

f0_c - подпрограмма пользователя, вычисляющая w:  = (A - s*I)v - w
       для a = diag (0., 1.e-5, 2.e-5, 3.e-5, 0.1, 1.)

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 aes0r_c(int *, int *, int *, int *, int *, int *, int *,
                       float *, U_fp, float *, int *, float *, float *,
                       float *, float *, int *);
    static int i__;
    static float y[12] /* was [6][2] */;
    extern int f0_c();
    static float er[2];
    static int ni, ip;
    static float ev[2];
    static int kl0, kr0;
    static float rab[211];

    kl0 = 0;
    kr0 = 0;
    ip = 0;
    ni = 1;
    aes0r_c(&c__6, &c__0, &c__2, &kl0, &kr0, &c__5, &ni, vj, (U_fp)f0_c, &c_b5,
            &c__211, rab, ev, y, er, &ip);

    printf("\n %16.7e %16.7e \n", ev[0], ev[1]);
    printf("\n %13.2e %13.2e \n", er[0], er[1]);
    printf("\n %5i %5i \n", ip, ni);
    for (i__ = 0; i__ <= 9; i__+=3) {
        printf("\n %10.2e %10.2e %10.2e \n", y[i__], y[i__+1], y[i__+2]);
    }
    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 */


Результаты:

       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) .