nfft.F90

Go to the documentation of this file.

!! Copyright (C) 2011 U. De Giovannini
!!
!! 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 nfft_oct_m
  use fftw_params_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_c_binding
  use loct_math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use varinfo_oct_m
  implicit none
 
  private
 
  public ::          &
    nfft_t,          &
    nfft_copy_info,  &
    nfft_init,       &
    nfft_end,        &
    nfft_precompute, &
    nfft_write_info, &
    nfft_guru_options, &
    znfft_forward,   &
    znfft_backward,  &
    dnfft_forward,   &
    dnfft_backward
 
  ! global constants
  integer, public, parameter ::         &
    nfft_real    = 0,                   &
    nfft_complex = 1
 
  !NFFT flags
  integer, public, parameter ::        &
    NFFT_PRE_PHI_HUT       =        0, &
    nfft_fg_psi            =        2, &
    nfft_pre_lin_psi       =        4, &
    nfft_pre_fg_psi        =        8, &
    nfft_pre_psi           =       16, &
    nfft_pre_full_psi      =       32, &
    nfft_malloc_x          =       64, &
    nfft_malloc_f_hat      =      128, &
    nfft_malloc_f          =      256, &
    nfft_fft_out_of_place  =      512, &
    nfft_fftw_init         =     1024
 
 
  type nfft_t
    private
 
    integer           :: N(3)
    integer           :: M(3)
    integer           :: dim
    integer           :: fftN(3)
    real(real64), public     :: norm
 
    ! Guru options
    logical, public   :: set_defaults = .false. 
    logical, public   :: guru
    integer, public   :: precompute
    integer, public   :: mm
    real(real64),   public   :: sigma
 
    type(c_ptr)       :: plan
 
  end type nfft_t
 
 
contains
 
 
 
  ! ---------------------------------------------------------
  ! GURU options
  subroutine nfft_guru_options(nfft, namespace)
    type(nfft_t),      intent(inout) :: nfft
    type(namespace_t), intent(in)    :: namespace
 
    push_sub(nfft_guru_options)
 
    !%Variable NFFTGuruInterface
    !%Type logical
    !%Default false
    !%Section Mesh::FFTs
    !%Description
    !% Perform NFFT with guru interface. This permits the fine tuning of several critical parameters.
    !%End
    call parse_variable(namespace, 'NFFTGuruInterface',  .false., nfft%guru)
 
 
    !%Variable NFFTCutoff
    !%Type integer
    !%Default 6
    !%Section Mesh::FFTs
    !%Description
    !% Cut-off parameter of the window function.
    !% See NFFT manual for details.
    !%End
    call parse_variable(namespace, 'NFFTCutoff', 6, nfft%mm)
 
 
    !%Variable NFFTOversampling
    !%Type float
    !%Default 2
    !%Section Mesh::FFTs
    !%Description
    !% NFFT oversampling factor (sigma). This will rule the size of the FFT under the hood.
    !%End
    call parse_variable(namespace, 'NFFTOversampling', m_two, nfft%sigma)
 
    !%Variable NFFTPrecompute
    !%Type integer
    !%Default NFFT_PRE_PSI
    !%Section Mesh::FFTs
    !%Description
    !% NFFT precomputation strategy.
    !%Option NFFT_PRE_LIN_PSI 4
    !% This method implements a linear interpolation from a lookup table.
    !%Option NFFT_PRE_PSI 16
    !% This method uses a medium amount of memory to store d*(2*m+1)*M real numbers and requires at most
    !% 2(2m + 1)d extra multiplications for each node.
    !% This is the default option.
    !%Option NFFT_PRE_FULL_PSI 32
    !% Is the fastest method but requires a large amount of memory as it requires to store (2*m+1)^d*M
    !% real numbers. No extra operations are needed during matrix vector multiplication.
    !%End
    call parse_variable(namespace, 'NFFTPrecompute', nfft_pre_psi, nfft%precompute)
    if (.not. varinfo_valid_option('NFFTPrecompute', nfft%precompute)) call messages_input_error(namespace, 'NFFTPrecompute')
!    call messages_print_var_option("NFFTPrecompute", nfft%precompute, namespace=namespace)
 
!     if (.not. varinfo_valid_option('NFFTPrecompute', nfft%precompute, is_flag=.true.)) then
!       call messages_input_error('NFFTPrecompute')
!     end if
 
 
    pop_sub(nfft_guru_options)
  end subroutine nfft_guru_options
 
 
 
  ! ---------------------------------------------------------
  subroutine nfft_init(nfft, nfft_options, N, dim, M, optimize)
    type(nfft_t),      intent(inout) :: nfft
    type(nfft_t),      intent(in)    :: nfft_options
    integer,           intent(inout) :: n(3)
    integer,           intent(inout) :: m(3)
    integer,           intent(in)    :: dim
    logical, optional, intent(in)    :: optimize
 
    integer :: ii, my_n(3)
    logical :: optimize_
    integer :: nfft_flags
 
 
    push_sub(nfft_init)
 
    optimize_ = optional_default(optimize, .true.)
 
    nfft%dim = dim
    nfft%M(:) = m(:)
    nfft%N(:) = n(:)
 
    if (.not. nfft%set_defaults) then
      nfft%guru = nfft_options%guru
      nfft%mm = nfft_options%mm
      nfft%sigma = nfft_options%sigma
      nfft%precompute = nfft_options%precompute
    end if
 
    ! set unused dimensions to 1
    nfft%M(dim+1:3) = 1
 
    my_n = 0
    do ii = 1, dim
      my_n(ii) = n(ii)*int(nfft%sigma)
      if (optimize_ .or. (.not. nfft%guru)) call loct_fft_optimize(my_n(ii), 1) ! ask for an odd number
    end do
 
    nfft%fftN(1:dim) = my_n(1:dim)
 
    if (nfft%guru) then ! Guru interface
      nfft_flags =  nfft_pre_phi_hut  + nfft_malloc_x +nfft_malloc_f_hat +&
        nfft_malloc_f + nfft_fftw_init + nfft_fft_out_of_place
 
      nfft_flags = nfft_flags + nfft%precompute
 
      call oct_nfft_init_guru(nfft%plan, dim, n, m(1)*m(2)*m(3), my_n, nfft%mm, &
        nfft_flags, fftw_measure + fftw_destroy_input)
 
    else ! Default interfaces
 
      select case (dim)
      case (3)
        call oct_nfft_init_3d(nfft%plan, n(1), n(2),n(3), m(1)*m(2)*m(3))
      case (2)
        call oct_nfft_init_2d(nfft%plan, n(1), n(2), m(1)*m(2))
      case (1)
        call oct_nfft_init_1d(nfft%plan,n(1),m(1))
      end select
 
    end if
 
 
    pop_sub(nfft_init)
  end subroutine nfft_init
 
  ! ---------------------------------------------------------
  subroutine nfft_write_info(nfft)
    type(nfft_t), intent(inout) :: nfft
 
    integer :: idir
!    integer :: mm
 
    push_sub(nfft_write_info)
 
    call messages_write("Info: NFFT parameters")
    call messages_new_line()
    call messages_write("      Fourier coefficients      N = ")
    do idir = 1,  nfft%dim
      call messages_write(nfft%N(idir))
      if (idir < nfft%dim) call messages_write(" x ")
    end do
    call messages_new_line()
 
    call messages_write("      Spatial nodes             M = ")
 
!     mm = nfft%M(1)*nfft%M(2)*nfft%M(3)
!
!     call messages_write(mm)
!     call messages_new_line()
    do idir = 1,  nfft%dim
      call messages_write(nfft%M(idir))
      if (idir < nfft%dim) call messages_write(" x ")
    end do
    call messages_new_line()
 
 
    call messages_write("      Oversampling factor   sigma = ")
    call messages_write(nfft%sigma)
    call messages_new_line()
 
    call messages_write("      FFT grid size             n = ")
    do idir = 1,  nfft%dim
      call messages_write(nfft%fftN(idir))
      if (idir < nfft%dim) call messages_write(" x ")
    end do
    call messages_new_line()
 
    call messages_write("      Real space cutoff         m = ")
    call messages_write(nfft%mm)
    call messages_new_line()
 
    call messages_write("      Pre-computation strategy    = ")
    select case (nfft%precompute)
    case (nfft_pre_lin_psi)
      call messages_write(" NFFT_PRE_LIN_PSI")
    case (nfft_pre_psi)
      call messages_write(" NFFT_PRE_PSI")
    case (nfft_pre_full_psi)
      call messages_write(" NFFT_PRE_FULL_PSI")
    end select
 
    call messages_info()
 
 
    pop_sub(nfft_write_info)
  end subroutine nfft_write_info
 
 
  ! ---------------------------------------------------------
  subroutine nfft_end(nfft)
    type(nfft_t), intent(inout) :: nfft
 
    push_sub(nfft_end)
 
    call oct_nfft_finalize(nfft%plan)
 
    pop_sub(nfft_end)
  end subroutine nfft_end
 
  ! ---------------------------------------------------------
  ! This routine is intend to copy the configuration parameters
  ! rather the whole structure.
  ! ---------------------------------------------------------
  subroutine nfft_copy_info(in, out)
    type(nfft_t), intent(in)  :: in
    type(nfft_t), intent(out) :: out
 
 
    push_sub(nfft_copy_info)
 
    out%N = in%N
    out%M = in%M
    out%dim = in%dim
    out%fftN = in%fftN
    out%norm = in%norm
 
    out%guru = in%guru
    out%precompute = in%precompute
    out%mm = in%mm
    out%sigma = in%sigma
 
 
    pop_sub(nfft_copy_info)
  end subroutine nfft_copy_info
 
 
  !----------------------------------------------------------
  ! Precompute the plan according to the position the grid nodes in real space
  ! x axis is X1, y axis is X2, z axis is X3
  ! NOTE: We only allow different spacing for each direction x,y,z
  ! the NFFT interface however is more general
  ! ---------------------------------------------------------
  subroutine nfft_precompute(nfft, X1, X2, X3)
    type(nfft_t),    intent(inout) :: nfft
    real(real64),    intent(in)    :: x1(:)
    real(real64), optional, intent(in)    :: x2(:)
    real(real64), optional, intent(in)    :: x3(:)
 
 
 
    real(real64)   :: x1_(1:nfft%m(1)), x2_(1:nfft%m(2)), x3_(1:nfft%m(3))
    real(real64)   :: length, cc, eps, dx1(1:nfft%m(1)-1),  dx2(1:nfft%m(2)-1), dx3(1:nfft%m(3)-1)
 
    integer :: ii
 
    push_sub(nfft_precompute)
 
!     eps = 1.000001 ! the sample nodes must be in [0.5,0.5)
    eps = m_one+m_epsilon ! the sample nodes must be in [0.5,0.5)
 
    select case (nfft%dim)
    case (3)
      length = (maxval(x1)-minval(x1))*eps
      cc = (minval(x1)+maxval(x1))/m_two
      x1_ =(x1-cc)/length
      length = (maxval(x2)-minval(x2))*eps
      cc = (minval(x2)+maxval(x2))/m_two
      x2_ =(x2-cc)/length
      length = (maxval(x3)-minval(x3))*eps
      cc = (minval(x3)+maxval(x3))/m_two
      x3_ =(x3-cc)/length
      call oct_nfft_precompute_one_psi_3d(nfft%plan, nfft%M, x1_, x2_, x3_)
 
      ! Set the normalization factor
      do ii = 1, nfft%M(1)-1
        dx1(ii)= abs(x1_(ii+1)-x1_(ii))
      end  do
      do ii = 1, nfft%M(2)-1
        dx2(ii)= abs(x2_(ii+1)-x2_(ii))
      end do
      do ii = 1, nfft%M(3)-1
        dx3(ii)= abs(x3_(ii+1)-x3_(ii))
      end do
      nfft%norm = m_one/(minval(dx1(:)) * minval(dx2(:)) * minval(dx3(:)))
 
    case (2)
      length = (maxval(x1)-minval(x1))*eps
      cc = (minval(x1)+maxval(x1))/m_two
      x1_ =(x1-cc)/length
      length = (maxval(x2)-minval(x2))*eps
      cc = (minval(x2)+maxval(x2))/m_two
      x2_ =(x2-cc)/length
      call oct_nfft_precompute_one_psi_2d(nfft%plan, nfft%M, x1_, x2_)
 
      ! Set the normalization factor
      do ii = 1, nfft%M(1)-1
        dx1(ii)= abs(x1_(ii+1)-x1_(ii))
      end do
      do ii = 1, nfft%M(2)-1
        dx2(ii)= abs(x2_(ii+1)-x2_(ii))
      end do
      nfft%norm = m_one/(minval(dx1(:)) * minval(dx2(:)))
 
 
    case (1)
      length = (maxval(x1)-minval(x1))*eps
      cc = (minval(x1)+maxval(x1))/m_two
      x1_ =(x1-cc)/length
      call oct_nfft_precompute_one_psi_1d(nfft%plan,nfft%M(1),x1_)
 
      ! Set the normalization factor
      do ii = 1, nfft%M(1)-1
        dx1(ii)= abs(x1_(ii+1)-x1_(ii))
      end do
      nfft%norm = m_one/(minval(dx1(:)))
 
    end select
 
    ! check the plan
    call oct_nfft_check(nfft%plan)
 
 
    write(message(1), '(a)') "Info: NFFT plan precomputed."
    call messages_info(1)
 
 
    pop_sub(nfft_precompute)
  end subroutine nfft_precompute
 
#include "undef.F90"
#include "real.F90"
#include "nfft_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "nfft_inc.F90"
 
end module nfft_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: