Подпрограмма:  aes1r_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 

Подпрограмма aes1r_c использует следующее свойство метода Ланцоша: в спектре трехдиагональной матрицы Т_j , полученной при выполнении достаточно большого числа j (j ~ 2N,...,6N) шагов алгоритма (1), появляются приближения ко всем различным собственным значениям исходной матрицы А.

В подпрограмме aes1r_c на шагах с номерами

       p = mi = M1 + (i - 1)*MH ,   i = 1, 2, ...,   (М1, МН - заданные числа) 

процесс (1) прерывается и вычисляются собственные значения матриц Т_p. Для разделения собственных значений матриц Т_p на "истинные" и "лишние" применяется критерий, использующий собственные значения матрицы Т_p^{^} , получаемой из Т_p вычеркиванием первой строки и первого столбца. Для истинных собственных значений матриц Т_p вычисляются оценки точности.

Пусть К_j - число истинных собственных значений матрицы Т_j , К'_j - число истинных значений Т_j , являющихся приближениями к собственным значениям матрицы А с точностью max { ЕРS, 10*j*ε }*|| А ||₂, где ЕРS - заданная пользователем величина, ε - относительная машинная точность, || А ||₂ - спектральная норма матрицы А.

Выход из подпрограммы aes1r_c происходит

1) при выполнении одного из соотношений

        K'p = Kp = Kp1 ,   K'p = N ;   ( где p1 = mi-1 ) 

2) по заданному ограничению на число шагов процесса (1);

3) если очередное вычисленное значение β_j+1 очень мало;

4) если при работе подпрограммы aee2r_c, используемой для вычисления собственных значений трехдиагональных матриц, произошла фатальная ошибка.

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

Сullum J., Donath W.Е. Lanczos and thе computation in specified intervals of thе spectrum of large, sparse real symmetric matrices. Sparse Мatrix Рroc., Кnoxville, Тen., 1978. Рhiladelphia, Pа, 1979, 220-255.

Х.Д. Икрамов. Разреженные матрицы. Математический анализ. Т. 20. (Итоги науки и техн. ВИНИТИ АН СССР). М., 1982.

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

    int aes1r_c (integer *n, integer *m0, integer *m1, integer *mh,
                 integer *m2, real *eps, S_fp ff, real *vj, real *w, real *e, real *
                 d, integer *l, real *ev, real *er, integer *ir, integer *ip)

   Параметры 

   Версии :   нет 

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

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

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

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

        
        int i__1;

        
    static float teta[28];
    extern int aes1r_c (int *, int *, int *, int *, int *, float *, U_fp,
                       float *, float *, float *, float *, int *, float *,
                       float *, int *, int *);
    static float a[28], b[57];
    static int i__, k;
    static float w[6];
    extern int f0_c ();
    static int m0, ip, ir[28];
    static float rr[28], eps0;

    eps0 = 1e-6f;
    m0 = 0;
    aes1r_c (&c__6, &m0, &c__4, &c__4, &c__28, &eps0, (U_fp)f0_c, vj, w, b, a,
            &k, teta, rr, ir, &ip);

    printf("\n %5i %5i %5i \n", k, ip, m0);
    for (i__ = 0; i__ <= 3; i__+=3) {
        printf("\n %15.7e %15.7e %15.7e \n",
                teta[i__], teta[i__+1], teta[i__+2]);
    }
    printf("\n %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e \n",
                rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
    printf("\n %5i %5i %5i %5i %5i %5i \n",
                ir[0], ir[1], ir[2], ir[3], ir[4], ir[5]);
    return 0;
}     

    int f0_c (float *v, float *w, int *n, float *s)
{
        
        int i__1;

        
    static int i__;

        
    --w;
    --v;

        
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        w[i__] = (i__ - *s) * v[i__] - w[i__];
    
    }
    return 0;
}     

Результаты:

      k  =  6,  ip  =  0
      teta  =   (1., 2., 3., 4., 5., 6.)
      rr  =   (1., 1., 1., 1., 1., 1.)
      ir  =   (-12, -12, -12, -12, -12, -12)