hipFFTW examples

hipFFTW examples#

hipFFTW is the FFTW3-compatible interface shipped with hipFFT. The routine names, planner flags and calling sequence are those of FFTW3, so existing FFTW code moves across with little change. hipFORT exposes it through the hipfort_hipfftw module.

The one difference that matters: the in and out arguments are device pointers. FFTW declares them void*, so a pointer from hipMalloc passes straight through, but host arrays do not work.

Every program on this page is a complete, self-contained example that is built and run as part of the hipFORT test suite. The tests live in test/f2003/hipfftw. Unlike the other FFT libraries there is no Fortran 2008 variant, because the FFTW API is pointer-based throughout and gains nothing from Fortran array pointers.

For the cuFFT-style interface to the same library, see hipFFT examples.

Transform workflow#

A hipFFTW transform follows the FFTW3 sequence:

  1. Allocate device memory with hipMalloc, or host-accessible memory with fftw_alloc_real and fftw_alloc_complex.

  2. Build a plan with fftw_plan_dft_1d, fftw_plan_dft_r2c_1d, fftw_plan_many_dft or fftw_plan_guru_dft.

  3. Run it with the matching fftw_execute_dft, fftw_execute_dft_r2c or fftw_execute_dft_c2r.

  4. Release the plan with fftw_destroy_plan.

Keep the following conventions in mind:

  • FFTW transforms are unnormalized. A forward transform followed by an inverse transform of length N returns N times the original data.

  • Planner flags are the standard FFTW values. FFTW_ESTIMATE is not emitted by the generated enums module, so declare it yourself as integer(c_int), parameter :: FFTW_ESTIMATE = 64.

  • Real forward transforms produce Hermitian-symmetric output, so only N/2 + 1 complex values are stored.

  • Multi-dimensional transforms use C row-major order, so the last dimension is contiguous.

  • The double precision routines are named fftw_* and the single precision ones fftwf_*.

Building an example#

The examples only need the hipfftw and hip hipFORT components:

find_package(hipfort REQUIRED COMPONENTS hip hipfftw)

add_executable(my_fft hipfftw_c2c.f03)
target_link_libraries(my_fft PRIVATE hipfort::hipfftw hipfort::hip)

See Using hipFORT in your application for the full set of build options.

Complex-to-complex transform#

A one-dimensional complex-to-complex transform. The plan is built with fftw_plan_dft_1d over two device pointers and executed once. The input is a sum of two harmonics, so the output has energy in exactly two bins.

program hipfftw_c2c_test
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftw

  implicit none

  integer(c_int), parameter :: N = 16
  integer(c_size_t), parameter :: Nbytes = N * 16  ! sizeof(double complex) = 16
  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-12
  ! FFTW_ESTIMATE is not emitted by the generated enums module; use the
  ! standard FFTW planner-flag value.
  integer(c_int), parameter :: FFTW_ESTIMATE = 64

  complex(c_double_complex), allocatable, target, dimension(:) :: hx, hresult
  type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr
  type(c_ptr) :: plan = c_null_ptr
  integer :: j
  double precision :: error, max_error
  complex(c_double_complex) :: w

  write(*,"(a)",advance="no") "-- Running test 'hipfftw_c2c' (Fortran 2003) - "

  ! Signal: x[j] = exp(2*pi*i*1*j/N) + 2*exp(2*pi*i*5*j/N)
  ! Expected DFT: X[1]=N, X[5]=2*N, all other bins zero (0-indexed)
  w = cmplx(0.0d0, 2.0d0 * pi / dble(N), kind=c_double_complex)
  allocate(hx(N), hresult(N))
  do j = 0, N-1
    hx(j+1) = exp(w * dble(j)) + 2.0d0 * exp(5.0d0 * w * dble(j))
  end do

  call hipCheck(hipMalloc(dx, Nbytes))
  call hipCheck(hipMalloc(dy, Nbytes))
  call hipCheck(hipMemcpy(dx, c_loc(hx(1)), Nbytes, hipMemcpyHostToDevice))

  ! Device pointers pass directly to the FFTW C API (in/out are void*).
  plan = fftw_plan_dft_1d(N, dx, dy, FFTW_FORWARD, FFTW_ESTIMATE)
  call fftw_execute_dft(plan, dx, dy)
  call fftw_destroy_plan(plan)

  call hipCheck(hipMemcpy(c_loc(hresult(1)), dy, Nbytes, hipMemcpyDeviceToHost))

  max_error = 0.0d0
  do j = 0, N-1
    if (j == 1) then
      error = abs(hresult(j+1) - cmplx(dble(N), 0.0d0, kind=c_double_complex))
    else if (j == 5) then
      error = abs(hresult(j+1) - cmplx(2.0d0 * N, 0.0d0, kind=c_double_complex))
    else
      error = abs(hresult(j+1))
    end if
    max_error = max(max_error, error)
  end do

  if (max_error > tol) then
    write(*,*) "FAILED! max error = ", max_error
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  deallocate(hx, hresult)

  write(*,*) "PASSED!"
