БЧА НИВЦ МГУ. FTA6R. Ускоренные дискретные преобразования Фурье

Текст подпрограммы и версий ( Фортран ) fta6r.zip , fta6d.zip	Тексты тестовых примеров ( Фортран ) tfta6r.zip , tfta6d.zip
Текст подпрограммы и версий ( Си ) fta6r_c.zip , fta6d_c.zip	Тексты тестовых примеров ( Си ) tfta6r_c.zip , tfta6d_c.zip
Текст подпрограммы и версий ( Паскаль ) fta6r_p.zip , fta6e_p.zip	Тексты тестовых примеров ( Паскаль ) tfta6r_p.zip , tfta6e_p.zip

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

Назначение

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

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

Пусть известны значения h_{k1,k2,...,kNDIM} = h (t_1,k1, t_2,k2, ..., t_NDIM,kNDIM) комплексно - значной функции h (t₁, t₂, ..., t_NDIM) от NDIM вещественных переменных на многомерной сетке размерности NDIM, узлы которой по каждой размерности (т.е. по каждой переменной) определяются следующим образом:

t_1,k1 = k1 Δt₁ ( k1 = 0, 1, 2, ..., N₁ - 1; Δt₁ - шаг сетки по первой размерности; N₁ равняется целой степени двух и определяет число узлов сетки по первой размерности)

t_2,k2 = k2 Δt₂ ( k2 = 0, 1, 2, ..., N₂ - 1; Δt₂ - шаг сетки по второй размерности; N₂ равняется целой степени двух и определяет число узлов сетки по второй размерности )
   .
   .
   .

t_NDIM,kNDIM = kNDIM Δt_NDIM ( kNDIM = 0, 1, 2, ..., N_NDIM ; Δt_NDIM - шаг сетки по последней размерности NDIM; N_NDIM равняется целой степени двух и определяет число узлов сетки по последней размерности NDIM )

Эти комплексные значения h_{k1,k2,...,kNDIM} задаются в виде одномерного вещественного массива ARRAY длины 2 * N₁ * N₂ * ... * N_NDIM (здесь * означает знак умножения) по следующему правилу (I)
(Re означает вещественную часть, а Im - мнимую часть комплексного числа):

   Re(h_{k1,k2,...,kNDIM}) ® ARRAY( 1 + 2(k1 + N₁*k2 + N₁*N₂*k3 + ... +
                                                     + N₁*N₂*...*N_NDIM-1*kNDIM) )

   Im(h_{k1,k2,...,kNDIM}) ® ARRAY( 2 + 2(k1 + N₁*k2 + N₁*N₂*k3 + ... +
                                                     + N₁*N₂*...*N_NDIM-1*kNDIM) )

Таким образом, исходные комплексные значения представляются в виде вещественного одномерного массива по правилам Фортрана, определенным для представления комплексных многомерных массивов в виде вещественных одномерных массивов (см. также "Замечания по использованию").

Информация о количестве узлов NDIM - мерной сетки по каждой размерности указывается в компонентах целого массива MESH длины NDIM следующим образом:

               MESH  =  (N₁, N₂, ..., N_NDIM)

Подпрограмма FTA6R имеет два режима работы, задаваемых при обращении к ней значением параметра IREG. В первом режиме (IREG = 1) выполняется NDIM - мерное прямое быстрое дискретное преобразование Фурье заданного указанным выше образом массива ARRAY, состоящее в получении N₁ * N₂ * ... * N комплексных чисел

H_{n1,n2,...,nNDIM} ( n1 = 0, 1, 2, ..., N₁ - 1, n2 = 0, 1, 2, ..., N₂ - 1, ..., nNDIM = 0, 1, 2, ..., N_NDIM - 1 ) по формуле:

H_{n1,n2,...,nNDIM}  =

                         N_NDIM-1        N₁-1
                     =       ∑          ...    ∑       exp(2πi * kNDIM * nNDIM / N_NDIM) *
                           ^kNDIM=0       ^k1=0

                                                * ... * exp(2πi * k1 * n1 / N₁) h_{k1,k2,...,kNDIM}

