Fortran 中具有 FFTW 的复杂数组的 DCT：如何指向虚部数组？
DCT of complex arrays with FFTW in Fortran: How to point to the imaginary part array?
我正在用 Fortran 编写伪谱 CFD 代码，它本质上是平面层中 Navier-Stokes 方程的时间步进器。在我的案例中，它确实是一个 3d 代码，但问题在 2d 中可以很好地理解，所以我会坚持这种情况。在几何上，我的 2d 平面层以 y=0 和 y=1 为界，并且沿另一个方向 x 是周期性的。在不深入研究杂草的情况下，一种有效的离散化是沿着 y 分解切比雪夫多项式上的场（例如速度），沿着 x 分解傅里叶模式。切比雪夫多项式本质上是伪装在扭曲网格上的余弦。 Navier-Stokes 方程在谱 space 中具有简单形式，但非线性项除外。因此，大部分计算都是在频谱 space 中进行的，偶尔会在物理 space 中进行偏移以计算非线性项：需要将复杂的 Chebyshev-Fourier 系数的二维数组转换为相应的二维字段（即网格上的实数值数组）。在没有其他限制的情况下，这种变换相对容易实现。例如，从复谱系数开始——让我们称它们为 c_in(NX, NY/2+1)——人们可以对 y 的每个值沿 x 进行复数到实数的离散傅里叶变换以获得real Chebyshev 系数的二维数组。然后，可以对所有 x 和 voila 沿着 y 执行离散余弦变换（FFTW 中的 FFTW_REDFT），最终得到真实的场, r_out(NX,NY).
万事俱备的根源在于，由于某些原因，我需要先计算DCT。这是一个问题，因为余弦变换仅针对 FFTW 中的真实数据实现。各种限制导致我不想将复杂的频谱系数数组拆分为实部和虚部的两个实数数组。鉴于这些对数据结构的限制，问题是：如何有效地让 FFTW 沿着复数数组的第一个索引计算多个 DCT。
到目前为止，我的解决方案包括使用 plan_many_r2r 高级接口来定义一个跨越虚数值的变换：我将 idist 设置为 2。因此，如果我将此计划与与 c_in(1,1) 的实部关联的指针 ptr2real_in 一起使用，计算所有实部的余弦变换。然后，我使用与 c_in(1,1) 的 虚数 部分关联的指针 ptr2imag_in 重复执行计划。之后，沿第二个维度计算复数到实数的 DFT 就很容易了。
所以这个方法的关键是定义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
一个人总是可以 transfer() 指向整数的指针，进行移位并 transfer() 返回，如果这是你想要的。
cptr = transfer( transfer(cptr, 1_c_intptr_t) + c_sizeof(1._c_double) , 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
我正在用 Fortran 编写伪谱 CFD 代码，它本质上是平面层中 Navier-Stokes 方程的时间步进器。在我的案例中，它确实是一个 3d 代码，但问题在 2d 中可以很好地理解，所以我会坚持这种情况。在几何上，我的 2d 平面层以 y=0 和 y=1 为界，并且沿另一个方向 x 是周期性的。在不深入研究杂草的情况下，一种有效的离散化是沿着 y 分解切比雪夫多项式上的场（例如速度），沿着 x 分解傅里叶模式。切比雪夫多项式本质上是伪装在扭曲网格上的余弦。 Navier-Stokes 方程在谱 space 中具有简单形式，但非线性项除外。因此，大部分计算都是在频谱 space 中进行的，偶尔会在物理 space 中进行偏移以计算非线性项：需要将复杂的 Chebyshev-Fourier 系数的二维数组转换为相应的二维字段（即网格上的实数值数组）。在没有其他限制的情况下，这种变换相对容易实现。例如，从复谱系数开始——让我们称它们为 c_in(NX, NY/2+1)——人们可以对 y 的每个值沿 x 进行复数到实数的离散傅里叶变换以获得real Chebyshev 系数的二维数组。然后，可以对所有 x 和 voila 沿着 y 执行离散余弦变换（FFTW 中的 FFTW_REDFT），最终得到真实的场, r_out(NX,NY).
万事俱备的根源在于，由于某些原因，我需要先计算DCT。这是一个问题，因为余弦变换仅针对 FFTW 中的真实数据实现。各种限制导致我不想将复杂的频谱系数数组拆分为实部和虚部的两个实数数组。鉴于这些对数据结构的限制，问题是：如何有效地让 FFTW 沿着复数数组的第一个索引计算多个 DCT。
到目前为止，我的解决方案包括使用 plan_many_r2r 高级接口来定义一个跨越虚数值的变换：我将 idist 设置为 2。因此，如果我将此计划与与 c_in(1,1) 的实部关联的指针 ptr2real_in 一起使用，计算所有实部的余弦变换。然后，我使用与 c_in(1,1) 的 虚数 部分关联的指针 ptr2imag_in 重复执行计划。之后，沿第二个维度计算复数到实数的 DFT 就很容易了。
所以这个方法的关键是定义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
一个人总是可以 transfer() 指向整数的指针，进行移位并 transfer() 返回，如果这是你想要的。
cptr = transfer( transfer(cptr, 1_c_intptr_t) + c_sizeof(1._c_double) , 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