finufft.F90

Go to the documentation of this file.

!! Copyright (C) 2026 Octopus developers
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
 
#include "global.h"
 
module finufft_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_c_binding
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  implicit none
 
  private
 
  public ::                  &
    finufft_t,               &
    finufft_init,            &
    finufft_end,             &
    finufft_precompute,      &
    finufft_set_nodes,       &
    finufft_modes_to_nodes,  &
    finufft_nodes_to_modes,  &
    finufft_write_info,      &
    finufft_read_options,    &
    finufft_set_tolerance,   &
    zfinufft_forward,     &
    zfinufft_backward,    &
    dfinufft_forward,     &
    dfinufft_backward
 
#ifdef HAVE_NFFT
 
  integer, parameter :: FINUFFT_TYPE2 = 2
  integer, parameter :: FINUFFT_IFLAG = -1
#endif
 
  type finufft_t
    private
 
    integer               :: N(3)
    integer               :: M(3)
    integer               :: dim
    real(real64), public  :: norm
 
    real(real64)          :: tol
 
    type(c_ptr)           :: plan = c_null_ptr  
 
    real(real64), pointer :: xj(:) => null()
    real(real64), pointer :: yj(:) => null()
    real(real64), pointer :: zj(:) => null()
 
  end type finufft_t
 
#ifdef HAVE_NFFT
 
  interface
 
    type(c_ptr) function oct_finufft_opts_init(nthreads) bind(c, name='oct_finufft_opts_init')
      import :: c_ptr, c_int
      integer(c_int), value :: nthreads
    end function oct_finufft_opts_init
 
    subroutine oct_finufft_opts_end(opts) bind(c, name='oct_finufft_opts_end')
      import :: c_ptr
      type(c_ptr), value :: opts
    end subroutine oct_finufft_opts_end
 
    integer(c_int) function finufft_makeplan(type_, dim, n_modes, iflag, n_transf, tol, plan, opts) &
      bind(c, name='finufft_makeplan')
      import :: c_int, c_int64_t, c_double, c_ptr
      integer(c_int),     value        :: type_
      integer(c_int),     value        :: dim
      integer(c_int64_t), intent(in)   :: n_modes(*)
      integer(c_int),     value        :: iflag
      integer(c_int),     value        :: n_transf
      real(c_double),     value        :: tol
      type(c_ptr),        intent(out)  :: plan
      type(c_ptr),        value        :: opts
    end function finufft_makeplan
 
    integer(c_int) function finufft_setpts(plan, M, xj, yj, zj, N, s, t, u) &
      bind(c, name='finufft_setpts')
      import :: c_int, c_int64_t, c_ptr
      type(c_ptr),        value :: plan
      integer(c_int64_t), value :: M
      type(c_ptr),        value :: xj, yj, zj
      integer(c_int64_t), value :: N
      type(c_ptr),        value :: s, t, u
    end function finufft_setpts
 
    integer(c_int) function finufft_execute(plan, cj, fk) bind(c, name='finufft_execute')
      import :: c_int, c_ptr, c_double_complex
      type(c_ptr),               value         :: plan
      complex(c_double_complex), intent(inout) :: cj(*)
      complex(c_double_complex), intent(in)    :: fk(*)
    end function finufft_execute
 
    integer(c_int) function finufft_execute_adjoint(plan, cj, fk) bind(c, name='finufft_execute_adjoint')
      import :: c_int, c_ptr, c_double_complex
      type(c_ptr),               value         :: plan
      complex(c_double_complex), intent(in)    :: cj(*)
      complex(c_double_complex), intent(inout) :: fk(*)
    end function finufft_execute_adjoint
 
    integer(c_int) function finufft_destroy(plan) bind(c, name='finufft_destroy')
      import :: c_int, c_ptr
      type(c_ptr), value :: plan
    end function finufft_destroy
 
  end interface
