Текст подпрограммы и версий
ags0r_c.zip  ags0d_c.zip  		Тексты тестовых примеров
tags0r_c.zip  tags0d_c.zip 

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

Назначение

Вычисление методом Ланцоша с выборочной ортогонализацией нескольких крайних собственных значений и соответствующих собственных векторов обобщенной проблемы собственных значений Ax = λ Bx с вещественными симметричными матрицами  A и  B и положительно определенной матрицей  B.

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

Метод Ланцоша для обобщенной проблемы собственных значений 

(1)       A x  =  λ B x ,

где  A и  B - заданные вещественные симметричные матрицы порядка  N,  B > 0, и заданного начального вектора  v ∈ EN,  v ≠ 0, состоит в построении последовательности векторов Ланцоша  v1, v2, ..., vj и последовательности симметричных трехдиагональных матриц  T1, T2, ..., Tj по следующим формулам: 

         v1  =  v / (vT Bv)1/2 ,    β1  =  0 ,
         ui   =  Avi - βi Bvi - 1 ,
(2)     αi  =  uiT vi ,

         T1  =  α1 ,                           i = 1 ,

                 |               |  0 |
                 |               |  .  |
         Ti =  |    Ti-1     |  0 |  ,          i > 1 ,
                 | _______| βi |
                 | 0 ... 0 βi   αi |

         wi + 1  =  ui - αi Bvi ,

         решение системы    Bui +1  =  wi + 1 ,
         βi + 1  =  (uTi + 1 wi + 1)1/2 ,
         vi + 1  =  ui + 1 / βi + 1             i = 1, ..., j .

Метод Ланцоша (2) может быть рассмотрен как реализация для задачи (1) метода Ритца, где в качестве подпространства аппроксимации используется подпространство  Крылова   K j = span { v1, B- 1 Av1, ..., (B- 1 A)j - 1 v1 }, в котором задано скалярное произведение (x, y)B = xT By. При этом, если (θi j, si j),  1 ≤ i ≤ j, - собственная пара симметричной трехдиагональной матрицы  Tj, получаемой на  j - ом шаге метода Ланцоша (2), а  Vj - матрица, столбцами которой являются векторы Ланцоша  v1, v2, ..., vj, то (θi j, Vj si j) будет соответствующей парой Ритца для пучка (A, B).

Полезной особенностью метода Ланцоша является то, что часто уже при  j ≈ 2 √N крайние пары Ритца оказываются хорошими приближениями к искомым крайним собственным парам пучка (A, B).

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

Как видно из формул (2), метод Ланцоша использует исходную матрицу  A только в операциях умножения матрицы на вектор, а матрицу  B - в операциях умножения матрицы на вектор и решения системы с матрицей  B, поэтому всю необходимую информацию о матрицах  A и  B удобно задавать в виде трех соответствующих подпрограмм, в которых пользователь может максимально учесть специфику матриц. Особенно удобен метод Ланцоша для разреженных или имеющих регулярную структуру матриц.

В подпрограмме ags0r_c гарантированной точностью по невязке для вычисленных собственных пар является величина  k √ε  || B- 1 A ||, где  k ~ 1,  ε - относительная машинная точность, но часто возможно получение лучшей точности. Под невязкой для собственной пары (θ, y) здесь понимается величина  || B- 1 Ay - θy ||B / || y ||B.

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

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

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

    int ags0r_c (integer *n, integer *kl, integer *kr, U_fp fa,
                 U_fp fb, U_fp sb, 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)

   Параметры
