hipFFT examples

hipFFT examples#

hipFFT is a thin layer over rocFFT with an API that follows cuFFT. hipFORT exposes it through the hipfort_hipfft module, and its Xt and multi-process APIs through the hipfort_hipfftxt module, as CUDA Fortran does with cufft and cufftXt. hipfort_hipfftxt uses hipfort_hipfft, so use hipfort_hipfftxt alone gives the whole hipFFT API.

Every program on this page is complete and self-contained, and is built and run as part of the hipFORT test suite. The Fortran 2008 version of each program lives in test/f2008/hipfft, and the equivalent Fortran 2003 version, which uses type(c_ptr) device pointers and explicit byte counts instead of Fortran array pointers, lives in test/f2003/hipfft. The Xt programs exist only in the Fortran 2003 form, in test/f2003/hipfftxt.

If you want direct access to rocFFT rather than a cuFFT-style interface, see the rocFFT examples. For FFTW3-compatible code, see the hipFFTW examples.

Transform workflow#

A hipFFT transform follows the same sequence as cuFFT:

  1. Create a plan with hipfftPlan1d, hipfftPlan2d, hipfftPlan3d or hipfftPlanMany, passing the transform lengths and the transform type. hipfftPlan1d and hipfftPlanMany also take a batch count; hipfftPlan2d and hipfftPlan3d always describe a single transform.

  2. Run the transform with the hipfftExec routine matching the plan type: hipfftExecZ2Z and hipfftExecC2C for complex-to-complex, hipfftExecD2Z and hipfftExecR2C for real-to-complex, hipfftExecZ2D and hipfftExecC2R for complex-to-real.

  3. Release the plan with hipfftDestroy.

Keep the following conventions in mind:

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

  • The transform type encodes the precision: HIPFFT_Z2Z and HIPFFT_D2Z are double precision, HIPFFT_C2C and HIPFFT_R2C are single.

  • Complex-to-complex transforms take a direction, HIPFFT_FORWARD or HIPFFT_BACKWARD. Real transforms take their direction from the type.

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

  • Multi-dimensional plans take lengths in C order, with the last dimension contiguous. This is the opposite of rocFFT, which takes the fastest-varying dimension first.

  • Every hipFFT call returns a status code. The programs wrap them in hipfftCheck from the hipfort_check module, which aborts on failure.

Building and running#

The programs only need the hipfft and hip hipFORT components:

find_package(hipfort REQUIRED COMPONENTS hip hipfft)

add_executable(my_fft hipfft_c2c_1d_z.f08)
target_link_libraries(my_fft PRIVATE hipfort::hipfft hipfort::hip)

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

Complex-to-complex transform#

The simplest case: an in-place, single-batch, one-dimensional complex-to-complex transform in double precision. The program runs a forward transform followed by an inverse transform and checks that the result is N times the input, which demonstrates that hipFFT does not normalize.

program hipfft_c2c_1d_z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  integer(c_int), parameter :: N = 16

  complex(8), allocatable, target, dimension(:) :: hx, hx_input
  complex(8), pointer, dimension(:) :: dx => null()
  type(c_ptr) :: plan = c_null_ptr
  integer(c_int) :: direction
  integer :: i
  double precision :: error
  double precision, parameter :: error_max = 1.0d-8

  write(*,"(a)",advance="no") "-- Running test 'hipFFT C2C 1D double (z)' (Fortran 2008 interfaces) - "

  allocate(hx(N))
  allocate(hx_input(N))
  do i = 1, N
     hx(i) = cmplx(dble(i), dble(N - i), kind=8)
  end do
  hx_input(:) = hx(:)

  call hipCheck(hipMalloc(dx, source=hx))

  ! A single hipFFT plan serves both directions: unlike rocFFT, the direction is
  ! an argument of hipfftExecZ2Z rather than a property of the plan.
  call hipfftCheck(hipfftPlan1d(plan, N, HIPFFT_Z2Z, 1))

  direction = HIPFFT_FORWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  direction = HIPFFT_BACKWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  call hipfftCheck(hipfftDestroy(plan))

  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx))

  ! hipFFT is unnormalized, so forward+inverse yields N times the original input.
  do i = 1, N
     error = abs(hx(i) - N * hx_input(i))
     if (error > error_max * N) then
        write(*,*) "FAILED! i=", i, " error=", error, " got=", hx(i), " expected=", N * hx_input(i)
        STOP 1
     end if
  end do

  deallocate(hx)
  deallocate(hx_input)

  write(*,*) "PASSED!"

end program hipfft_c2c_1d_z

Use HIPFFT_C2C and complex(4) host data for a single-precision transform, as in test/f2008/hipfft/hipfft_c2c_1d_c.f08.

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

Real transforms use HIPFFT_D2Z and HIPFFT_Z2D. Because the spectrum of real data is Hermitian symmetric, the complex buffer holds N/2 + 1 elements.

