DCT сложных массивов с FFTW в Фортране: как указать на массив мнимой части?

Я пишу псевдоспектральный код CFD на Фортране, который, по сути, является шаговым шагом по уравнениям Навье-Стокса в плоском слое. В моем случае это действительно трехмерный код, но проблему можно хорошо понять в 2d, поэтому я буду придерживаться этого случая. Геометрически мой 2-мерный плоский слой ограничен y=0 а также y=1и является периодическим вдоль другого направления x, Не вдаваясь слишком много в сорняки, эффективная дискретизация состоит в том, чтобы разложить поля (например, скорость) на полиномы Чебышева вдоль y и моды Фурье вдоль x, Чебышевские многочлены по существу являются замаскированными косинусами на искаженной сетке. Уравнения Навье-Стокса имеют простую форму в спектральном пространстве, за исключением нелинейного члена. Таким образом, большинство вычислений выполняется в спектральном пространстве со случайными отклонениями в физическом пространстве для вычисления нелинейного члена: необходимо преобразовать двумерный массив комплексных коэффициентов Чебышева-Фурье в соответствующее двумерное поле (то есть массив вещественных значений на сетке). Без других ограничений это преобразование относительно легко реализовать. Например, начиная со сложных спектральных коэффициентов - назовем их c_in(NX, NY/2+1) - можно взять комплекс к реальному, дискретному преобразованию Фурье по x для каждого значения y получить двумерный массив вещественных чебышевских коэффициентов. Затем можно выполнить дискретное косинусное преобразование (FFTW_REDFT в FFTW) вместе y для всех x и, вуаля, наконец-то получается настоящее поле, r_out(NX,NY),

Корень всех проблем в том, что по некоторым причинам мне нужно сначала вычислить DCT. Это проблема, потому что косинусные преобразования реализуются только для реальных данных в FFTW. Различные ограничения приводят к тому, что я не хочу разбивать свой сложный массив спектральных коэффициентов на два реальных массива для действительной и мнимой частей. Учитывая эти ограничения на структуру данных, возникает вопрос: как эффективно заставить FFTW вычислять несколько DCT вдоль первого индекса массива комплексных чисел.

Пока что мое решение состоит в использовании plan_many_r2r расширенный интерфейс для определения преобразования, которое перепрыгивает через мнимые значения: я установил idist 2. В результате, если я использую этот план с указателем ptr2real_in связано с реальной частью c_in(1,1)вычисляется косинус-преобразование всех вещественных частей. Затем я повторяю выполнение плана с указателем ptr2imag_in связано с мнимой частью c_in(1,1), После этого вычисление комплексного и реального ДПФ по второму измерению становится простым.

Так что суть этого подхода заключается в определении ptr2imag_in, который на самом деле является адресом памяти c_in(1,1) смещен на размер в памяти C_double, Ниже приведен минимальный пример, который работает, но выглядит для меня неуклюже. В частности, я определяю указатель на мнимую часть сложного массива таким образом

cptr = C_loc(c_in(1,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])
cptr = C_loc(ptr2imag_in(2,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])

Мне кажется, что все, что мне нужно сделать, это сдвинуться cptr на 8 байт. Как я мог это сделать? Сбой следующего кода:

cptr = C_loc(c_in(1,1))
cptr = cptr + 8 
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])

Ниже приведен полный минимальный пример взятия DCT, сопровождаемого комплексным и реальным DFT:

Program monkeying_with_dct 
Use, Intrinsic :: ISO_C_BINDING
Implicit None
include 'fftw3.f03'

Integer, Parameter :: dp = C_Double
Complex(C_double), Pointer :: c_in (:,:)
Complex(C_double), Pointer :: c_out(:,:)
Real(C_Double),    Pointer :: r_out(:,:)
Real(C_Double),    Pointer :: ptr2real_in (:,:)
Real(C_Double),    Pointer :: ptr2real_out(:,:)
Real(C_Double),    Pointer :: ptr2imag_in (:,:)
Real(C_Double),    Pointer :: ptr2imag_out(:,:)
Type(C_ptr) :: cptr
Type(C_ptr) :: plan_IDCT
Type(C_ptr) :: plan_C2R 
Type(C_ptr) :: pdum


Integer, Parameter :: NX = 512
Integer, Parameter :: NY = 1024