n - порядок исходной матрицы (тип: целый);
kl, kr - задаваемые соответственно число наименьших и число наибольших собственных значений и соответствующих собственных векторов, которые требуется вычислить (тип: целый);
fa - имя подпрограммы вычисления z: = Ax - z; первый оператор подпрограммы должен иметь вид:
int fa (float *x, float *z, int *n)
Здесь:
x - заданный вещественный вектор длины  n;
z - вещественный вектор длины  n, содержащий на входе в подпрограмму заданный вектор, а на выходе из подпрограммы - вычисленный вектор Ax - z;
n - порядок матрицы  A (тип: целый);
fb - имя подпрограммы умножения матрицы B на вектор; первый оператор подпрограммы должен иметь вид:
int fb (float *x, float *z, int *n)
Здесь:
x - заданный вещественный вектор длины  n;
z - вещественный вектор длины  n, содержащий на выходе из подпрограммы произведение B x;
n - порядок матрицы  B (тип: целый);
sb - имя подпрограммы решения системы с матрицей  B; первый оператор подпрограммы должен иметь вид:
int sb (float *x, float *z, int *n)
Здесь:
x - вещественный вектор длины  n, содержащий правую часть системы Bz = x;
z - вещественный вектор длины  n, содержащий на выходе из подпрограммы sb вычисленное решение системы Bz = x;
n - порядок матрицы  B (тип: целый)
eps - заданная требуемая точность по невязке для вычисляемых собственных значений и соответствующих собственных векторов (тип: вещественный);
m - целочисленная переменная, содержащая на входе в подпрограмму задаваемое максимальное число шагов алгоритма, а на выходе из подпрограммы - число выполненных шагов;
vj - вещественный вектор длины  n, содержащий на входе в подпрограмму заданный начальный вектор; если  vjT * B * vj = 0, то в качестве начального будет взят вектор (1, ..., 1)T;
k - заданное максимальное число вычисляемых собственных значений и соответствующих собственных векторов матриц  Tj,  k ≥ kl + kr + 2 (тип: целый);
l - длина рабочего вектора  r,  l ≥ m * (k + 9) + 8 * k (тип: целый);
r - вещественный вектор длины  l, используемый как рабочий;
r1, r2 - вещественные векторы длины  n, используемые как рабочие;
ev - вещественный вектор длины (kl + kr), в первых kl компонентах которого на выходе из подпрограммы содержатся в порядке неубывания вычисленные kl наименьших собственных значений, а в последних kr компонентах в порядке невозрастания - kr наибольших собственных значений;
er - вещественный вектор длины (kl + kr), содержащий невязки для вычисленных собственных пар;
v - вещественный двумерный массив размерности n * m, в столбцах которого в ходе работы подпрограммы располагаются вычисленные векторы Ланцоша;
y - вещественный двумерный массив размерности n * k, в  i - ом,  1 ≤ i ≤ kl + kr, столбце которого на выходе из подпрограммы содержится вычисленный  b - нормированный собственный вектор, соответствующий выходному собственному значению ev ( i );
ip - простая целочисленная переменная, содержащая на входе в подпрограмму признак решаемой задачи, при этом:
ip= 1 - если требуется выполнить ровно  m шагов алгоритма;
ip= 0 - если допустим выход из подпрограммы при  j < m,  j - число выполненных шагов; в этом случае вычисления будут прекращены после того, как на некотором шаге  j будет получено заданное число приближений к собственным парам задачи (1) с заданной точностью по невязке;
  на выходе из подпрограммы ip содержит признак выхода, причем:
ip= -1 - если заданное значение параметра  l не удовлетворяет неравенству  l ≥ m (9 + k) + 8k;
ip= -2 - если заданное значение параметра  k не удовлетворяет неравенству  k ≥ kl + kr + 2;
ip= -3 - если обе описанные выше ошибки имели место;
ip=  0 - если требуемые собственные значения и соответствующие собственные векторы вычислены с заданной точностью;
ip=  1 - если выполнено заданное число  m шагов алгоритма, а заданная точностью еще не достигнута;
ip=  2 - если на некотором шаге  j,  j < m, алгоритма потребовалось вычислить более  k собственных пар матрицы  Tj; в этом случае для продолжения процесса не хватает памяти;
ip=  3 - если заданная точность не может быть достигнута;
ip=  4 - если получено слишком малое значение  βj + 1;
ip=  5 - если резко возрастает объем работ по проведению выборочной ортогонализации. 
   Версии
ags0d_c - вычисление методом Ланцоша с выборочной ортогонализацией нескольких крайних собственных пар обобщенной проблемы собственных значений Ax = λ Bx, где  A и  B - вещественные симметричные матрицы,  B > 0, в случае, когда информация о матрицах  A и  B задана с двойной точностью. 
   Вызываемые подпрограммы
       avz6r_c -
       avz6d_c   		упорядочивание вектора по возрастанию его компонент с запоминанием произведенных перестановок.
       avz2r_c -
       avz2d_c   		нахождение максимальной по модулю компоненты и ее индекса из всего вектора или из заданного подмножества компонент этого вектора.
       am05r_c -
       am05d_c   		вычисление произведения матрицы на вектор. 
   Замечания по использованию
  1. 

Наименование подпрограмм, соответствующих формальным параметрам fa, fb, sb, должны быть указаны в операторе extern в программе пользователя, вызывающей ags0r_c;

  2. 

Подпрограммы fa, fb, sb должны сохранять вектор  x;

  3. 