program hipfft_r2c_c2r_1d_d
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  integer(c_int), parameter :: N = 16
  integer(c_int), parameter :: Ncomplex = N/2 + 1

  real(8), allocatable, target, dimension(:) :: hr, hr_input
  complex(8), allocatable, target, dimension(:) :: hc
  real(8), pointer, dimension(:) :: dr => null()
  complex(8), pointer, dimension(:) :: dc => null()
  type(c_ptr) :: plan_fwd = c_null_ptr
  type(c_ptr) :: plan_bwd = c_null_ptr
  integer :: i
  double precision :: error
  double precision, parameter :: error_max = 1.0d-8

  write(*,"(a)",advance="no") "-- Running test 'hipFFT R2C/C2R 1D double (d)' (Fortran 2008 interfaces) - "

  allocate(hr(N))
  allocate(hr_input(N))
  do i = 1, N
     hr(i) = dble(i) + dble(mod(i,3)) - dble(mod(i,7))
  end do
  hr_input(:) = hr(:)

  ! Device buffers: real input of length N, complex output of length N/2+1.
  ! The complex side holds N/2+1 elements due to Hermitian symmetry.
  call hipCheck(hipMalloc(dr, source=hr))
  allocate(hc(Ncomplex))
  call hipCheck(hipMalloc(dc, source=hc))

  ! Forward real-to-complex (out-of-place): dr -> dc.
  call hipfftCheck(hipfftPlan1d(plan_fwd, N, HIPFFT_D2Z, 1))
  call hipfftCheck(hipfftExecD2Z(plan_fwd, dr, dc))
  call hipCheck(hipDeviceSynchronize())
  call hipfftCheck(hipfftDestroy(plan_fwd))

  ! Inverse complex-to-real (out-of-place): dc -> dr. hipFFT unnormalized -> N*input.
  call hipfftCheck(hipfftPlan1d(plan_bwd, N, HIPFFT_Z2D, 1))
  call hipfftCheck(hipfftExecZ2D(plan_bwd, dc, dr))
  call hipCheck(hipDeviceSynchronize())
  call hipfftCheck(hipfftDestroy(plan_bwd))

  call hipCheck(hipMemcpy(hr, dr, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dr))
  call hipCheck(hipFree(dc))

  ! After forward+inverse the real data should equal N times the original input.
  do i = 1, N
     error = abs(hr(i) - N * hr_input(i))
     if (error > error_max * N) then
        write(*,*) "FAILED! i=", i, " error=", error, " got=", hr(i), " expected=", N * hr_input(i)
        STOP 1
     end if
  end do

  deallocate(hr)
  deallocate(hr_input)
  deallocate(hc)

  write(*,*) "PASSED!"

end program hipfft_r2c_c2r_1d_d

test/f2008/hipfft/hipfft_r2c_c2r_1d_s.f08 is the single-precision equivalent, using HIPFFT_R2C and HIPFFT_C2R.

Multi-dimensional transforms#

A two-dimensional transform uses hipfftPlan2d. The last dimension is contiguous, so a plan created as hipfftPlan2d(plan, Nx, Ny, type) expects Ny to vary fastest in memory.

program hipfft_c2c_2d_z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  integer(c_int), parameter :: Nx = 4, Ny = 8
  integer(c_int), parameter :: Ntot = Nx * Ny

  complex(8), allocatable, target, dimension(:) :: hx, hx_input
  complex(8), pointer, dimension(:) :: dx => null()
  type(c_ptr) :: plan = c_null_ptr
  integer(c_int) :: direction
  integer :: i
  double precision :: error
  double precision, parameter :: error_max = 1.0d-8

  write(*,"(a)",advance="no") "-- Running test 'hipFFT C2C 2D double (z)' (Fortran 2008 interfaces) - "

  allocate(hx(Ntot))
  allocate(hx_input(Ntot))
  do i = 1, Ntot
     hx(i) = cmplx(dble(i), dble(Ntot - i), kind=8)
  end do
  hx_input(:) = hx(:)

  call hipCheck(hipMalloc(dx, source=hx))

  ! In-place 2D C2C plan over an Nx*Ny complex double grid; the direction is
  ! supplied at exec time rather than at plan creation.
  call hipfftCheck(hipfftPlan2d(plan, Nx, Ny, HIPFFT_Z2Z))

  direction = HIPFFT_FORWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  direction = HIPFFT_BACKWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  call hipfftCheck(hipfftDestroy(plan))

  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx))

  ! hipFFT is unnormalized, so forward+inverse yields Ntot times the original input.
  do i = 1, Ntot
     error = abs(hx(i) - Ntot * hx_input(i))
     if (error > error_max * Ntot) then
        write(*,*) "FAILED! i=", i, " error=", error, " got=", hx(i), " expected=", Ntot * hx_input(i)
        STOP 1
     end if
  end do

  deallocate(hx)
  deallocate(hx_input)

  write(*,*) "PASSED!"

end program hipfft_c2c_2d_z

test/f2008/hipfft/hipfft_c2c_3d_z.f08 extends the same pattern to three dimensions with hipfftPlan3d.

Batched transforms#

To transform many signals with one plan, use hipfftPlanMany and pass the batch count. The inembed and onembed arrays describe the memory layout; passing null pointers selects the contiguous default.

