Текст подпрограммы и версий ( Фортран )
ams1r.zip , ams1c.zip 		Тексты тестовых примеров ( Фортран )
tams1r.zip , tams1c.zip
Текст подпрограммы и версий ( Си )
ams1r_c.zip , ams1c_c.zip 		Тексты тестовых примеров ( Си )
tams1r_c.zip , tams1c_c.zip
Текст подпрограммы и версий ( Паскаль )
ams1r_p.zip , ams1c_p.zip 		Тексты тестовых примеров ( Паскаль )
tams1r_p.zip , tams1c_p.zip

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

Назначение

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

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

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 = (y0, y1,..., yN - 1) , где вектор  x = (x0, x1,..., xN - 1) также задан.
Эта задача эквивалентна вычислению (вообще говоря) непериодической свертки двух последовательностей  as,  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 до циркулянтной матрицы А0 порядка М, которая определяется М элементами своего первого столбца: 

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

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

Далее определим вектор  x0 длины М, получаемый из заданного вектора  x добавлением К + N - 1 нулей: x0 = (x0, x1,..., xN - 1, 0, 0,..., 0). Тогда исходная задача сводится к задаче умножения циркулянтной матрицы А0 на вектор  x0, которая решается в подпрограмме АМС1R на основе алгоритма Быстрого Преобразования Фурье [2], при этом N компонент вычисленного вектора y0 = А0 x0 = (y00, y10,..., y0M - 1) будут совпадать с искомым вектором  y, а именно:  y0N - 1 + s = ys ,  s = 0, 1, ..., N - 1.

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

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

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 = (y0, y1,..., yN - 1) эквивалентна вычислению взаимной корреляции двух последовательностей  as,  xj,  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

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

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

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

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

1.  Б.Голд, Ч.Pейдер. Цифровая обработка сигналов. М., "Советское радио", 1973.
2.  В.А.Морозов, Н.Н.Кирсанова, А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов, сб. "Численный анализ на ФОPТPАНе", вып.15, Изд - во МГУ, М., 1976.
3.  В.В.Бадаева, В.А.Морозов. Алгоритмы быстрого и ускоренного решения некоторых специальных систем линейных алгебраических уравнений, Сб. "Численный анализ на ФОPТPАНе", вып.20, Изд - во МГУ, М., 1977.

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

    SUBROUTINE  AMS1R (N, M, Q, B, X) 
   Параметры
N - размерность исходной матрицы А, N ≥ 2 (тип: целый);
M - минимальная целая степень числа два, удовлетворяющая условию М ≥ 2N - 1 (тип: целый);
Q - управляющий параметр (тип: вещественный):
Q= 0. - если исходная матрица А теплицева,
Q= 1. - если исходная матрица А ханкелева;
B - вещественный одномерный массив длины М, в первых 2N - 1 элементах которого содержатся заданные определяющие элементы теплицевой матрицы А: a- (N - 1),  a- (N - 2),..., a- 1, a0, a1, ..., aN - 1 или ханкелевой матрицы А: a0, a1, ..., aN - 1, aN, ..., a2N - 3, a2N - 2;
X - вещественный одномерный массив длины М, в первых N элементах которого содержатся заданные компоненты вектора x = (x0, x1,..., xN - 1) на входе и вычисленные компоненты вектора y = (y0, y1,..., yN - 1), на выходе. 
   Версии
AMS1C - вычисление произведения комплексной теплицевой или ханкелевой матрицы А на комплексный вектор  x. Для подпрограммы АМS1С параметры B, X имеют тип СОМРLЕХ. 
   Вызываемые подпрограммы
 AMC1R -
 AMC1C   		подпрограммы умножения вещественной циркулянтной матрицы на вещественный вектор и компексной циркулянтной матрицы на комплексный вектор на основе алгоритма Быстрого Преобразования Фурье (используются в АМS1R и АМS1С соответственно). 
   Замечания по использованию
  1. 

После выхода из подпрограмм АМS1R и АМS1С значения определяющих элементов исходной матрицы А в массиве B не сохраняются.

  2. 

Если матрица А циркулянтная (частный случай теплицевой),то целесообразно вместо подпрограмм АМS1R, АМS1С воспользоваться соответствующими подпрограммами АМС1R и АMC1C.

  3. 

Если теплицева матрица А имеет вид 

          | 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     |
и задается 2N - 1 элементами первого столбца и первой строки:  a- (N - 1), a- 1, a0, a1, a2,..., aN - 1; то до обращения к подпрограммам АМS1R, АМS1С эти элементы следует упорядочить в массиве B по убыванию индексов:  aN - 1, ..., a1, a0, a- 1, a- 2, ..., a- (N - 1);

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

   1.
      DIMENSION  B(16), X(16)
      DATA  B /6., 5., 4., 3., 2., 1., -1., 0., 2., -3., 7., 0., 0., 0., 0., 0./
      DATA  X /1., -1., -2., 3., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0./
      N = 6
      M = 16
      Q = 0.
      CALL  AMS1R (N, M, Q, B, X)

Результат:     y  =  ( 17., 13., 13., 13., -4., 8. )

   2.
      DIMENSION  B(16), X(16)
      DATA  B /1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 0., 0., 0., 0., 0./
      DATA  X /1., -1., -2., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0./
      N = 6
      M = 16
      Q = 1.
      CALL  AMS1R (N, M, Q, B, X)

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

   3.
      COMPLEX  B(8), X(8)
      DATA  B /(4., -1.), (3., 2.), (2., 0.), (1., 1.), (-1., 0.), (0., -1.), (2., 1.), (0., 0.)/
      DATA  X /(1., -1.), (-1., 0.), (-2., 1.), (3., -2.), (0., 0.), (0., 0.), (0., 0.), (0., 0.)/
      N = 4
      M = 8
      Q = 0.
      CALL  AMS1C (N, M, Q, B, X)

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

   4.
      COMPLEX  B(8), X(8)
      DATA  B /(1., 1.), (2., 0.), (3., 2.), (4., -1.), (5., 0.), (6., 1.), (7., 0.), (0., 0.)/
      DATA  X /(1., -1.), (0., 1.), (1., 0.), (-1., 1.), (0., 0.), (0., 0.), (0., 0.), (0., 0.)/
      N = 4
      M = 8
      Q = 1.
      CALL  AMS1C (N, M, Q, B, X)

Результат:     y  =  ( (2., 9.), (-1., 5.), (4., 8.), (2., 8.) )