#endif
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine finufft_read_options(finufft, namespace)
    type(finufft_t),   intent(inout) :: finufft
    type(namespace_t), intent(in)    :: namespace
 
    push_sub(finufft_read_options)
 
    !%Variable FINUFFTTolerance
    !%Type float
    !%Default 1e-8
    !%Section Mesh::FFTs
    !%Description
    !% Requested relative tolerance of the non-equispaced FFTs. Tighter
    !% tolerances widen the spreading kernel and are more expensive.
    !%End
    call parse_variable(namespace, 'FINUFFTTolerance', 1.0e-8_real64, finufft%tol)
 
    if (finufft%tol <= m_zero) then
      call messages_input_error(namespace, 'FINUFFTTolerance')
    end if
 
    pop_sub(finufft_read_options)
  end subroutine finufft_read_options
 
  ! ---------------------------------------------------------
  subroutine finufft_set_tolerance(finufft, tol)
    type(finufft_t), intent(inout) :: finufft
    real(real64),    intent(in)    :: tol
 
    push_sub(finufft_set_tolerance)
 
    finufft%tol = tol
 
    pop_sub(finufft_set_tolerance)
  end subroutine finufft_set_tolerance
 
  ! ---------------------------------------------------------
  subroutine finufft_init(finufft, finufft_options, N, dim, M, nthreads)
    type(finufft_t), intent(inout) :: finufft
    type(finufft_t), intent(in)    :: finufft_options
    integer,         intent(in)    :: n(3)
    integer,         intent(in)    :: dim
    integer,         intent(in)    :: m(3)
    integer, intent(in), optional  :: nthreads
 
#ifdef HAVE_NFFT
    integer(c_int64_t) :: n_modes(3)
    integer(c_int)     :: ier
    type(c_ptr)        :: finufft_opts
    integer(c_int)     :: n_threads
#endif
 
    push_sub(finufft_init)
 
    call profiling_in("FINUFFT_INIT")
 
    finufft%dim = dim
    finufft%N(:) = n(:)
    finufft%M(:) = m(:)
    ! unused directions hold a single mode and a single node
    finufft%N(dim+1:3) = 1
    finufft%M(dim+1:3) = 1
 
    finufft%tol = finufft_options%tol
 
#ifdef HAVE_NFFT
    n_modes(1:3) = int(finufft%N(1:3), c_int64_t)
    n_threads = optional_default(nthreads, 0_c_int)
 
    ! Selects the FINUFFT defaults, including the mode ordering relied on here
    ! (modeord = 0, modes -N/2 .. N/2-1)
    finufft_opts = oct_finufft_opts_init(n_threads)
    if (.not. c_associated(finufft_opts)) then
      message(1) = "Could not allocate the FINUFFT options."
      call messages_fatal(1)
    end if
 
    ier = finufft_makeplan(int(finufft_type2, c_int), int(dim, c_int), n_modes, &
      int(finufft_iflag, c_int), 1_c_int, real(finufft%tol, c_double), finufft%plan, finufft_opts)
    call finufft_check(ier, 'finufft_makeplan')
 
    call oct_finufft_opts_end(finufft_opts)
#endif
 
    call profiling_out("FINUFFT_INIT")
 
    pop_sub(finufft_init)
  end subroutine finufft_init
 
  ! ---------------------------------------------------------
  subroutine finufft_precompute(finufft, X1, X2, X3, period)
    type(finufft_t), intent(inout) :: finufft
    real(real64), contiguous, intent(in) :: x1(:)
    real(real64), contiguous, optional, intent(in) :: x2(:)
    real(real64), contiguous, optional, intent(in) :: x3(:)
    real(real64), optional, intent(in) :: period(:)
 
    real(real64)       :: v1(1:finufft%m(1)), v2(1:finufft%m(2)), v3(1:finufft%m(3))
    real(real64)       :: inv_spacing, period_(3)
    integer(c_int64_t) :: nj
    integer            :: i1, i2, i3, jj
#ifdef HAVE_NFFT
    integer(c_int)     :: ier
    type(c_ptr)        :: yj_ptr, zj_ptr
#endif
 
    push_sub(finufft_precompute)
 
    call profiling_in("FINUFFT_PRECOMPUTE")
 
    finufft%norm = m_one
 
    period_ = m_zero
    if (present(period)) period_(1:3) = period(1:3)
 
    assert(size(x1) == finufft%M(1))
    call finufft_scale_nodes(x1, v1, inv_spacing, period_(1))
    finufft%norm = finufft%norm * inv_spacing
 
    v2 = m_zero
    if (finufft%dim >= 2) then
      assert(present(x2))
      assert(size(x2) == finufft%M(2))
      call finufft_scale_nodes(x2, v2, inv_spacing, period_(2))
      finufft%norm = finufft%norm * inv_spacing
    end if
 
    v3 = m_zero
    if (finufft%dim >= 3) then
      assert(present(x3))
      assert(size(x3) == finufft%M(3))
      call finufft_scale_nodes(x3, v3, inv_spacing, period_(3))
      finufft%norm = finufft%norm * inv_spacing
    end if
 
    nj = int(finufft%M(1), c_int64_t)*int(finufft%M(2), c_int64_t)*int(finufft%M(3), c_int64_t)
 
    safe_deallocate_p(finufft%xj)
    safe_deallocate_p(finufft%yj)
    safe_deallocate_p(finufft%zj)
    safe_allocate(finufft%xj(1:nj))
    if (finufft%dim >= 2) then
      safe_allocate(finufft%yj(1:nj))
    end if
    if (finufft%dim >= 3) then
      safe_allocate(finufft%zj(1:nj))
    end if
 
    ! Fortran order, matching the transform arrays, so execute needs no repacking
    jj = 0
    do i3 = 1, finufft%M(3)
      do i2 = 1, finufft%M(2)
        do i1 = 1, finufft%M(1)
          jj = jj + 1
          finufft%xj(jj) = m_two*m_pi*v1(i1)
          if (finufft%dim >= 2) finufft%yj(jj) = m_two*m_pi*v2(i2)
          if (finufft%dim >= 3) finufft%zj(jj) = m_two*m_pi*v3(i3)
        end do
      end do
    end do
 