program hipfft_c2c_1d_batched_z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  integer(c_int), parameter :: N = 16
  integer(c_int), parameter :: Nbatch = 4
  integer(c_int), parameter :: Ntot = N * Nbatch

  ! one_i avoids ambiguity: the raw hipfftPlanMany_ interface requires integer(c_int)
  ! VALUE arguments, and a bare literal '1' would have the default integer kind.
  integer(c_int), parameter :: one_i = 1

  complex(8), allocatable, target, dimension(:) :: hx, hx_input
  complex(8), pointer, dimension(:) :: dx => null()
  type(c_ptr) :: plan = c_null_ptr
  integer(c_int) :: direction
  ! nlen is the per-dimension transform length array passed by address to
  ! hipfftPlanMany (it cannot be called "n": Fortran is case-insensitive and N
  ! is already the transform length parameter).
  integer(c_int), target :: nlen(1)
  integer :: i
  double precision :: error
  double precision, parameter :: error_max = 1.0d-8

  write(*,"(a)",advance="no") "-- Running test 'hipFFT C2C 1D batched double (z)' (Fortran 2008 interfaces) - "

  nlen(1) = N

  allocate(hx(Ntot))
  allocate(hx_input(Ntot))
  do i = 1, Ntot
     hx(i) = cmplx(dble(i), dble(Ntot - i), kind=8)
  end do
  hx_input(:) = hx(:)

  call hipCheck(hipMalloc(dx, source=hx))

  ! NULL inembed/onembed selects the simple contiguous batched layout: each
  ! transform spans N consecutive complex doubles (istride=1, idist=N,
  ! ostride=1, odist=N). Passing c_loc(nlen) and c_null_ptr routes to the raw
  ! hipfftPlanMany_ C interface rather than the array
  ! overloads (which expect integer(c_int) arrays for n, inembed, onembed).
  ! Unlike rocFFT, a single plan handle serves both directions; the direction
  ! is supplied at exec time via hipfftExecZ2Z.
  call hipfftCheck(hipfftPlanMany(plan, one_i, c_loc(nlen), &
                                  c_null_ptr, one_i, N, &
                                  c_null_ptr, one_i, N, &
                                  HIPFFT_Z2Z, Nbatch))

  direction = HIPFFT_FORWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  direction = HIPFFT_BACKWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  call hipfftCheck(hipfftDestroy(plan))

  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx))

  ! hipFFT is unnormalized: forward+inverse yields N (per-transform length,
  ! not Ntot) times the original input for each element.
  do i = 1, Ntot
     error = abs(hx(i) - N * hx_input(i))
     if (error > error_max * N) then
        write(*,*) "FAILED! i=", i, " error=", error, " got=", hx(i), " expected=", N * hx_input(i)
        STOP 1
     end if
  end do

  deallocate(hx)
  deallocate(hx_input)

  write(*,*) "PASSED!"

end program hipfft_c2c_1d_batched_z

Advanced data layout#

hipfftPlanMany also describes strided and interleaved data. The stride is the distance between consecutive elements of one transform, and the distance is the gap between the start of consecutive transforms. This program batches a two-dimensional transform.

program hipfft_planmany_2d_z2z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  ! Advanced-interface 2D batched C2C transform (Fortran 2008 interfaces).
  ! Last dimension (Ny) is contiguous, matching C row-major order.
  integer(c_int), parameter :: Nx = 4, Ny = 5
  integer(c_int), parameter :: howmany = 3
  integer(c_int), parameter :: Ntot = Nx * Ny

  ! nlen, inembed_a, onembed_a: target arrays passed directly to the rank_1 wrapper.
  integer(c_int), target :: nlen(2)
  integer(c_int), target :: inembed_a(2)
  integer(c_int), target :: onembed_a(2)
  ! Scalar plan-many arguments: non-VALUE in the rank_1 wrapper, so use variables.
  integer(c_int) :: fft_rank, istride_v, idist_v, ostride_v, odist_v, howmany_v

  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-8

  complex(8), allocatable, target, dimension(:) :: hx
  complex(8), pointer, dimension(:) :: dx => null()
  type(c_ptr) :: plan = c_null_ptr
  integer(c_int) :: direction
  integer :: b, i, j, kx, ky, pos
  complex(8) :: expected
  double precision :: error, max_error

  write(*,"(a)",advance="no") &
    "-- Running test 'hipFFT PlanMany C2C 2D batched (z)' (Fortran 2008 interfaces) - "

  nlen      = [Nx, Ny]
  inembed_a = [Nx, Ny]
  onembed_a = [Nx, Ny]
  fft_rank  = 2
  istride_v = 1
  ostride_v = 1
  idist_v   = Nx * Ny
  odist_v   = Nx * Ny
  howmany_v = howmany

  allocate(hx(Ntot * howmany))

  ! Batch b holds exp(2*pi*i*(kx*ii/Nx + ky*jj/Ny)) with kx=b mod Nx, ky=(b+1) mod Ny.
  ! Forward transform must place Nx*Ny in bin (kx,ky) and zero everywhere else.
  ! A wrong stride, embed or batch distance moves energy into other bins and fails.
  do b = 0, howmany - 1
     kx = mod(b, Nx)
     ky = mod(b + 1, Ny)
     do i = 0, Nx - 1
        do j = 0, Ny - 1
           pos = b * idist_v + i * inembed_a(2) + j
           hx(pos + 1) = exp(cmplx(0.0d0, &
                2.0d0 * pi * (dble(kx * i) / dble(Nx) + dble(ky * j) / dble(Ny)), &
                kind=8))
        end do
     end do
  end do

  call hipCheck(hipMalloc(dx, source=hx))

  ! The rank_1 wrapper calls c_loc(nlen), c_loc(inembed_a), c_loc(onembed_a) internally;
  ! no c_loc here. Non-VALUE scalars must be variables (not literals).
  call hipfftCheck(hipfftPlanMany(plan, fft_rank, nlen, inembed_a, &
       istride_v, idist_v, onembed_a, ostride_v, odist_v, &
       HIPFFT_Z2Z, howmany_v))

  direction = HIPFFT_FORWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  call hipfftCheck(hipfftDestroy(plan))

  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx))

  max_error = 0.0d0
  do b = 0, howmany - 1
     kx = mod(b, Nx)
     ky = mod(b + 1, Ny)
     do i = 0, Nx - 1
        do j = 0, Ny - 1
           pos = b * odist_v + i * onembed_a(2) + j
           if (i == kx .and. j == ky) then
              expected = cmplx(dble(Ntot), 0.0d0, kind=8)
           else
              expected = (0.0d0, 0.0d0)
           end if
           error = abs(hx(pos + 1) - expected)
           max_error = max(max_error, error)
        end do
     end do
  end do

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

  deallocate(hx)

  write(*,*) "PASSED!"