end program hipfftw_c2c_test

Real-to-complex and complex-to-real transforms#

fftw_plan_dft_r2c_1d and fftw_plan_dft_c2r_1d build the two halves of a real round trip. The complex side holds N/2 + 1 elements.

program hipfftw_r2c_c2r_test
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftw

  implicit none

  integer(c_int), parameter :: N = 16
  integer(c_int), parameter :: Nc = N / 2 + 1
  integer(c_size_t), parameter :: Nbytes_r = N * 8    ! sizeof(double) = 8
  integer(c_size_t), parameter :: Nbytes_c = Nc * 16  ! sizeof(double complex) = 16
  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-12
  ! FFTW_ESTIMATE is not emitted by the generated enums module; use the
  ! standard FFTW planner-flag value.
  integer(c_int), parameter :: FFTW_ESTIMATE = 64

  real(c_double), allocatable, target, dimension(:) :: hx, hresult_r
  complex(c_double_complex), allocatable, target, dimension(:) :: hresult_c
  type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr, dz = c_null_ptr
  type(c_ptr) :: plan = c_null_ptr
  integer :: j
  double precision :: error, max_error

  write(*,"(a)",advance="no") "-- Running test 'hipfftw_r2c_c2r' (Fortran 2003) - "

  ! Real signal: x[j] = 1 + 2*cos(2*pi*j/N) + 3*cos(2*pi*3*j/N)
  ! R2C output (Nc = N/2+1 complex values, 0-indexed):
  !   X[0] = N, X[1] = N, X[3] = 3*N/2, rest zero
  allocate(hx(N), hresult_c(Nc), hresult_r(N))
  do j = 0, N-1
    hx(j+1) = 1.0d0 + 2.0d0 * cos(2.0d0 * pi * j / dble(N)) &
                     + 3.0d0 * cos(6.0d0 * pi * j / dble(N))
  end do

  call hipCheck(hipMalloc(dx, Nbytes_r))
  call hipCheck(hipMalloc(dy, Nbytes_c))
  call hipCheck(hipMemcpy(dx, c_loc(hx(1)), Nbytes_r, hipMemcpyHostToDevice))

  ! Forward R2C. Device pointers pass directly to the FFTW C API.
  plan = fftw_plan_dft_r2c_1d(N, dx, dy, FFTW_ESTIMATE)
  call fftw_execute_dft_r2c(plan, dx, dy)
  call fftw_destroy_plan(plan)

  call hipCheck(hipMemcpy(c_loc(hresult_c(1)), dy, Nbytes_c, hipMemcpyDeviceToHost))

  max_error = 0.0d0
  do j = 0, Nc-1
    if (j == 0 .or. j == 1) then
      error = abs(hresult_c(j+1) - cmplx(dble(N), 0.0d0, kind=c_double_complex))
    else if (j == 3) then
      error = abs(hresult_c(j+1) - cmplx(1.5d0 * N, 0.0d0, kind=c_double_complex))
    else
      error = abs(hresult_c(j+1))
    end if
    max_error = max(max_error, error)
  end do

  if (max_error > tol) then
    write(*,*) "FAILED! R2C: max error = ", max_error
    call exit(1)
  end if

  ! Backward C2R: c2r(r2c(x)) should equal N * x
  call hipCheck(hipMemcpy(dy, c_loc(hresult_c(1)), Nbytes_c, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dz, Nbytes_r))

  plan = fftw_plan_dft_c2r_1d(N, dy, dz, FFTW_ESTIMATE)
  call fftw_execute_dft_c2r(plan, dy, dz)
  call fftw_destroy_plan(plan)

  call hipCheck(hipMemcpy(c_loc(hresult_r(1)), dz, Nbytes_r, hipMemcpyDeviceToHost))

  max_error = 0.0d0
  do j = 1, N
    error = abs(hresult_r(j) - dble(N) * hx(j))
    max_error = max(max_error, error)
  end do

  if (max_error > tol) then
    write(*,*) "FAILED! C2R round-trip: max error = ", max_error
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  call hipCheck(hipFree(dz))
  deallocate(hx, hresult_c, hresult_r)

  write(*,*) "PASSED!"
