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

Назначение

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

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

1. Задана теплицева матрица А размерности N*N:

                   | a0      a-1    a-2    ...  a-(N-2)   a-(N-1)  |
                   | a1      a0     a-1    ...  a-(N-3)   a-(N-2)  |
                   | a2      a1      a0    ...  a-(N-4)   a-(N-3)  |
(1)              | . . . . . . . . . . . . . . . . . . . . . . . . . . . .           ,
                   | aN-2  aN-3  aN-4  ...  a0         a-1        |
                   | aN-1  aN-2  aN-3  ...  a1         a0         |

которая определяется 2n - 1 элементами первой строки и первого столбца:

         a -(N-1), a -(N-2), ... , a -2, a -1, a0, a1, a2, ... , aN-2, aN-1 ; 

все элементы этой матрицы на любой ее диагонали, параллельной главной, равны между собой. Требуется вычислить вектор y = Ax = (y₀, y₁,..., y_N - 1) , где вектор  x = (x₀, x₁,..., x_N - 1) также задан.
Эта задача эквивалентна вычислению (вообще говоря) непериодической свертки двух последовательностей  a_s,  x _j,  s = - (N - 1), ..., - 1, 0, 1, ..., N - 1,   j = 0, 1,..., N - 1 разной длины 2N - 1 и N соответственно [1], т.е.

                            N-1 
            ys      =      ∑     as -j x j  ,     s = 0, 1, ... , N-1 ,
                             j=0 

Дополним теплицеву матрицу А порядка N до циркулянтной матрицы А⁰ порядка М, которая определяется М элементами своего первого столбца:

         a -(N-1), a -(N-2), ... , a -1, a0, a1, ... , aN-1, 0, 0, ... , 0 , 

где М - наименьшая целая степень числа два, удовлетворяющая условию М ≥ 2N - 1, а число нулей K равно разности: К = М - (2N - 1) ≥  0.

Далее определим вектор  x⁰ длины М, получаемый из заданного вектора  x добавлением К + N - 1 нулей: x⁰ = (x₀, x₁,..., x_N - 1, 0, 0,..., 0). Тогда исходная задача сводится к задаче умножения циркулянтной матрицы А⁰ на вектор  x⁰, которая решается в подпрограмме amc1r_c на основе алгоритма Быстрого Преобразования Фурье [2], при этом N компонент вычисленного вектора y⁰ = А⁰ x⁰ = (y₀⁰, y₁⁰,..., y⁰_M - 1) будут совпадать с искомым вектором  y, а именно:  y⁰_{N - 1 + s} = y_s ,  s = 0, 1, ..., N - 1.

Аналогичный метод умножения теплицевой матрицы на вектор изложен в работе [3].

Реализованный в подпрограмме алгоритм требует памяти ЭВМ не менее 4N и не более 8N (из которых 3N занимает исходная информация). Время вычислений пропорционально М log₂ М, т.е. по порядку не превосходит 4N (2 + log₂N).

2. Пусть А - ханкелева матрица вида

               | a0         a1        a2     ...  aN-2    aN-1  |
               | a1         a2        a3     ...  aN-1    aN     |
               | a2         a3        a4     ...  aN       aN+1 |
(1)          | . . . . . . . . . . . . . . . . . . . . . . . . . . .                   ,
               | aN-2      aN-1    aN     ... a2N-4   a2N-3 |
               | aN-1      aN       aN+1 ... a2N-3   a2N-2 |

которая определяется 2N - 1 элементами первой строки и последнего столбца:

         a0, a1, ... , aN-1, aN, ... , a2N-3, a2N-2 ; 

все элементы этой матрицы на любой ее диагонали, параллельной побочной, равны между собой. В данном случае задача вычисления вектора y = Ax = (y₀, y₁,..., y_N - 1) эквивалентна вычислению взаимной корреляции двух последовательностей  a_s,  x_j,  s = 0, 1,..., 2N - 2,  j = 0, 1,..., N - 1 длины 2N - 1 и N соответственно [1], т.е.

                            N-1 
            ys      =      ∑     as+j x j ,    s = 0, 1, ... , N-1 ,
                             j=0 

в частности, при  x_j = a_j,  j = 0, 1,..., N - 1 вектор  y = Аx определяет автокорреляцию последовательности  a_s,  s = 0, 1, ..., 2N - 2.

Легко видеть, что ханкелева матрица А переходит в теплицеву T путем умножения слева на матрицу перестановок U вида

                               |                   1 |
                               |                1    |
                               |             .        |
                     U =    |          .           |    =   U -1 ,
                               |       .              |
                               |    1                |
                               | 1                   | 

т.е. Т = UА. Поэтому для получения вектора  y = Аx сначала вычисляется вектор Tx (по описанному выше алгоритму) и затем делается соответствующая перестановка его компонент, а именно:  y = U(Тx).

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

    int ams1r_c (integer *n, integer *m, real *q, real *b, real *x)

   Параметры 

   Версии 

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

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

          | a0         a1        a2         ...   aN-1  |
          | a-1        a0        a1         ...   aN-2  |
          | a-2        a-1       a0         ...   aN-3  |
          | . . . . . . . . . . . . . . . . . . . . . . . . .
          | a-(N-1)   a-(N-2)  a-(N-3)  ...  a0     | 

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

  
1.
int main(void)
{
    /* Initialized data */
    static float b[16] = { 1.f,2.f,3.f,4.f,5.f,6.f,7.f,8.f,9.f,10.f,11.f,0.f,
                           0.f,0.f,0.f,0.f };
    static float x[16] = { 1.f,-1.f,-2.f,0.f,1.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
                           0.f,0.f,0.f,0.f };

    /* Local variables */
    extern int ams1r_c(int *, int *, float *, float *, float *);
    static int i__, m, n;
    static float q;

    n = 6;
    m = 16;
    q = 1.f;
    ams1r_c(&n, &m, &q, b, x);

    for (i__ = 1; i__ <= 6; ++i__) {
         printf("\n  %16.7e \n", x[i__ - 1]);
    }
/* L1: */
    return 0;
} /* main */


Результат:     y  =  ( -2., -3., -4., -5., -6., -7. )

2.
int main(void)
{
    /* Initialized data */
    static complex b[8] = { {4.f,-1.f},{3.f,2.f},{2.f,0.f},{1.f,1.f},
                            {-1.f,0.f},{0.f,-1.f},{2.f,1.f},{0.f,0.f} };
    static complex x[8] = { {1.f,-1.f},{-1.f,0.f},{-2.f,1.f},{3.f,-2.f},
                            {0.f,0.f},{0.f,0.f},{0.f,0.f},{0.f,0.f} };

    /* Local variables */
    extern int ams1c_c(int *, int *, float *, complex *, complex *);
    static int i__, m, n;
    static float q;

    n = 4;
    m = 8;
    q = 0.f;
    ams1c_c(&n, &m, &q, b, x);

    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n  %20.12e %20.12e \n", x[i__ - 1].r, x[i__ - 1].i);
    }
/* L1: */
    return 0;
} /* main */

Результат:     y  =  ( (2., -12.), (7., 2.), (3., -6.), (10., 0.) )