end program hipfft_planmany_2d_z2z

Querying the work area size#

hipFFT needs scratch memory whose size depends on the transform. There are two ways to ask about it. hipfftEstimate1d and friends give a heuristic upper bound before a plan exists, which is useful for budgeting. hipfftGetSize1d and friends refine that estimate for a given handle and parameter set, and hipfftGetSize reports the exact requirement of a plan that has already been made.

program hipfft_estimate_getsize_d
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  ! hipfftEstimate1d/2d/3d vs hipfftGetSize1d/2d/3d vs hipfftGetSize.
  ! Estimate is a pre-plan heuristic upper bound; GetSize* is accurate.
  ! (Fortran 2008 interfaces)

  ! Sizes for the three dimension tests.
  integer(c_int), parameter :: N1 = 32, Nx2 = 8, Ny2 = 4
  integer(c_int), parameter :: Nx3 = 4, Ny3 = 4, Nz3 = 2

  ! D2Z real input / complex output for the functional transform check.
  integer(c_int), parameter :: Ncomplex = N1 / 2 + 1

  double precision, parameter :: tol = 1.0d-8

  real(8), allocatable, target, dimension(:) :: hr
  complex(8), allocatable, target, dimension(:) :: hc
  real(8), pointer, dimension(:) :: dr => null()
  complex(8), pointer, dimension(:) :: dc => null()
  type(c_ptr) :: plan = c_null_ptr
  integer(c_size_t) :: estSz, mkSz, gsSz
  integer :: i
  double precision :: dc_sum, dc_expected

  write(*,"(a)",advance="no") &
    "-- Running test 'hipFFT Estimate/GetSize double (d)' (Fortran 2008 interfaces) - "

  ! -----------------------------------------------------------------------
  ! 1-D: HIPFFT_D2Z, N=32, batch=1
  ! -----------------------------------------------------------------------
  call hipfftCheck(hipfftEstimate1d(N1, HIPFFT_D2Z, 1, estSz))

  call hipfftCheck(hipfftCreate(plan))
  ! hipfftGetSize needs a made plan; hipfftGetSize1d does not configure the handle.
  call hipfftCheck(hipfftMakePlan1d(plan, N1, HIPFFT_D2Z, 1, mkSz))
  call hipfftCheck(hipfftGetSize1d(plan, N1, HIPFFT_D2Z, 1, gsSz))
  call hipfftCheck(hipfftGetSize(plan, mkSz))
  call hipfftCheck(hipfftDestroy(plan))

  ! Estimate is an upper bound; strict equality is never guaranteed.
  if (estSz < gsSz) then
     write(*,*) "FAILED! 1D: Estimate < GetSize1d:", estSz, "<", gsSz
     STOP 1
  end if

  ! -----------------------------------------------------------------------
  ! 2-D: HIPFFT_Z2Z, Nx=8, Ny=4
  ! -----------------------------------------------------------------------
  call hipfftCheck(hipfftEstimate2d(Nx2, Ny2, HIPFFT_Z2Z, estSz))

  call hipfftCheck(hipfftCreate(plan))
  call hipfftCheck(hipfftMakePlan2d(plan, Nx2, Ny2, HIPFFT_Z2Z, mkSz))
  call hipfftCheck(hipfftGetSize2d(plan, Nx2, Ny2, HIPFFT_Z2Z, gsSz))
  call hipfftCheck(hipfftGetSize(plan, mkSz))
  call hipfftCheck(hipfftDestroy(plan))

  if (estSz < gsSz) then
     write(*,*) "FAILED! 2D: Estimate < GetSize2d:", estSz, "<", gsSz
     STOP 1
  end if

  ! -----------------------------------------------------------------------
  ! 3-D: HIPFFT_Z2Z, Nx=4, Ny=4, Nz=2
  ! -----------------------------------------------------------------------
  call hipfftCheck(hipfftEstimate3d(Nx3, Ny3, Nz3, HIPFFT_Z2Z, estSz))

  call hipfftCheck(hipfftCreate(plan))
  call hipfftCheck(hipfftMakePlan3d(plan, Nx3, Ny3, Nz3, HIPFFT_Z2Z, mkSz))
  call hipfftCheck(hipfftGetSize3d(plan, Nx3, Ny3, Nz3, HIPFFT_Z2Z, gsSz))
  call hipfftCheck(hipfftGetSize(plan, mkSz))
  call hipfftCheck(hipfftDestroy(plan))

  if (estSz < gsSz) then
     write(*,*) "FAILED! 3D: Estimate < GetSize3d:", estSz, "<", gsSz
     STOP 1
  end if

  ! -----------------------------------------------------------------------
  ! Functional check: 1-D D2Z, N=32. DC bin must equal sum(input).
  ! -----------------------------------------------------------------------
  allocate(hr(N1))
  allocate(hc(Ncomplex))

  do i = 1, N1
     hr(i) = dble(i) + dble(mod(i, 5)) * 0.5d0
  end do

  dc_expected = sum(hr)

  call hipCheck(hipMalloc(dr, source=hr))
  call hipCheck(hipMalloc(dc, source=hc))

  call hipfftCheck(hipfftPlan1d(plan, N1, HIPFFT_D2Z, 1))
  call hipfftCheck(hipfftExecD2Z(plan, dr, dc))
  call hipCheck(hipDeviceSynchronize())
  call hipfftCheck(hipfftDestroy(plan))

  call hipCheck(hipMemcpy(hc, dc, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dr))
  call hipCheck(hipFree(dc))

  ! DC bin (Fortran index 1) = sum of all input values.
  dc_sum = dble(real(hc(1), kind=8))
  if (abs(dc_sum - dc_expected) > tol * N1 * abs(dc_expected)) then
     write(*,*) "FAILED! DC bin =", dc_sum, " expected =", dc_expected
     STOP 1
  end if
  if (abs(aimag(hc(1))) > tol * N1 * abs(dc_expected)) then
     write(*,*) "FAILED! DC imaginary part non-zero:", aimag(hc(1))
     STOP 1
  end if

  deallocate(hr)
  deallocate(hc)

  write(*,*) "PASSED!"