end program hipfftw_r2c_c2r_test

Multi-dimensional transforms#

fftw_plan_dft_2d takes the dimensions in C order, so the second argument varies fastest in memory.

program hipfftw_dft_2d_test
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftw

  implicit none

  integer(c_int), parameter :: NX = 4, NY = 8
  integer(c_int), parameter :: kx = 1, ky = 3
  integer(c_size_t), parameter :: Nbytes = NX * NY * 16  ! sizeof(double complex) = 16
  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-10
  ! FFTW_ESTIMATE is not emitted by the generated enums module; use the
  ! standard FFTW planner-flag value.
  integer(c_int), parameter :: FFTW_ESTIMATE = 64

  complex(c_double_complex), allocatable, target, dimension(:) :: hx, hresult
  type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr
  type(c_ptr) :: plan = c_null_ptr
  integer :: ii, jj, kxi, kyi
  double precision :: error, max_error
  complex(c_double_complex) :: wx, wy, expected

  write(*,"(a)",advance="no") "-- Running test 'hipfftw_dft_2d' (Fortran 2003) - "

  ! Signal: h[ii,jj] = exp(2*pi*i*(kx*ii/NX + ky*jj/NY)), C row-major (jj contiguous).
  ! Linear index: ii*NY + jj (0-indexed).
  ! Expected 2D DFT: output bin (kx,ky) = NX*NY, all others zero.
  wx = cmplx(0.0d0, 2.0d0 * pi * dble(kx) / dble(NX), kind=c_double_complex)
  wy = cmplx(0.0d0, 2.0d0 * pi * dble(ky) / dble(NY), kind=c_double_complex)
  allocate(hx(NX*NY), hresult(NX*NY))
  do ii = 0, NX-1
    do jj = 0, NY-1
      hx(ii*NY + jj + 1) = exp(wx*dble(ii)) * exp(wy*dble(jj))
    end do
  end do

  call hipCheck(hipMalloc(dx, Nbytes))
  call hipCheck(hipMalloc(dy, Nbytes))
  call hipCheck(hipMemcpy(dx, c_loc(hx(1)), Nbytes, hipMemcpyHostToDevice))

  ! n0=NX (slow), n1=NY (fast); C row-major, no dimension reversal needed.
  plan = fftw_plan_dft_2d(NX, NY, dx, dy, FFTW_FORWARD, FFTW_ESTIMATE)
  call fftw_execute_dft(plan, dx, dy)
  call fftw_destroy_plan(plan)

  call hipCheck(hipMemcpy(c_loc(hresult(1)), dy, Nbytes, hipMemcpyDeviceToHost))

  ! Verify: output bin (kx,ky) should equal NX*NY; all others zero.
  max_error = 0.0d0
  do kxi = 0, NX-1
    do kyi = 0, NY-1
      if (kxi == kx .and. kyi == ky) then
        expected = cmplx(dble(NX*NY), 0.0d0, kind=c_double_complex)
      else
        expected = cmplx(0.0d0, 0.0d0, kind=c_double_complex)
      end if
      error = abs(hresult(kxi*NY + kyi + 1) - expected)
      max_error = max(max_error, error)
    end do
  end do

  if (max_error > tol * dble(NX*NY)) then
    write(*,*) "FAILED! max error = ", max_error
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  deallocate(hx, hresult)

  write(*,*) "PASSED!"
end program hipfftw_dft_2d_test

test/f2003/hipfftw/hipfftw_dft_3d.f03 extends the same pattern to three dimensions with fftw_plan_dft_3d.

Batched transforms#

The _many planners transform a batch of signals with one plan. The inembed and onembed arguments describe the memory layout, stride is the gap between elements of one transform and dist the gap between the start of consecutive transforms. This example covers the complex-to-complex, real-to-complex and complex-to-real cases.

