Текст подпрограммы и версий
aeb5r_c.zip  aeb5d_c.zip 
Тексты тестовых примеров
taeb5r_c.zip  taeb5d_c.zip 

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

Назначение

Вычисление методом бисекции по заданным номерам группы подряд идущих собственных значений вещественной симметричной ленточной матрицы, заданной в компактной форме.

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

Пусть все собственные значения заданной симметричной ленточной матрицы порядка N занумерованы в порядке неубывания

              λ1 ≤ λ2 ≤ ... ≤ λN-1 ≤ λN . 

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

Подпрограмма aeb5r_c позволяет вычислять различные собственные значения с различной точностью. Требуемая точность задается пользователем в специальном векторе.

Уилкинсон, Райнш. "Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра". М.: "Машиностроение", 1976.

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

    int aeb5r_c (integer *nm, integer *n, integer *m1, real *b,
            integer *l1, integer *l2, real *ev, real *eps, real *c)

Параметры

nm - число строк двумерного массива b, указанное при описании этого массива в вызывающей подпрограмме (тип: целый);
n - порядок исходной матрицы, n ≤ nm (тип: целый);
m - число нижних ненулевых кодиагоналей исходной матрицы, включая главную диагональ (тип: целый);
b - вещественный двумерный массив размерности nm * m, в первых n строках которого задана исходная симметричная ленточная матрица B компактной форме (см. Организация Библиотеки. Способы представления матриц специального вида);
l1, l2 - заданные минимальный и максимальный номера искомых собственных значений (тип: целый); возможны оба варианта: l1 ≤ l2 и l1 ≥ l2 (см. описание параметра ev);
ev - вещественный вектор длины | l2 - l1 | + 1, содержащий на выходе из подпрограммы вычисленные собственные значения; при этом, если l1 < l2, то собственные значения располагаются в порядке неубывания, а если l2 < l1, то - в порядке невозрастания;
eps - вещественный вектор длины | l2 - l1 | + 1, в компонентах которого пользователем задаются значения допустимых погрешностей для соответствующих вычисляемых собственных значений;
c - вещественный двумерный рабочий массив размерности m * min (2m - 1, n).

Версии

aeb5d_c - вычисление методом бисекции по заданным номерам группы подряд идущих собственных значений симметричной ленточной матрицы, заданной в компактной форме и с двойной точностью.

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

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

  1. 

В подпpогpамме aeb5d_c паpаметpы b, ev, eps, c имеют тип double;

  2. 

Подпpогpамма aeb5r_c (aeb5d_c) сохpаняет исходный массив b;

  3. 

При работе подпрограммы слишком малые компоненты вектора eps будут увеличены;

  4. 

Так как подпрограмма aeb5r_c применяет метод бисекции непосредственно к ленточной матрице, то ее рекомендуется использовать только, если ширина ленты невелика и число искомых собственных значений достаточно мало, точнее, если выполнно неравенство

                                            l
           n  >  0.25 (m - 1)    ∑    r i ,
                                          i =1 
где l = | l2 - l1 | + 1,  ri = - log2 (eps ( i ) / || A || ),  || A || - спектральная норма исходной матрицы A. В противном случае выгоднее исходную матрицу привести к трехдиагональному виду и применять метод бисекции к трехдиагональной матрице.

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

int main(void)
{
    /* Initialized data */
    static float b[15] /* was [5][3] */ = { 0.f,0.f,4.f,0.f,4.f,0.f,3.f,3.f,
                                            3.f,3.f,5.f,3.25f,1.f,4.25f,6.f };
    static float eps[2] = { 1e-5f,1e-5f };

    /* Local variables */
    extern int aeb5r_c(int *, int *, int *, float *, int *, int *,
                       float *, float *, float *);
    static float c__[15] /* was [3][5] */;
    static int m, n, l1, l2;
    static float ev[2];
    static int nm;

    n = 5;
    nm = 5;
    m = 3;
    l1 = 2;
    l2 = 3;
    aeb5r_c(&nm, &n, &m, b, &l1, &l2, ev, eps, c__);

    printf("\n %15.8e %15.8e \n", ev[0], ev[1]);
    printf("\n %15.2e %15.2e \n", eps[0], eps[1]);
    return 0;
} /* main */


Результаты:

               ev  =  ( .10000000e+01, .20000000e+01 )

              eps  =  ( .10e-04, .10e-04 )