end program hipfft_estimate_getsize_d

Managing the work area#

By default a plan allocates its own work area. Call hipfftSetAutoAllocation with 0 before building the plan to turn that off, then supply your own buffer with hipfftSetWorkArea. This lets several plans share one allocation, or lets the application control when the memory is reserved. Plans built this way use hipfftCreate and hipfftMakePlanMany rather than hipfftPlanMany, because auto-allocation has to be turned off between the two calls, and the required work area size is only known once hipfftMakePlanMany has returned it.

program hipfft_makeplanmany_z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  ! hipfftCreate / hipfftSetAutoAllocation / hipfftMakePlanMany /
  ! hipfftGetSizeMany / explicit work buffer / hipfftExecZ2Z
  ! (Fortran 2008 interfaces)
  integer(c_int), parameter :: N = 16
  integer(c_int), parameter :: batch_count = 3
  integer(c_int), parameter :: Ntot = N * batch_count

  ! nlen, inembed_1d, onembed_1d are target so the rank_1 wrapper can c_loc them.
  integer(c_int), target :: nlen(1), inembed_1d(1), onembed_1d(1)
  ! Non-VALUE scalars for the rank_1 wrappers.
  integer(c_int) :: fft_rank, istride_v, idist_v, ostride_v, odist_v, batch_v
  integer(kind(HIPFFT_Z2Z)) :: fft_type

  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-8

  complex(8), allocatable, target, dimension(:) :: hx
  complex(8), pointer, dimension(:) :: dx => null()
  type(c_ptr) :: plan = c_null_ptr
  type(c_ptr) :: workbuf = c_null_ptr
  integer(c_size_t) :: workSzMake, gSzMany
  integer(c_int) :: direction
  integer :: b, j, k_exp, pos
  complex(8) :: expected
  double precision :: error, max_error

  write(*,"(a)",advance="no") &
    "-- Running test 'hipFFT MakePlanMany Z2Z 1D batched (z)' (Fortran 2008 interfaces) - "

  nlen(1)      = N
  inembed_1d(1) = N
  onembed_1d(1) = N
  fft_rank  = 1
  istride_v = 1
  idist_v   = N
  ostride_v = 1
  odist_v   = N
  batch_v   = batch_count
  fft_type  = HIPFFT_Z2Z

  allocate(hx(Ntot))

  ! Batch b (0-indexed) holds exp(2*pi*i*(b+1)*j/N), j=0..N-1.
  ! Forward transform must place N in bin b+1 (0-indexed) and zero elsewhere.
  do b = 0, batch_count - 1
     do j = 0, N - 1
        hx(b * N + j + 1) = exp(cmplx(0.0d0, &
             2.0d0 * pi * dble((b + 1) * j) / dble(N), kind=8))
     end do
  end do

  call hipCheck(hipMalloc(dx, source=hx))

  call hipfftCheck(hipfftCreate(plan))
  call hipfftCheck(hipfftSetAutoAllocation(plan, 0_c_int))

  ! MakePlanMany configures the plan and returns the required work area size.
  call hipfftCheck(hipfftMakePlanMany(plan, fft_rank, nlen, inembed_1d, &
       istride_v, idist_v, onembed_1d, ostride_v, odist_v, &
       fft_type, batch_v, workSzMake))

  ! GetSizeMany queries the configured plan; must agree with MakePlanMany.
  call hipfftCheck(hipfftGetSizeMany(plan, fft_rank, nlen, inembed_1d, &
       istride_v, idist_v, onembed_1d, ostride_v, odist_v, &
       fft_type, batch_v, gSzMany))

  if (gSzMany /= workSzMake) then
     write(*,*) "FAILED! GetSizeMany=", gSzMany, " MakePlanMany=", workSzMake
     STOP 1
  end if

  ! Provide the work buffer explicitly (auto-alloc was disabled above).
  if (workSzMake > 0_c_size_t) then
     call hipCheck(hipMalloc(workbuf, workSzMake))
     call hipfftCheck(hipfftSetWorkArea(plan, workbuf))
  end if

  direction = HIPFFT_FORWARD
  call hipfftCheck(hipfftExecZ2Z(plan, dx, dx, direction))
  call hipCheck(hipDeviceSynchronize())

  call hipfftCheck(hipfftDestroy(plan))

  if (c_associated(workbuf)) call hipCheck(hipFree(workbuf))

  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx))

  max_error = 0.0d0
  do b = 0, batch_count - 1
     k_exp = b + 1
     do j = 0, N - 1
        pos = b * N + j
        if (j == k_exp) then
           expected = cmplx(dble(N), 0.0d0, kind=8)
        else
           expected = (0.0d0, 0.0d0)
        end if
        error = abs(hx(pos + 1) - expected)
        max_error = max(max_error, error)
     end do
  end do

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

  deallocate(hx)

  write(*,*) "PASSED!"