! =============================================================================
! GPU hipfftw test for the *_many interfaces
! =============================================================================
!
! This program validates hipfort's fftw_plan_many_dft, fftw_plan_many_dft_r2c,
! and fftw_plan_many_dft_c2r wrappers using GPU device memory.
!
! KEY CONVENTION: These are the generated bindings to the FFTW C API, which
! takes dimensions in C order (last index fastest, row-major) and does NOT
! reverse them. For a Fortran array z(NX, NY) with NX fastest, pass the
! reversed dimension list n = [NY, NX] (and likewise for the embed arrays).
! =============================================================================
program hipfftw_many_test
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftw
  implicit none

  double precision, parameter :: tol = 1.0d-12
  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  ! FFTW_ESTIMATE is not emitted by the generated enums module; use the
  ! standard FFTW planner-flag value.
  integer(c_int), parameter :: FFTW_ESTIMATE = 64
  integer :: nfail, seed_size
  integer, allocatable :: seed(:)

  call random_seed(size=seed_size)
  allocate(seed(seed_size))
  seed = 42
  call random_seed(put=seed)
  deallocate(seed)

  nfail = 0
  write(*,'(a)') "=== hipfftw many GPU Tests ==="

  call test_1d_c2c_interleaved(nfail)
  call test_2d_c2c_padded_embed(nfail)
  call test_1d_c2r_roundtrip(nfail)
  call test_many_vs_individual_c2c(nfail)

  write(*,*)
  if (nfail > 0) then
    write(*,'(i0,a)') nfail, " test(s) FAILED"
    call exit(1)
  else
    write(*,*) "PASSED!"
  end if