print *,'... allocating memory ...'
pdum = fftw_alloc_complex(int((NY/2+1)*NX, C_size_T))
Call C_F_Pointer(pdum, c_in , [NX, NY/2+1])
pdum = fftw_alloc_complex(int((NY/2+1)*NX, C_size_T))
Call C_F_Pointer(pdum, c_out, [NX, NY/2+1])
pdum = fftw_alloc_real(int(NY*NX, C_size_T))
Call C_F_Pointer(pdum, r_out, [NX, NY])

print *,'... initializing data ...'
c_in      = Cmplx(0._dp, 0._dp,  Kind=dp)
c_in(2,3) = Cmplx(1._dp, 0.5_dp, Kind=dp) 

print *, '... defining a pointer to the real part of input complex data ...'                
cptr = C_loc(c_in(1,1))
Call C_F_Pointer(cptr, ptr2real_in, [2*NX, (NY/2+1)])
print *, '... defining a pointer to the imag part of input complex data ...'                
cptr = C_loc(c_in(1,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])
cptr = C_loc(ptr2imag_in(2,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])

print *, '... defining a pointer to the real part of transformed complex data ...'                
cptr = C_loc(c_out(1,1))
Call C_F_Pointer(cptr, ptr2real_out, [2*NX, (NY/2+1)])
print *, '... defining a pointer to the imag part of transformed complex data ...'                
cptr = C_loc(c_out(1,1))
Call C_F_Pointer(cptr, ptr2imag_out, [2*NX, (NY/2+1)])
cptr = C_loc(ptr2imag_out(2,1))
Call C_F_Pointer(cptr, ptr2imag_out, [2*NX, (NY/2+1)])


print*, '... planning IDCT ...'
plan_IDCT = fftw_plan_many_r2r(1, [NX], (NY/2+1), &
                        ptr2real_in,  [2*NX], 2, 2*NX, &
                        ptr2real_out, [2*NX], 2, 2*NX, &
                        [FFTW_REDFT01]  , FFTW_MEASURE) 

print*, '... planning C2R DFT ...' 
plan_C2R   = fftw_plan_many_dft_c2r(1, [NY], NX, & 
                        c_out, [NY/2+1], NX, 1, & 
                        r_out, [NY], NX, 1, & 
                        FFTW_MEASURE)  

print*, '... DCT of the real part ...' 
Call fftw_execute_r2r(plan_IDCT, ptr2real_in, ptr2real_out) 
print*, '... DCT of the imaginary part ...' 
Call fftw_execute_r2r(plan_IDCT, ptr2imag_in, ptr2imag_out) 
print*, '... DFT Complex to real ...' 
Call fftw_execute_dft_c2r(plan_C2R, c_out,r_out) 


End Program monkeying_with_dct
Источник

2 ответа

Решение

Можно уже transfer() указатель на целое число, сделать сдвиг и transfer() назад, если это то, что вы хотите.

cptr = transfer( transfer(cptr, 1_c_intptr_t) + 8_c_intptr_t  , cptr)

или можно просто вызвать маленький C, который делает арифметику указателей более управляемым способом. Но я не уверен, что это действительно то, что вам нужно.

В Fortran 2008 можно просто использовать %im синтаксис, но поддержка компилятора пока не очень хорошая, gfortran его вообще не поддерживает..

Источник

Вы можете попробовать что-то в смысле следующего пробного 1-dim примера (я не уверен, насколько эффективен скомпилированный код...).

program mapctor
 implicit none
 integer, parameter :: n = 10
 integer            :: i

 complex            :: cx(n)   ! the couplex array
 real               :: rx(2,n) ! mapped to real

 EQUIVALENCE(cx(1), rx(1,1))   ! some old fortran stuff !

! fill in some test values
 do i=1,n
    cx(i) = cmplx(i,-i)
 enddo


 write(6,*)"REAL PART:"
 call xfft(rx(1,:),n)   

 write(6,*)"IMAG PART:"
 call xfft(rx(2,:),n)

end program mapctor


subroutine xfft(x,n)  ! mock-up fft routine, or whatever
  implicit none
  real,    intent(in) :: x(n)
  integer, intent(in) :: n

  write(6,'(f12.6)') x

end subroutine xfft
Источник
Другие вопросы по тегам