end program hipfft_makeplanmany_z

Running on HIP streams#

By default hipFFT executes on the null stream. Bind a plan to an application-owned stream with hipfftSetStream to overlap independent transforms. Each stream has to be synchronized before its results are read back. This program runs two plans, each on its own stream, with different input harmonics so that a swapped stream would be visible in the output.

program hipfft_setstream_z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfft

  implicit none

  ! Two 1-D Z2Z plans bound to different HIP streams, different input harmonics.
  ! A swapped stream would place energy in the wrong bin and fail verification.
  ! (Fortran 2008 interfaces)
  integer(c_int), parameter :: N = 16

  ! Distinct harmonics so a stream swap is observable in both buffers.
  integer, parameter :: k1 = 3, k2 = 7

  double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  double precision, parameter :: tol = 1.0d-8

  complex(8), allocatable, target, dimension(:) :: hx1, hx2
  complex(8), pointer, dimension(:) :: dx1 => null(), dx2 => null()
  type(c_ptr) :: plan1 = c_null_ptr, plan2 = c_null_ptr
  type(c_ptr) :: stream1 = c_null_ptr, stream2 = c_null_ptr
  integer(c_int) :: direction
  integer :: j
  double precision :: error, max_error

  write(*,"(a)",advance="no") &
    "-- Running test 'hipFFT SetStream Z2Z 1D (z)' (Fortran 2008 interfaces) - "

  allocate(hx1(N))
  allocate(hx2(N))

  ! Buffer 1: pure harmonic at bin k1 (0-indexed).
  do j = 0, N - 1
     hx1(j + 1) = exp(cmplx(0.0d0, 2.0d0 * pi * dble(k1 * j) / dble(N), kind=8))
  end do

  ! Buffer 2: pure harmonic at bin k2 (0-indexed).
  do j = 0, N - 1
     hx2(j + 1) = exp(cmplx(0.0d0, 2.0d0 * pi * dble(k2 * j) / dble(N), kind=8))
  end do

  call hipCheck(hipMalloc(dx1, source=hx1))
  call hipCheck(hipMalloc(dx2, source=hx2))

  call hipCheck(hipStreamCreate(stream1))
  call hipCheck(hipStreamCreate(stream2))

  call hipfftCheck(hipfftPlan1d(plan1, N, HIPFFT_Z2Z, 1))
  call hipfftCheck(hipfftPlan1d(plan2, N, HIPFFT_Z2Z, 1))

  call hipfftCheck(hipfftSetStream(plan1, stream1))
  call hipfftCheck(hipfftSetStream(plan2, stream2))

  direction = HIPFFT_FORWARD
  call hipfftCheck(hipfftExecZ2Z(plan1, dx1, dx1, direction))
  call hipfftCheck(hipfftExecZ2Z(plan2, dx2, dx2, direction))

  ! Synchronize each stream independently (not hipDeviceSynchronize).
  call hipCheck(hipStreamSynchronize(stream1))
  call hipCheck(hipStreamSynchronize(stream2))

  call hipfftCheck(hipfftDestroy(plan1))
  call hipfftCheck(hipfftDestroy(plan2))

  call hipCheck(hipMemcpy(hx1, dx1, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hx2, dx2, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx1))
  call hipCheck(hipFree(dx2))

  call hipCheck(hipStreamDestroy(stream1))
  call hipCheck(hipStreamDestroy(stream2))

  ! Verify buffer 1: energy only at Fortran index k1+1 (0-indexed bin k1).
  max_error = 0.0d0
  do j = 1, N
     if (j == k1 + 1) then
        error = abs(hx1(j) - cmplx(dble(N), 0.0d0, kind=8))
     else
        error = abs(hx1(j))
     end if
     max_error = max(max_error, error)
  end do

  if (max_error > tol * N) then
     write(*,*) "FAILED! buffer1 max error = ", max_error
     STOP 1
  end if

  ! Verify buffer 2: energy only at Fortran index k2+1 (0-indexed bin k2).
  max_error = 0.0d0
  do j = 1, N
     if (j == k2 + 1) then
        error = abs(hx2(j) - cmplx(dble(N), 0.0d0, kind=8))
     else
        error = abs(hx2(j))
     end if
     max_error = max(max_error, error)
  end do

  if (max_error > tol * N) then
     write(*,*) "FAILED! buffer2 max error = ", max_error
     STOP 1
  end if

  deallocate(hx1)
  deallocate(hx2)

  write(*,*) "PASSED!"