contains

  subroutine report(max_err, tolerance, nfail)
    double precision, intent(in) :: max_err, tolerance
    integer, intent(inout) :: nfail
    if (max_err > tolerance) then
      write(*,'(a,es12.4,a)') "FAILED (max_err=", max_err, ")"
      nfail = nfail + 1
    else
      write(*,'(a,es12.4,a)') "PASSED (max_err=", max_err, ")"
    end if
  end subroutine

  ! ===========================================================================
  ! 1D C2C, interleaved batches (batch-major layout)
  ! N=8, howmany=3, istride=howmany=3, idist=1
  ! Data layout: in(j*3+b+1) = signal(j) for transform b, element j
  ! Signal: x[j] = exp(2*pi*i*j/N) + 2*exp(2*pi*i*3*j/N)
  ! Expected DFT: X[1]=N, X[3]=2N, rest 0
  ! ===========================================================================
  subroutine test_1d_c2c_interleaved(nfail)
    integer, intent(inout) :: nfail
    integer(c_int), parameter :: N = 8, howmany = 3
    integer(c_size_t), parameter :: total_bytes = N * howmany * 16
    complex(c_double_complex), allocatable, target :: hx(:), hresult(:)
    type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr
    type(c_ptr) :: plan = c_null_ptr
    integer(c_int), target :: n_arr(1), ie(1), oe(1)
    integer :: j, b
    double precision :: max_err, err
    complex(c_double_complex) :: w, expected

    write(*,'(a)',advance="no") "  1D C2C interleaved batches:    "
    allocate(hx(N*howmany), hresult(N*howmany))
    hx = cmplx(0d0, 0d0, kind=c_double_complex)

    w = cmplx(0d0, 2d0*pi/dble(N), kind=c_double_complex)
    do b = 0, howmany-1
      do j = 0, N-1
        hx(j*howmany + b + 1) = exp(w*dble(j)) + 2d0*exp(3d0*w*dble(j))
      end do
    end do

    call hipCheck(hipMalloc(dx, total_bytes))
    call hipCheck(hipMalloc(dy, total_bytes))
    call hipCheck(hipMemcpy(dx, c_loc(hx(1)), total_bytes, hipMemcpyHostToDevice))

    n_arr = [N]; ie = [N]; oe = [N]
    plan = fftw_plan_many_dft(1, c_loc(n_arr), howmany, &
        dx, c_loc(ie), howmany, 1, dy, c_loc(oe), howmany, 1, FFTW_FORWARD, FFTW_ESTIMATE)
    call fftw_execute_dft(plan, dx, dy)
    call fftw_destroy_plan(plan)

    call hipCheck(hipMemcpy(c_loc(hresult(1)), dy, total_bytes, hipMemcpyDeviceToHost))

    max_err = 0d0
    do b = 0, howmany-1
      do j = 0, N-1
        if (j == 1) then
          expected = cmplx(dble(N), 0d0, kind=c_double_complex)
        else if (j == 3) then
          expected = cmplx(2d0*dble(N), 0d0, kind=c_double_complex)
        else
          expected = cmplx(0d0, 0d0, kind=c_double_complex)
        end if
        err = abs(hresult(j*howmany + b + 1) - expected)
        max_err = max(max_err, err)
      end do
    end do

    call report(max_err, tol, nfail)
    call hipCheck(hipFree(dx))
    call hipCheck(hipFree(dy))
    deallocate(hx, hresult)
  end subroutine

  ! ===========================================================================
  ! 2D C2C with padded embed (embed != n)
  ! NX=4 (fast in Fortran), NY=6 (slow), LDX=8 (padded), howmany=2
  ! Fortran array: z(LDX, NY, howmany), transform z(1:NX, 1:NY, :)
  !
  ! Generated bindings call the FFTW C API directly (no dimension reversal),
  ! so pass dimensions in C order: n=[NY,NX], inembed=[NY,LDX].
  !
  ! Signal: z(ix,iy) = exp(2*pi*i*(ix-1)/NX) * exp(2*pi*i*2*(iy-1)/NY)
  ! Expected 2D DFT (0-indexed): Z(kx=1,ky=2) = NX*NY, rest 0
  ! ===========================================================================
  subroutine test_2d_c2c_padded_embed(nfail)
    integer, intent(inout) :: nfail
    integer(c_int), parameter :: NX = 4, NY = 6, LDX = 8, howmany = 2
    integer(c_size_t), parameter :: total_bytes = LDX * NY * howmany * 16
    complex(c_double_complex), allocatable, target :: hx(:), hresult(:)
    type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr
    type(c_ptr) :: plan = c_null_ptr
    integer(c_int), target :: n_arr(2), ie(2), oe(2)
    integer :: ix, iy, kx, ky, b
    double precision :: max_err, err
    complex(c_double_complex) :: wx, wy, expected

    write(*,'(a)',advance="no") "  2D C2C padded embed:           "
    allocate(hx(LDX*NY*howmany), hresult(LDX*NY*howmany))
    hx = cmplx(0d0, 0d0, kind=c_double_complex)

    wx = cmplx(0d0, 2d0*pi/dble(NX), kind=c_double_complex)
    wy = cmplx(0d0, 2d0*pi/dble(NY), kind=c_double_complex)
    do b = 0, howmany-1
      do iy = 0, NY-1
        do ix = 0, NX-1
          hx(b*LDX*NY + iy*LDX + ix + 1) = exp(wx*dble(ix)) * exp(2d0*wy*dble(iy))
        end do
      end do
    end do

    call hipCheck(hipMalloc(dx, total_bytes))
    call hipCheck(hipMalloc(dy, total_bytes))
    call hipCheck(hipMemcpy(dx, c_loc(hx(1)), total_bytes, hipMemcpyHostToDevice))

    ! C order (reversed from Fortran): fastest dimension last.
    n_arr = [NY, NX]
    ie = [NY, LDX]
    oe = [NY, LDX]
    plan = fftw_plan_many_dft(2, c_loc(n_arr), howmany, &
        dx, c_loc(ie), 1, LDX*NY, dy, c_loc(oe), 1, LDX*NY, FFTW_FORWARD, FFTW_ESTIMATE)
    call fftw_execute_dft(plan, dx, dy)
    call fftw_destroy_plan(plan)

    call hipCheck(hipMemcpy(c_loc(hresult(1)), dy, total_bytes, hipMemcpyDeviceToHost))

    max_err = 0d0
    do b = 0, howmany-1
      do ky = 0, NY-1
        do kx = 0, NX-1
          if (kx == 1 .and. ky == 2) then
            expected = cmplx(dble(NX*NY), 0d0, kind=c_double_complex)
          else
            expected = cmplx(0d0, 0d0, kind=c_double_complex)
          end if
          err = abs(hresult(b*LDX*NY + ky*LDX + kx + 1) - expected)
          max_err = max(max_err, err)
        end do
      end do
    end do

    call report(max_err, tol, nfail)
    call hipCheck(hipFree(dx))
    call hipCheck(hipFree(dy))
    deallocate(hx, hresult)
  end subroutine

  ! ===========================================================================
  ! 1D C2R round-trip: c2r(r2c(x)) = N * x
  ! Tests both fftw_plan_many_dft_r2c and fftw_plan_many_dft_c2r.
  ! Uses random real input (any real signal is valid for R2C).
  ! ===========================================================================
  subroutine test_1d_c2r_roundtrip(nfail)
    integer, intent(inout) :: nfail
    integer(c_int), parameter :: N = 16, Nc = N/2+1, howmany = 3
    integer(c_size_t), parameter :: rbytes = N * howmany * 8
    integer(c_size_t), parameter :: cbytes = Nc * howmany * 16
    real(c_double), allocatable, target :: hx(:), hresult(:)
    type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr, dz = c_null_ptr
    type(c_ptr) :: plan = c_null_ptr
    integer(c_int), target :: n_arr(1), ie_r(1), oe_c(1), ie_c(1), oe_r(1)
    integer :: j, b
    double precision :: max_err, err

    write(*,'(a)',advance="no") "  1D R2C/C2R round-trip:         "
    allocate(hx(N*howmany), hresult(N*howmany))

    call random_number(hx)

    call hipCheck(hipMalloc(dx, rbytes))
    call hipCheck(hipMalloc(dy, cbytes))
    call hipCheck(hipMalloc(dz, rbytes))
    call hipCheck(hipMemcpy(dx, c_loc(hx(1)), rbytes, hipMemcpyHostToDevice))

    n_arr = [N]; ie_r = [N]; oe_c = [Nc]; ie_c = [Nc]; oe_r = [N]

    plan = fftw_plan_many_dft_r2c(1, c_loc(n_arr), howmany, &
        dx, c_loc(ie_r), 1, N, dy, c_loc(oe_c), 1, Nc, FFTW_ESTIMATE)
    call fftw_execute_dft_r2c(plan, dx, dy)
    call fftw_destroy_plan(plan)

    plan = fftw_plan_many_dft_c2r(1, c_loc(n_arr), howmany, &
        dy, c_loc(ie_c), 1, Nc, dz, c_loc(oe_r), 1, N, FFTW_ESTIMATE)
    call fftw_execute_dft_c2r(plan, dy, dz)
    call fftw_destroy_plan(plan)

    call hipCheck(hipMemcpy(c_loc(hresult(1)), dz, rbytes, hipMemcpyDeviceToHost))

    max_err = 0d0
    do b = 0, howmany-1
      do j = 0, N-1
        err = abs(hresult(b*N + j + 1) - dble(N) * hx(b*N + j + 1))
        max_err = max(max_err, err)
      end do
    end do

    call report(max_err, tol, nfail)
    call hipCheck(hipFree(dx))
    call hipCheck(hipFree(dy))
    call hipCheck(hipFree(dz))
    deallocate(hx, hresult)
  end subroutine

  ! ===========================================================================
  ! 1D C2C: many_dft vs individual dft_1d per batch on GPU
  ! Uses random complex input (any complex signal is valid for C2C).
  ! ===========================================================================
  subroutine test_many_vs_individual_c2c(nfail)
    integer, intent(inout) :: nfail
    integer(c_int), parameter :: N = 16, howmany = 5
    integer(c_size_t), parameter :: total_bytes = N * howmany * 16
    integer(c_size_t), parameter :: slice_bytes = N * 16
    complex(c_double_complex), allocatable, target :: hx(:), hout_many(:), hout_ind(:)
    type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr
    type(c_ptr) :: d_sin = c_null_ptr, d_sout = c_null_ptr
    type(c_ptr) :: plan = c_null_ptr
    integer(c_int), target :: n_arr(1), ie(1), oe(1)
    integer :: j
    double precision :: max_err, err
    real(c_double) :: rr(N*howmany), ri(N*howmany)
    integer :: b

    write(*,'(a)',advance="no") "  1D C2C many vs individual:     "
    allocate(hx(N*howmany), hout_many(N*howmany), hout_ind(N*howmany))

    call random_number(rr)
    call random_number(ri)
    do j = 1, N*howmany
      hx(j) = cmplx(rr(j), ri(j), kind=c_double_complex)
    end do

    call hipCheck(hipMalloc(dx, total_bytes))
    call hipCheck(hipMalloc(dy, total_bytes))
    call hipCheck(hipMemcpy(dx, c_loc(hx(1)), total_bytes, hipMemcpyHostToDevice))

    n_arr = [N]; ie = [N]; oe = [N]
    plan = fftw_plan_many_dft(1, c_loc(n_arr), howmany, &
        dx, c_loc(ie), 1, N, dy, c_loc(oe), 1, N, FFTW_FORWARD, FFTW_ESTIMATE)
    call fftw_execute_dft(plan, dx, dy)
    call fftw_destroy_plan(plan)
    call hipCheck(hipMemcpy(c_loc(hout_many(1)), dy, total_bytes, hipMemcpyDeviceToHost))
    call hipCheck(hipFree(dx))
    call hipCheck(hipFree(dy))

    call hipCheck(hipMalloc(d_sin, slice_bytes))
    call hipCheck(hipMalloc(d_sout, slice_bytes))

    do b = 0, howmany-1
      call hipCheck(hipMemcpy(d_sin, c_loc(hx(b*N+1)), slice_bytes, hipMemcpyHostToDevice))
      plan = fftw_plan_dft_1d(N, d_sin, d_sout, FFTW_FORWARD, FFTW_ESTIMATE)
      call fftw_execute_dft(plan, d_sin, d_sout)
      call fftw_destroy_plan(plan)
      call hipCheck(hipMemcpy(c_loc(hout_ind(b*N+1)), d_sout, slice_bytes, hipMemcpyDeviceToHost))
    end do
    call hipCheck(hipFree(d_sin))
    call hipCheck(hipFree(d_sout))

    max_err = 0d0
    do j = 1, N*howmany
      err = abs(hout_many(j) - hout_ind(j))
      max_err = max(max_err, err)
    end do

    call report(max_err, tol, nfail)
    deallocate(hx, hout_many, hout_ind)
  end subroutine