Вычисленные значения H_{n1,n2,...,nNDIM} располагаются в том же массиве ARRAY по правилу (I).

Во втором режиме (IREG = - 1) выполняется NDIM - мерное обратное быстрое дискретное преобразование Фурье, состоящее в восстановлении из значений H_{n1,n2,...,nNDIM}, заданных по правилу (I) в одномерном вещественном массиве ARRAY длины 2 * N₁ * N₂ * ... * N_NDIM, значений h_{k1,k2,...,kNDIM} по формуле

h_{k1,k2,...,kNDIM}  =  1 / (N₁*N₂*...*N_NDIM) *

                          N_NDIM-1      N₁-1
                      *      ∑          ...   ∑      exp(- 2πi * nNDIM * kNDIM / N_NDIM) *
                         ^nNDIM=0        ⁿ¹⁼⁰

                                               * ... * exp(- 2πi * n1 * k1 / N₁) H_{n1,n2,...,nNDIM}

Вычисленные значения h_{k1,k2,...,kNDIM} располагаются в том же массиве ARRAY по правилу (I).

Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков. Численные методы. Изд - во "Наука", 1987.

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

    SUBROUTINE  FTA6R (ARRAY, NDIM, MESH, IREG)

   Параметры

ARRAY -	вещественный одномерный массив длины 2 * N₁ * N₂ * ... * N_NDIM, содержащий на входе h_{k1,k2,...,kNDIM} (IREG = 1) или H_{n1,n2,...,nNDIM} (IREG = - 1), а на выходе значения H_{n1,n2,...,nNDIM} (IREG = 1) или h_{k1,k2,...,kNDIM} (IREG = - 1);
NDIM -	число размерностей заданной многомерной сетки (тип: целый);
MESH -	целый массив, i - ая компоненты которого содержит число узлов заданной NDIM - мерной сетки по i - му направлению ( i = 1, 2, 3, ..., NDIM), т.е. MESH = (N₁, N₂, ..., N_NDIM); каждое значение N_i должно быть целой степенью двух;
IREG -	задает режим работы подпрограммы (тип: целый); при этом:

IREG= 1 -	когда выполняется прямое преобразование;
IREG= -1 -	когда выполняется обратное преобразование.

   Версии

FTA6D -

выполнение прямого и обратного быстрого дискретного преобразования Фурье многомерного массива комплексных чисел, длина каждой размерности которого равняется целой степени двух, в режиме удвоенной точности; при этом параметр ARRAY должен иметь тип DOUBLE PRECISION.

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

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

	1.	В подпрограммах FTA6R и FTA6D проверка того, что каждое из значений N₁, N₂, ..., N_NDIM должно быть степенью двух, не производится.
	2.	В правиле (I) предполагается, что k_i меняются от 0 до N_i - 1 ( i = 0, 1, 2, ..., NDIM). Поскольку в Фортране индексы массивов не могут начинаться с 0, а только с 1, то при реальной засылке по правилу (I) значений h_{k1,k2,...,kNDIM} в массив ARRAY значения k_i должны меняться от 1 до N_i. Это означает, что при такой реальной засылке правило (I) может быть непосредственно применено если при каждом использовании из k_i вычитать 1 (см. также "Пример использования").

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

       DIMENSION  ARRAY(32), MESH(2) 
       NDIM = 2 
       N1 = 4 
       N2 = 4 
       MESH(1) = N1 
       MESH(2) = N2 
       DO 1  K2 = 1, N1 
       DO 1  K1 = 1, N2 
       I = 1 + 2*( (K1 - 1) + N1*(K2 - 1) )
       ARRAY(I) = 4.0*K1 + K2 - 4.0
       ARRAY(I + 1) = 0.0
    1 CONTINUE
       CALL  FTA6R(ARRAY, NDIM, MESH, 1)

Результаты:

    ARRAY = (136, 0, - 32, - 32, - 32, 0, - 32, 32, - 8, - 8, 0, 0, 0, 0,
                       0, 0, - 8, 0, 0, 0, 0, 0, 0, 0, - 8, 8, 0, 0, 0, 0, 0, 0)