end program hipfft_setstream_z

Transforms on hipDataType with the Xt API#

hipfftXtMakePlanMany describes a batched transform like hipfftMakePlanMany, but takes 64-bit lengths, strides and distances, and names the input, output and execution data types with hipDataType values (HIP_C_32F, HIP_C_64F, HIP_R_32F, …) instead of a transform type. hipfftXtGetSizeMany returns the work area size for the same parameters, and hipfftXtExec runs the plan whatever its data types, taking the direction for complex transforms.

program hipfftxt_makeplanmany_c
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftxt

  implicit none

  integer(c_int), parameter :: N = 32, batch = 3
  integer(c_size_t), parameter :: Nbytes = int(N,c_size_t) * int(batch,c_size_t) * 8
  real, parameter :: error_max = 1.0e-5

  complex(c_float_complex), target :: hx(N, batch), hx_input(N, batch)
  integer(c_int64_t), target :: nn(1)
  type(c_ptr) :: plan = c_null_ptr
  type(c_ptr) :: dx = c_null_ptr
  integer(c_size_t) :: work, work_query
  integer :: i, b
  real :: error

  write(*,"(a)",advance="no") "-- Running test 'hipFFT XtMakePlanMany C2C' (Fortran 2003 interfaces) - "

  do b = 1, batch
     do i = 1, N
        hx(i, b) = cmplx(real(i + b), real(N - i), c_float_complex)
     end do
  end do
  hx_input = hx

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

  ! The Xt plan names its data types (hipDataType) instead of an hipfftType.
  nn(1) = N
  call hipfftCheck(hipfftCreate(plan))
  call hipfftCheck(hipfftXtGetSizeMany(plan, 1, c_loc(nn), c_null_ptr, 1_c_int64_t, &
       int(N,c_int64_t), HIP_C_32F, c_null_ptr, 1_c_int64_t, int(N,c_int64_t), HIP_C_32F, &
       int(batch,c_int64_t), work_query, HIP_C_32F))
  call hipfftCheck(hipfftXtMakePlanMany(plan, 1, c_loc(nn), c_null_ptr, 1_c_int64_t, &
       int(N,c_int64_t), HIP_C_32F, c_null_ptr, 1_c_int64_t, int(N,c_int64_t), HIP_C_32F, &
       int(batch,c_int64_t), work, HIP_C_32F))

  call hipfftCheck(hipfftXtExec(plan, dx, dx, HIPFFT_FORWARD))
  call hipfftCheck(hipfftXtExec(plan, dx, dx, HIPFFT_BACKWARD))
  call hipCheck(hipDeviceSynchronize())
  call hipfftCheck(hipfftDestroy(plan))

  call hipCheck(hipMemcpy(c_loc(hx), dx, Nbytes, hipMemcpyDeviceToHost))
  call hipCheck(hipFree(dx))

  ! hipFFT is unnormalized, so forward+inverse yields N times the original input.
  do b = 1, batch
     do i = 1, N
        error = abs(hx(i, b) - N * hx_input(i, b)) / (N * max(1.0, abs(hx_input(i, b))))
        if (.not. (error <= error_max)) then
           write(*,*) "FAILED! i=", i, " b=", b, " error=", error
           call exit(1)
        end if
     end do
  end do

  write(*,*) "PASSED!"

end program hipfftxt_makeplanmany_c

Multi-GPU transforms#

A plan can be spread over several GPUs of the same process. The data then lives in a descriptor that hipFFT allocates and distributes, rather than in a buffer of yours:

  1. Create an empty plan with hipfftCreate, and list the GPUs with hipfftXtSetGPUs before making it.

  2. Make the plan with hipfftMakePlan2d, hipfftMakePlan3d or hipfftMakePlanMany. A multi-GPU plan writes one work area size per GPU, so pass c_loc of an integer(c_size_t), target array with one slot per GPU, through the type(c_ptr) form. The form that takes an integer(c_size_t) variable is for a plan on one GPU.

  3. Allocate the data with hipfftXtMalloc. It returns the address of a hipLibXtDesc in a type(c_ptr).

  4. Copy the host data in and out with hipfftXtMemcpy, passing that type(c_ptr) on the device side, and run the transform with the hipfftXtExecDescriptor routine matching the plan type (hipfftXtExecDescriptorZ2Z, hipfftXtExecDescriptorC2C, …).

  5. Release the descriptor with hipfftXtFree, then the plan with hipfftDestroy.

The descriptor is owned by hipFFT. Associate a type(hipLibXtDesc), pointer with the address hipfftXtMalloc returned, using c_f_pointer, and pass that pointer to hipfftXtExecDescriptor* and hipfftXtFree. Never pass a copy of the type: hipfftXtFree would then release the wrong address. The descriptor field leads in the same way to a hipXtDesc, whose nGPUs, GPUs, data and size fields say where each part of the data is.

The hipLibXtDesc and hipXtDesc types are in the hipfort_hipfftxt_types module, and the Xt constants (HIPFFT_XT_FORMAT_INPLACE, HIPFFT_COPY_HOST_TO_DEVICE, …) in hipfort_hipfftxt_enums; hipfort_hipfftxt uses both.

On AMD GPUs a descriptor needs a plan over two or more GPUs, and unbatched one-dimensional transforms are not supported. This program runs a two-dimensional transform on two GPUs, forward then inverse, and checks the result against the input; on a machine with fewer than two GPUs it reports so and the test is skipped.