end program hipfftw_many_test

The guru interface#

The guru interface describes a transform as arrays of fftw_iodim descriptors, one per dimension, each giving a length and its input and output strides. It expresses layouts the simpler planners cannot.

Note how the arrays are passed: the dummy arguments are scalar type(fftw_iodim), so the example passes the first element of each array and the callee receives the base address of the contiguous struct array.

program hipfftw_guru_test
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftw
  use hipfort_hipfftw_types

  implicit none

  integer(c_int), parameter :: N = 16, HOWMANY = 3
  integer(c_size_t), parameter :: Nbytes = N * HOWMANY * 16  ! sizeof(double complex) = 16
  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-10
  ! FFTW_ESTIMATE is not emitted by the generated enums module; use the
  ! standard FFTW planner-flag value.
  integer(c_int), parameter :: FFTW_ESTIMATE = 64

  complex(c_double_complex), allocatable, target, dimension(:) :: hx, hresult
  type(c_ptr) :: dx = c_null_ptr, dy = c_null_ptr
  type(c_ptr) :: plan = c_null_ptr
  type(fftw_iodim), target :: dims(1), howmany_dims(1)
  integer :: b, j, k
  double precision :: error, max_error
  complex(c_double_complex) :: w, expected

  write(*,"(a)",advance="no") "-- Running test 'hipfftw_guru' (Fortran 2003) - "

  ! 1D batched C2C: HOWMANY batches of length N, contiguous (stride=1, dist=N).
  ! Batch b (0-indexed): x[j] = exp(2*pi*i*(b+1)*j/N).
  ! Expected forward DFT: output bin b+1 of batch b = N, all others zero.
  allocate(hx(N*HOWMANY), hresult(N*HOWMANY))
  do b = 0, HOWMANY-1
    w = cmplx(0.0d0, 2.0d0 * pi * dble(b+1) / dble(N), kind=c_double_complex)
    do j = 0, N-1
      hx(b*N + j + 1) = exp(w * dble(j))
    end do
  end do

  call hipCheck(hipMalloc(dx, Nbytes))
  call hipCheck(hipMalloc(dy, Nbytes))
  call hipCheck(hipMemcpy(dx, c_loc(hx(1)), Nbytes, hipMemcpyHostToDevice))

  dims(1) = fftw_iodim(N, 1, 1)         ! n=N, is=1, os=1
  howmany_dims(1) = fftw_iodim(HOWMANY, N, N)  ! n=HOWMANY, is=N, os=N
  ! The dims/howmany_dims dummies are declared scalar type(fftw_iodim); pass the
  ! first element of each array so the callee receives the base address of
  ! the contiguous struct array — looks like an element but acts as a pointer.
  plan = fftw_plan_guru_dft(1, dims(1), 1, howmany_dims(1), &
      dx, dy, FFTW_FORWARD, FFTW_ESTIMATE)
  call fftw_execute_dft(plan, dx, dy)
  call fftw_destroy_plan(plan)

  call hipCheck(hipMemcpy(c_loc(hresult(1)), dy, Nbytes, hipMemcpyDeviceToHost))

  ! Verify: for batch b, output bin b+1 should equal N; all others zero.
  max_error = 0.0d0
  do b = 0, HOWMANY-1
    do k = 0, N-1
      if (k == b+1) then
        expected = cmplx(dble(N), 0.0d0, kind=c_double_complex)
      else
        expected = cmplx(0.0d0, 0.0d0, kind=c_double_complex)
      end if
      error = abs(hresult(b*N + k + 1) - expected)
      max_error = max(max_error, error)
    end do
  end do

  if (max_error > tol * dble(N)) then
    write(*,*) "FAILED! max error = ", max_error
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  deallocate(hx, hresult)

  write(*,*) "PASSED!"