Если начальный вектор был выбран неудачно, то возможно, что вычисленные собственные пары не обязательно будут соответствовать первым kl и последним kr собственным парам задачи (1). Поэтому, вообще говоря, для проверки необходимо обратиться к ags0r_c с разными начальными векторами. Особенно это необходимо, если на выходе из подпрограммы ip = 4 или мало выходное значение параметра  m;

  4. 

Неудачный выбор начального вектора может привести к выходу из подпрограммы на шаге с номером  j < kl + kr по признаку ip = 4. В этом случае будут вычислены все  j собственных пар, причем все  j собственных значений будут расположены в массиве ev в порядке неубывания. Понятно, что в этом случае повторное обращение к подпрограмме совершенно необходимо;

  5. 

Если на выходе из подпрограммы  ip < 0, то никакой информации об искомых собственных парах получено не будет; если же  ip > 0, то на выходе из подпрограммы будут получены kl первых и kr последних пар Ритца, вычисленных на шаге прерывания (номер этого шага - выходное значение  m), т.е. будут получены приближения к искомым собственным парам, но точность их еще не является удовлетворительной. Эти приближения могут быть использованы для построения нового начального вектора и повторного обращения к ags0r_c;

  6. 

Если на выходе из подпрограммы ip = 2, то необходимо увеличить значение параметра  k и соответственно  l. Рекомендации для выбора значения  k, которые были бы удовлетворительными для любых матриц, не могут быть даны, т.к. требуемое значение  k зависит от спектра исходной матрицы и равно числу пар Ритца, невязки которых станут меньше  √ε  || B- 1 A || к тому моменту, когда невязка для заданного числа крайних пар Ритца станут меньше eps;

  7. 

В подпрограмме ags0d_c параметры eps, vj, r, r1, r2, ev, er, v, y имеют тип double;

  8. 

Подпрограммы ags0r1_c(ags0d1_c), aes0r0_c(aes0d0_c), av04r1_c(av04d1_c) используются в качестве рабочих;

  9. 

Для подпрограммы ags0d_c в подпрограммах fa, fb и sb параметры  x и  z, должны иметь тип double;

  10.  Если для матрицы  B легко выполнимо разложение Холецкого  B = LLT, то метод Ланцоша можно применять непосредственно к неявно заданной матрице  C = L- 1 A (LT) - 1, т.е. использовать не ags0r_c, а подпрограмму, реализующую метод Ланцоша для стандартной проблемы собственных значений.

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

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

    /* Local variables */
    extern int ags0r_c(int *, int *, int *, U_fp, U_fp, U_fp, float *, int *,
                       float *, int *, int *, float *, float *, float *,
                       float *, float *, float *, float *, int *);
    static int i__, m;
    static float r__[82], v[9] /* was [3][3] */,
                         y[15] /* was [3][5] */, r1[3], r2[3];
    extern int fa_c(), fb_c(), sb_c();
    static float er[3];
    static int ip;
    static float ev[3];

#define y_ref(a_1,a_2) y[(a_2)*3 + a_1 - 4]

    m = 3;
    ip = 0;
    ags0r_c(&c__3, &c__3, &c__0, (U_fp)fa_c, (U_fp)fb_c, (U_fp)sb_c, &c_b4,
            &m, vj, &c__5, &c__82, r__, r1, r2, ev, er, v, y, &ip);

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


int fa_c(float *x, float *z__, int *n)
{
    /* Initialized data */
    static float rab[3] = { 1.f,9.f,14.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] * x[i__] - z__[i__];
/* l10: */
    }
    return 0;
} /* fa_c */

int sb_c(float *x, float *z__, int *n)
{
    /* Initialized data */
    static float rab[3] = { 1.f,1.f,.5f };

    /* 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__] = x[i__] * rab[i__ - 1];
/* l10: */
    }
    return 0;
} /* sb_c */

int fb_c(float *x, float *z__, int *n)
{
    /* Initialized data */
    static float rab[3] = { 1.f,1.f,2.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] * x[i__];
/* l10: */
    }
    return 0;
} /* fb_c */


Результаты:

       ev = (1., 7., 9.) ,
       er = (1.d - 15, 4.d - 15, 9.d - 16) ,

                    |  - 1.                5.d - 17                    - 9.d - 18 |
       y_ref =  |    2.d - 16      2.d - 15                       1.           |
                    |    2.d - 16      0.70710678118655  - 3.d - 16 |

       ip = 0
       m = 3