Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/ag_htm_c/agh6r_c.htm
Дата изменения: Mon Dec 8 15:07:48 2014
Дата индексирования: Sun Apr 10 01:39:23 2016
Кодировка: Windows-1251
БЧА НИВЦ МГУ. AGH6R_C. Обобщенная проблема собственных значений для симметричных и эрмитовых матриц
Текст подпрограммы и версий
agh6r_c.zip  agh9r_c.zip  		Тексты тестовых примеров
tagh6r_c.zip  tagh9r_c.zip 

Подпрограмма:  agh6r_c (версия agh9r_c)

Назначение

Вычисление собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов в обобщенной проблеме ABx = λx (или BAx = λx) для вещественных симметрических матриц A и B.

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

Подпрограмма agh6r_c реализует алгоритм вычисления собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов уравнения вида АВx = λx, где A и B - вещественные симметрические матрицы, и матрица B положительно определена.

При помощи разложения Холецкого для матрицы B:  В = LLT уравнение ABx = λx приводится к стандартному виду Qy = λy,  где Q = LTAL,  y = LTx. Стандартная задача решается путем приведения симметрической матрицы Q к трехдиагональному виду Q1 и вычисления собственных значений, принадлежащих заданному интервалу, методом бисекций и соответствующих собственных вектоpов, используя метод обратных итераций.

Восстановление собственных вектоpов  x  исходного уравнения АВx = λx из соответствующих вектоpов  z  уравнения Q1z = λz осуществляется согласно соотношению: x = L - TPz, где P - матрица преобразования Q1 = PTQP, L - преобpазование Холецкого, причем xTВx = E,  где E - единичная матрица.

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

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

    int agh6r_c (integer *n, integer *mm, integer *m, real *rlb,
            real *rub, real *a, real *b, real *ev, real *v, integer *irab,
            real *rab, integer *ierr)

   Параметры
n - порядок исходных матриц (тип: целый);
mm - оценка свеpху числа собственных значений уpавнения ABx = λx, принадлежащих заданному интервалу mm ≤ n (тип: целый); если фактическое число собственных значений m, принадлежащих заданному интервалу, больше, чем mm, то собственные значения и соответствующие собственные векторы не вычисляются;
m - целая переменная, в которой запоминается вычисленное число собственных значений, принадлежащих заданному интервалу;
   rlb -
   rub   		заданные нижняя и верхняя границы интервала собственных значений (тип: вещественный); если rlb > rub, то собственные значения не вычисляются;
a, b - вещественные двумерные массивы размера n на n, содержащие исходные симметрические матрицы A и B соответственно;
ev - вещественный вектоp длины mm, содержаший вычисленные в возрастающем порядке собственные значения;
v - вещественный двумерный массив размера n на mm, содержащий в первых m столбцах вычисленные собственные векторы, нормированные так, что vTBv = E;
irab - целый вектоp длины mm, содержащий индексы расположенных в возрастающем порядке m собственных значений;
rab - вещественный вектоp длины 9 на n, используемый как рабочий;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом значение ierr:
 

- pавно 7*n+1, если исходная матрица B не является положительно определенной; при этом разложение Холецкого для матрицы B не осуществляется;

- pавно 3*n+1, если значение mm меньше истинного числа вычисленных собственных значений m на интервале; при этом собственные значения и соответствующие собственные векторы не вычисляются.

- pавно - k, если для вычисления собственного вектоpа с индексом k потребовалось более 5 итераций; при этом компоненты этого вектоpа полагаются равными нулю. Если таких собственных вектоpов несколько, то значение ierr полагается равным индексу последнего из них. 
   Версии
agh9r_c - вычисление собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов в обобщенной проблеме BAx = λx для вещественных симметричных матриц A и B. 
   Вызываемые подпрограммы
utag10_c - подпрограмма выдачи диагностических сообщений при работе подпрограмм agh6r_c и agh9r_c. 
   Замечания по использованию
 

Исходные матрицы A и B можно задавать лишь верхними треугольными половинами.

Подпрограмма agh6r_c сохраняет строгий верхний треугольник массива a и полный верхний треугольник массива b, остальные элементы массивов a и b используются как pабочие.

Программа agh9r_c нормирует собственные векторы  x  таким образом, что xTB - 1x = E.

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

int main(void)
{
    /* Initialized data */
    static float a[25]   /* was [5][5] */ = { 10.f,0.f,0.f,0.f,0.f,2.f,12.f,
            0.f,0.f,0.f,3.f,1.f,11.f,0.f,0.f,1.f,2.f,1.f,9.f,0.f,1.f,1.f,-1.f,
            1.f,15.f };
    static float b[25]   /* was [5][5] */ = { 12.f,0.f,0.f,0.f,0.f,1.f,14.f,
            0.f,0.f,0.f,-1.f,1.f,16.f,0.f,0.f,2.f,-1.f,-1.f,12.f,0.f,1.f,1.f,
            1.f,-1.f,11.f };

    /* Local variables */
    static int irab[5], ierr;
    extern int agh6r_c(int *, int *, int *, float *, float *, float *,
                       float *, float *, float *, int *, float *, int *);
    static int i__, m, n;
    static float v[25] /* was [5][5] */;
    static int mm;
    static float ev[5], rab[45], rlb, rub;

#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]
#define b_ref(a_1,a_2) b[(a_2)*5 + a_1 - 6]
#define v_ref(a_1,a_2) v[(a_2)*5 + a_1 - 6]

    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e %15.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
                a_ref(i__, 4), a_ref(i__, 5));
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e %15.7e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3),
                b_ref(i__, 4), b_ref(i__, 5));
    }
    n = 5;
    mm = 5;
    rlb = 70.f;
    rub = 300.f;
    agh6r_c(&n, &mm, &m, &rlb, &rub, a, b, ev, v, irab, rab, &ierr);

    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e %15.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
                a_ref(i__, 4), a_ref(i__, 5));
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e %15.7e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3),
                b_ref(i__, 4), b_ref(i__, 5));
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e \n", ev[i__-1]);
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e %15.7e \n",
                v_ref(i__, 1), v_ref(i__, 2), v_ref(i__, 3),
                v_ref(i__, 4), v_ref(i__, 5));
    }
    printf("\n %5i \n", ierr);
    return 0;
} /* main */


Результаты:

Собственные значения в интервале (70., 300.)

                 |  77.697191195 |
                 | 112.15419325  |
      ev  =  | 134.68646332  |
                 | 167.48487891  |
                 | 242.97727332  |

Собственные векторы, соответствующие вычисленным 
в интервале (70., 300.) собственным значениям:

                 |  0.2349114135 |                      |  0.1288556917 |
                 | -0.0410915167 |                     | -0.1193865988 |
      v1  =  | -0.0383075946 |    ,      v2  =  | -0.0282771880 |    , 
                 | -0.2059003675 |                     |  0.1923580004 |
                 | -0.0734707966 |                     | -0.0097623271 |
 
                 |  0.0042355205 |                     |  0.0183136812 |
                 | -0.1812063856 |                    | -0.0266749519 |
      v3  =  |  0.1210383986 |    ,      v4  =  |  0.1834456078 |    , 
                 | -0.0609182758 |                     |  0.0051904406 |
                 |  0.1690213925 |                     | -0.2218442867 |

                 | -0.1249195279 |
                 | -0.1535463561 |
      v5  =  | -0.1145245145 |    ,      ierr  =  0
                 | -0.0657938487 |
                 | -0.1010161054 |