program hipfftxt_descriptor_z
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipfftxt

  implicit none

  integer(c_int), parameter :: NX = 16, NY = 16
  double precision, parameter :: error_max = 1.0d-10

  complex(c_double_complex), target :: hx(NX, NY), hx_input(NX, NY)
  integer(c_int), target :: gpus(2)
  type(c_ptr) :: plan = c_null_ptr
  type(c_ptr) :: d = c_null_ptr
  type(hipLibXtDesc), pointer :: desc
  type(hipXtDesc), pointer :: xt
  integer(c_size_t), target :: work(2)   ! one work size per GPU
  integer(c_int) :: ndev
  integer :: i, j
  double precision :: error

  write(*,"(a)",advance="no") "-- Running test 'hipFFT Xt descriptor 2D Z2Z' (Fortran 2003 interfaces) - "

  ! The descriptor API needs a plan spread over at least two GPUs.
  call hipCheck(hipGetDeviceCount(ndev))
  ! Exit 77, which CTest reports as skipped: nothing was checked.
  if (ndev < 2) then
     write(*,*) "SKIPPED (needs two GPUs, found ", ndev, ")"
     stop 77
  end if

  do j = 1, NY
     do i = 1, NX
        hx(i, j) = cmplx(dble(i + j), dble(NX - i), c_double_complex)
     end do
  end do
  hx_input = hx

  gpus = [0, 1]
  call hipfftCheck(hipfftCreate(plan))
  call hipfftCheck(hipfftXtSetGPUs(plan, 2, c_loc(gpus)))
  ! A multi-GPU plan writes one work size per GPU: pass c_loc of an array with
  ! one slot per GPU to the type(c_ptr) form.
  call hipfftCheck(hipfftMakePlan2d(plan, NY, NX, HIPFFT_Z2Z, c_loc(work)))

  ! hipFFT allocates and owns the descriptor: keep its address, and pass the
  ! pointer that c_f_pointer associates with it, never a copy of the type.
  call hipfftCheck(hipfftXtMalloc(plan, d, HIPFFT_XT_FORMAT_INPLACE))
  call c_f_pointer(d, desc)
  call c_f_pointer(desc%descriptor, xt)
  if (xt%nGPUs /= 2) then
     write(*,*) "FAILED! descriptor spans ", xt%nGPUs, " GPUs (expected 2)"
     call exit(1)
  end if

  call hipfftCheck(hipfftXtMemcpy(plan, d, c_loc(hx), HIPFFT_COPY_HOST_TO_DEVICE))
  call hipfftCheck(hipfftXtExecDescriptorZ2Z(plan, desc, desc, HIPFFT_FORWARD))
  call hipfftCheck(hipfftXtExecDescriptorZ2Z(plan, desc, desc, HIPFFT_BACKWARD))
  call hipfftCheck(hipfftXtMemcpy(plan, c_loc(hx), d, HIPFFT_COPY_DEVICE_TO_HOST))
  call hipfftCheck(hipfftXtFree(desc))
  call hipfftCheck(hipfftDestroy(plan))

  ! hipFFT is unnormalized, so forward+inverse yields NX*NY times the input.
  do j = 1, NY
     do i = 1, NX
        error = abs(hx(i, j) - NX * NY * hx_input(i, j)) / (NX * NY * max(1.0d0, abs(hx_input(i, j))))
        if (.not. (error <= error_max)) then
           write(*,*) "FAILED! i=", i, " j=", j, " error=", error
           call exit(1)
        end if
     end do
  end do

  write(*,*) "PASSED!"

end program hipfftxt_descriptor_z

Multi-process transforms with MPI#

hipfftMpAttachComm, hipfftXtSetDistribution and hipfftXtSetSubformatDefault distribute a transform over the ranks of an MPI communicator. They are bound in hipfort_hipfftxt, but exist only in a hipFFT built with MPI support (HIPFFT_MPI_ENABLE), which the ROCm packages are not: with those, a program that calls them does not link.

hipfftMpAttachComm takes a type(c_ptr) to a C MPI_Comm, not a Fortran communicator handle. The size and the representation of MPI_Comm depend on the MPI implementation, so build it in C from the Fortran handle with MPI_Comm_f2c, for example with this helper:

#include <mpi.h>
#include <stdlib.h>

/* The C communicator for a Fortran one, in storage the caller frees. */
void* c_comm_from_fortran(MPI_Fint fcomm)
{
  MPI_Comm* comm = malloc(sizeof *comm);
  *comm = MPI_Comm_f2c(fcomm);
  return comm;
}

use iso_c_binding
use mpi_f08
use hipfort_check
use hipfort_hipfftxt

interface
  function c_comm_from_fortran(fcomm) bind(c, name="c_comm_from_fortran")
    import :: c_ptr, c_int
    integer(c_int), value :: fcomm
    type(c_ptr) :: c_comm_from_fortran
  end function
end interface

type(c_ptr) :: comm

comm = c_comm_from_fortran(MPI_COMM_WORLD%MPI_VAL)
call hipfftCheck(hipfftMpAttachComm(plan, HIPFFT_COMM_MPI, comm))

hipFFT keeps the pointer rather than a copy of the communicator, so the storage it points to has to outlive the plan. The plan must be attached to the communicator before it is made.

Contents