#ifdef HAVE_NFFT
    yj_ptr = c_null_ptr
    zj_ptr = c_null_ptr
    if (finufft%dim >= 2) yj_ptr = c_loc(finufft%yj)
    if (finufft%dim >= 3) zj_ptr = c_loc(finufft%zj)
 
    ! The trailing arguments only apply to type 3
    ier = finufft_setpts(finufft%plan, nj, c_loc(finufft%xj), yj_ptr, zj_ptr, &
      0_c_int64_t, c_null_ptr, c_null_ptr, c_null_ptr)
    call finufft_check(ier, 'finufft_setpts')
#endif
 
    call profiling_out("FINUFFT_PRECOMPUTE")
 
    write(message(1), '(a)') "Info: FINUFFT plan precomputed."
    call messages_info(1)
 
    pop_sub(finufft_precompute)
  end subroutine finufft_precompute
 
  ! ---------------------------------------------------------
  subroutine finufft_set_nodes(finufft, xj, yj, zj, period)
    type(finufft_t), intent(inout) :: finufft
    real(real64), contiguous, intent(in) :: xj(:)
    real(real64), contiguous, intent(in) :: yj(:)
    real(real64), contiguous, intent(in) :: zj(:)
    real(real64), contiguous, intent(in) :: period(:)
 
    integer(c_int64_t) :: nj
#ifdef HAVE_NFFT
    integer(c_int)     :: ier
#endif
 
    push_sub(finufft_set_nodes)
 
    call profiling_in("FINUFFT_SET_NODES")
 
    assert(finufft%dim == 3)
 
    nj = int(finufft%M(1), c_int64_t)
    assert(size(xj) == nj)
    assert(size(yj) == nj)
    assert(size(zj) == nj)
 
    finufft%norm = m_one
 
    safe_deallocate_p(finufft%xj)
    safe_deallocate_p(finufft%yj)
    safe_deallocate_p(finufft%zj)
    safe_allocate(finufft%xj(1:nj))
    safe_allocate(finufft%yj(1:nj))
    safe_allocate(finufft%zj(1:nj))
 
    finufft%xj = m_two*m_pi*xj/period(1)
    finufft%yj = m_two*m_pi*yj/period(2)
    finufft%zj = m_two*m_pi*zj/period(3)
 
#ifdef HAVE_NFFT
    ier = finufft_setpts(finufft%plan, nj, c_loc(finufft%xj), c_loc(finufft%yj), c_loc(finufft%zj), &
      0_c_int64_t, c_null_ptr, c_null_ptr, c_null_ptr)
    call finufft_check(ier, 'finufft_setpts')
#endif
 
    call profiling_out("FINUFFT_SET_NODES")
 
    pop_sub(finufft_set_nodes)
  end subroutine finufft_set_nodes
 
  ! ---------------------------------------------------------
  subroutine finufft_modes_to_nodes(finufft, modes, values)
    type(finufft_t),             intent(in)  :: finufft
    complex(real64), contiguous, intent(in)  :: modes(:,:,:)
    complex(real64), contiguous, intent(out) :: values(:)
 
#ifdef HAVE_NFFT
    integer(c_int) :: ier
#endif
 
    push_sub(finufft_modes_to_nodes)
 
    call profiling_in("FINUFFT_MODES_TO_NODES")
 
#ifdef HAVE_NFFT
    ier = finufft_execute(finufft%plan, values, modes)
    call finufft_check(ier, 'finufft_execute')
#else
    values = m_z0