end program hipfftw_guru_test

Allocating buffers#

fftw_alloc_real and fftw_alloc_complex, along with the fftwf_ single precision forms, return correctly aligned host-accessible buffers that can be handed straight to a plan. Use c_f_pointer to get a Fortran array view, and release them with fftw_free.

! hipFFTW allocation API: fftw_alloc_real / fftw_alloc_complex and their
! single-precision counterparts return host-accessible buffers that can be
! handed straight to a plan, mirroring the allocation tests in
! clients/tests/hipfftw_test.cpp.
program hipfftw_alloc_test
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftw

  implicit none

  integer(c_int), parameter :: N = 16
  integer(c_int), parameter :: Nc = N / 2 + 1
  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-12
  real, parameter :: tol_s = 1.0e-4
  ! FFTW_ESTIMATE is not emitted by the generated enums module; use the
  ! standard FFTW planner-flag value.
  integer(c_int), parameter :: FFTW_ESTIMATE = 64

  type(c_ptr) :: pr = c_null_ptr, pc = c_null_ptr, plan = c_null_ptr
  real(c_double), pointer :: x(:) => null()
  complex(c_double_complex), pointer :: y(:) => null()
  real(c_float), pointer :: xs(:) => null()
  complex(c_float_complex), pointer :: ys(:) => null()
  integer :: j
  double precision :: error, max_error

  write(*,"(a)",advance="no") "-- Running test 'hipfftw_alloc' (Fortran 2003) - "

  ! Real signal: x[j] = 1 + 2*cos(2*pi*j/N) + 3*cos(2*pi*3*j/N)
  ! R2C output (0-indexed): X[0] = N, X[1] = N, X[3] = 3*N/2, rest zero.
  pr = fftw_alloc_real(int(N, c_size_t))
  pc = fftw_alloc_complex(int(Nc, c_size_t))
  if (.not. c_associated(pr) .or. .not. c_associated(pc)) then
    write(*,*) "FAILED! fftw_alloc_real/fftw_alloc_complex returned a null pointer"
    call exit(1)
  end if

  call c_f_pointer(pr, x, [N])
  call c_f_pointer(pc, y, [Nc])
  do j = 0, N-1
    x(j+1) = 1.0d0 + 2.0d0 * cos(2.0d0 * pi * j / dble(N)) &
                   + 3.0d0 * cos(6.0d0 * pi * j / dble(N))
  end do

  plan = fftw_plan_dft_r2c_1d(N, pr, pc, FFTW_ESTIMATE)
  call fftw_execute_dft_r2c(plan, pr, pc)
  call fftw_destroy_plan(plan)

  max_error = 0.0d0
  do j = 0, Nc-1
    if (j == 0 .or. j == 1) then
      error = abs(y(j+1) - cmplx(dble(N), 0.0d0, kind=c_double_complex))
    else if (j == 3) then
      error = abs(y(j+1) - cmplx(1.5d0 * N, 0.0d0, kind=c_double_complex))
    else
      error = abs(y(j+1))
    end if
    max_error = max(max_error, error)
  end do

  call fftw_free(pr)
  call fftw_free(pc)

  if (max_error > tol) then
    write(*,*) "FAILED! double precision: max error = ", max_error
    call exit(1)
  end if

  pr = fftwf_alloc_real(int(N, c_size_t))
  pc = fftwf_alloc_complex(int(Nc, c_size_t))
  if (.not. c_associated(pr) .or. .not. c_associated(pc)) then
    write(*,*) "FAILED! fftwf_alloc_real/fftwf_alloc_complex returned a null pointer"
    call exit(1)
  end if

  call c_f_pointer(pr, xs, [N])
  call c_f_pointer(pc, ys, [Nc])
  do j = 0, N-1
    xs(j+1) = real(1.0d0 + 2.0d0 * cos(2.0d0 * pi * j / dble(N)) &
                         + 3.0d0 * cos(6.0d0 * pi * j / dble(N)), c_float)
  end do

  plan = fftwf_plan_dft_r2c_1d(N, pr, pc, FFTW_ESTIMATE)
  call fftwf_execute_dft_r2c(plan, pr, pc)
  call fftwf_destroy_plan(plan)

  max_error = 0.0d0
  do j = 0, Nc-1
    if (j == 0 .or. j == 1) then
      error = abs(ys(j+1) - cmplx(real(N), 0.0, kind=c_float_complex))
    else if (j == 3) then
      error = abs(ys(j+1) - cmplx(1.5 * N, 0.0, kind=c_float_complex))
    else
      error = abs(ys(j+1))
    end if
    max_error = max(max_error, error)
  end do

  call fftwf_free(pr)
  call fftwf_free(pc)

  if (max_error > tol_s) then
    write(*,*) "FAILED! single precision: max error = ", max_error
    call exit(1)
  end if

  write(*,*) "PASSED!"
end program hipfftw_alloc_test
Contents