!===============================================================================
! Copyright 2011-2017 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

! Content:
!       Example of using dfftw_plan_dft_2d function.
!
!*****************************************************************************

program dp_plan_dft_2d

  include 'fftw3.f'

  ! Sizes of 2D transform
  integer, parameter :: N1 = 3
  integer, parameter :: N2 = 3

  ! Arbitrary harmonic used to verify FFT
  !integer, parameter :: H1 = 1
  !integer, parameter :: H2 = -N2/2

  ! Need double precision
  integer, parameter :: WP = selected_real_kind(15,307)

  ! Execution status
  integer :: status = 0

  ! Data array
  complex(WP), allocatable :: x(:,:)

  ! FFTW plan
  integer*8 :: fwd = 0, bwd = 0

  print *,"Example dp_plan_dft_2d"
  print *,"Forward and backward 2D complex in-place FFT"
  print *,"Configuration parameters:"
  print '("  N = ["I0","I0"]")', N1, N2
  !print '("  H = ["I0","I0"]")', H1, H2

  print *,"Allocate array"
  allocate ( x(N1,N2), STAT = status)
  if (0 /= status) goto 999

  print *,"Create FFTW forward transform plan"
  call dfftw_plan_dft_2d(fwd, N1, N2, x, x, FFTW_FORWARD, FFTW_ESTIMATE)
  if (0 == fwd) goto 999

  print *,"Create FFTW backward transform plan"
  call dfftw_plan_dft_2d(bwd, N1, N2, x, x, FFTW_BACKWARD, FFTW_ESTIMATE)
  if (0 == bwd) goto 999

  print *,"Initialize data for forward transform"
  !call init(x, N1, N2, H1, H2)
  
  WRITE(*,*) shape(x)

  x = 1.0_WP

  WRITE(*,'(6F4.1)') x
    
  !x(1,1)= 0.1D0; x(1,2)= 0.2D0;  x(1,3)= 0.3D0; 
  !x(2,1)= 0.4D0; x(2,2)= 0.5D0;  x(2,3)= 0.6D0; 
  !x(3,1)= 0.7D0; x(3,2)= 0.8D0;  x(3,3)= 0.9D0; 
  
  !write(*,'(3F8.1)') x

  print *,"Compute forward transform"
  call dfftw_execute(fwd)
  
  write(*,'(6F4.1)') x

  !print *,"Verify the result of the forward transform"
  !status = verificate(x, N1, N2, H1, H2)
  !if (0 /= status) goto 999

  !print *,"Initialize data for backward transform"
  !call init(x, N1, N2, -H1, -H2)

  print *,"Compute backward transform"
  call dfftw_execute_dft(bwd, x, x)

  write(*,'(6F4.1)') x

  !print *,"Verify the result of the backward transform"
  !status = verificate(x, N1, N2, H1, H2)
  !if (0 /= status) goto 999

100 continue

  print *,"Destroy FFTW plans"
  call dfftw_destroy_plan(fwd)
  call dfftw_destroy_plan(bwd)

  print *,"Deallocate array"
  deallocate(x)

  if (status == 0) then
    print *, "TEST PASSED"
    call exit(0)
  else
    print *, "TEST FAILED"
    call exit(1)
  end if

999 print '("  Error, status = ",I0)', status
  status = 1
  goto 100

contains

  ! Compute mod(K*L,M) accurately
  pure integer*8 function moda(k,l,m)
    integer, intent(in) :: k,l,m
    integer*8 :: k8
    k8 = k
    moda = mod(k8*l,m)
  end function moda

  ! Initialize array with harmonic /H1, H2/
  subroutine init(x, N1, N2, H1, H2)
    integer N1, N2, H1, H2
    complex(WP) :: x(:,:)

    integer k1, k2
    complex(WP), parameter :: I_TWOPI = (0,6.2831853071795864769_WP)

    forall (k1=1:N1, k2=1:N2)
      x(k1,k2) = exp( I_TWOPI * ( &
        real(  moda(k1-1,H1, N1), WP)/N1 &
        + real(moda(k2-1,H2, N2), WP)/N2 )) / (N1*N2)
    end forall
  end subroutine init

  ! Verify that x(N1,N2) is unit peak at x(H1,H2)
  integer function verificate(x, N1, N2, H1, H2)
    integer N1, N2, H1, H2
    complex(WP) :: x(:,:)

    integer k1, k2
    real(WP) err, errthr, maxerr
    complex(WP) :: res_exp, res_got

    ! Note, this simple error bound doesn't take into account error of
    ! input data
    errthr = 5.0 * log(real(N1*N2,WP)) / log(2.0_WP) * EPSILON(1.0_WP)
    print '("  Check if err is below errthr " G10.3)', errthr

    maxerr = 0
    do k2 = 1, N2
      do k1 = 1, N1
        if (mod(k1-1-H1, N1)==0 .AND. mod(k2-1-H2, N2)==0) then
          res_exp = 1
        else
          res_exp = 0
        end if
        res_got = x(k1,k2)
        err = abs(res_got - res_exp)
        maxerr = max(err,maxerr)
        if (.not.(err < errthr)) then
          print '("  x("I0","I0"):"$)', k1, k2
          print '(" expected ("G24.17","G24.17"),"$)', res_exp
          print '(" got ("G24.17","G24.17"),"$)', x(k1,k2)
          print '(" err "G10.3)', err
          print *," Verification FAILED"
          verificate = 1
          return
        end if
      end do
    end do
    print '("  Verified,  maximum error was " G10.3)', maxerr
    verificate = 0
  end function verificate

end program dp_plan_dft_2d