#endif
 
    call profiling_out("FINUFFT_MODES_TO_NODES")
 
    pop_sub(finufft_modes_to_nodes)
  end subroutine finufft_modes_to_nodes
 
  ! ---------------------------------------------------------
  subroutine finufft_nodes_to_modes(finufft, values, modes)
    type(finufft_t),             intent(in)  :: finufft
    complex(real64), contiguous, intent(in)  :: values(:)
    complex(real64), contiguous, intent(out) :: modes(:,:,:)
 
#ifdef HAVE_NFFT
    integer(c_int) :: ier
#endif
 
    push_sub(finufft_nodes_to_modes)
 
    call profiling_in("FINUFFT_NODES_TO_MODES")
 
#ifdef HAVE_NFFT
    ier = finufft_execute_adjoint(finufft%plan, values, modes)
    call finufft_check(ier, 'finufft_execute_adjoint')
#else
    modes = m_z0
#endif
 
    call profiling_out("FINUFFT_NODES_TO_MODES")
 
    pop_sub(finufft_nodes_to_modes)
  end subroutine finufft_nodes_to_modes
 
  ! ---------------------------------------------------------
  subroutine finufft_scale_nodes(X, v, inv_spacing, period)
    real(real64), contiguous, intent(in) :: X(:)
    real(real64), contiguous, intent(out) :: v(:)
    real(real64), intent(out) :: inv_spacing
    real(real64), intent(in)  :: period
 
    real(real64) :: length, cc
    real(real64) :: spacing(size(v) - 1)
    integer :: i
 
    push_sub(finufft_scale_nodes)
 
    assert(size(x) == size(v))
    assert(size(x) > 1)
 
    if (period > m_zero) then
      v = x/period
      inv_spacing = real(size(v), real64)
    else
      ! the sample nodes must be in [-1/2, 1/2)
      length = (maxval(x) - minval(x))*(m_one + m_epsilon)
      cc = (minval(x) + maxval(x))/m_two
      v = (x - cc)/length
 
      spacing = [(v(i+1) - v(i), i = 1, size(v)-1)]
      inv_spacing = m_one/minval(abs(spacing))
    end if
 
    pop_sub(finufft_scale_nodes)
  end subroutine finufft_scale_nodes
 
  ! ---------------------------------------------------------
  subroutine finufft_write_info(finufft)
    type(finufft_t), intent(in) :: finufft
 
    integer :: idir
 
    push_sub(finufft_write_info)
 
    call messages_write("Info: FINUFFT parameters")
    call messages_new_line()
 
    call messages_write("      Fourier coefficients      N = ")
    do idir = 1, finufft%dim
      call messages_write(finufft%N(idir))
      if (idir < finufft%dim) call messages_write(" x ")
    end do
    call messages_new_line()
 
    call messages_write("      Spatial nodes             M = ")
    do idir = 1, finufft%dim
      call messages_write(finufft%M(idir))
      if (idir < finufft%dim) call messages_write(" x ")
    end do
    call messages_new_line()
 
    call messages_write("      Requested tolerance    tol = ")
    call messages_write(finufft%tol)
    call messages_new_line()
 
    call messages_info()
 
    pop_sub(finufft_write_info)
  end subroutine finufft_write_info
 
  ! ---------------------------------------------------------
  subroutine finufft_end(finufft)
    type(finufft_t), intent(inout) :: finufft
 
#ifdef HAVE_NFFT
    integer(c_int) :: ier
#endif
 
    push_sub(finufft_end)
 
#ifdef HAVE_NFFT
    if (c_associated(finufft%plan)) then
      ier = finufft_destroy(finufft%plan)
      call finufft_check(ier, 'finufft_destroy')
      finufft%plan = c_null_ptr
    end if
#endif
 
    safe_deallocate_p(finufft%xj)
    safe_deallocate_p(finufft%yj)
    safe_deallocate_p(finufft%zj)
 
    pop_sub(finufft_end)
  end subroutine finufft_end
 
  ! ---------------------------------------------------------
#ifdef HAVE_NFFT
 
  subroutine finufft_check(ier, routine)
    integer(c_int), intent(in) :: ier
    character(*),   intent(in) :: routine
 
    if (ier == 0) return
 
    write(message(1), '(3a,i0,a)') "FINUFFT routine ", trim(routine), &
      " failed with error code ", ier, "."
    call messages_fatal(1)
  end subroutine finufft_check
#endif
 
#include "undef.F90"
#include "real.F90"
#include "finufft_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "finufft_inc.F90"
 
end module finufft_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: