Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/fft/tfft.f90
03540012fedefeature/refactoring

Last Change
Raw File

tfft.f90
View Options

!>
!> @file tfft.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> spclibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> spclibs 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 Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
! Test the FFT routines exported by module fft
!
USE fft
IMPLICIT NONE
INTEGER :: nx=8
DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE :: a,b,c
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: rn1, rn2
DOUBLE PRECISION :: pi, argx
INTEGER :: ix, nx0=2
!
PRINT*, 'Enter array dim. nx'
READ(*,*) nx
ALLOCATE(a(nx,4), b(nx,4), c(4,nx))
ALLOCATE(rn1(nx,4), rn2(nx,4))
!
! Create initial array
pi = 4.0d0*ATAN(1.0d0)
WRITE(*,*) 'Enter mode nx0'
READ(*,*) nx0
DO ix=0,nx-1
argx = 2.0d0*pi/nx*nx0*ix
a(ix+1,:) = COS(argx)
END DO
!
WRITE(*,*) '-----------------------'
WRITE(*,*) 'Reals of original array'
WRITE(*,*) '-----------------------'
WRITE(*,'(10f10.4)') REAL(a)
!________________________________________________________________________________
!
WRITE(*,*) '-----------------'
WRITE(*,*) 'testing four1 ...'
WRITE(*,*) '-----------------'
!
! Forward transform
b=a
CALL four1(b(:,1), -1)
!
! Check
WRITE(*,'(a/(10f10.4))') 'Reals of FFT', REAL(b(:,1))
WRITE(*,'(a/(10f10.4))') 'Imag of FFT', AIMAG(b(:,1))
!
! Backward transform
CALL four1(b(:,1), 1)
b = b/REAL(nx)
!
! Check
WRITE(*,'(a/(10f10.4))') 'Reals of FFT', REAL(b(:,1))
WRITE(*,'(a/(10f10.4))') 'Imag of FFT', AIMAG(b(:,1))
!
b(:,1) = b(:,1)-a(:,1)
rn1(:,1) = REAL(b(:,1),8)
rn2(:,1) = AIMAG(b(:,1))
PRINT *, 'Min. max err of real', MINVAL(rn1(:,1)), MAXVAL(rn1(:,1))
PRINT *, 'Min. max err of imag', MINVAL(rn2(:,1)), MAXVAL(rn2(:,1))
!________________________________________________________________________________
!
WRITE(*,*) '-------------------'
WRITE(*,*) 'testing fourcol ...'
WRITE(*,*) '-------------------'
!
! Forward transform
b=a
CALL fourcol(b(:,:), -1)
!
! Check
WRITE(*,'(a/(10f10.4))') 'Reals of FFT', REAL(b(:,1))
WRITE(*,'(a/(10f10.4))') 'Imag of FFT', AIMAG(b(:,1))
!
! Backward transform
CALL fourcol(b(:,:), 1)
b = b/REAL(nx)
!
! Check
WRITE(*,'(a/(10f10.4))') 'Reals of FFT', REAL(b(:,1))
WRITE(*,'(a/(10f10.4))') 'Imag of FFT', AIMAG(b(:,1))
!
b = b-a
rn1 = REAL(b,8)
rn2 = AIMAG(b)
PRINT *, 'Min. max err of real', MINVAL(rn1), MAXVAL(rn1)
PRINT *, 'Min. max err of imag', MINVAL(rn2), MAXVAL(rn2)
!________________________________________________________________________________
!
WRITE(*,*) '-------------------'
WRITE(*,*) 'testing fourrow ...'
WRITE(*,*) '-------------------'
!
! Forward transform
b=a
c = TRANSPOSE(b)
CALL fourrow(c, -1)
!
! Check
WRITE(*,'(a/(10f10.4))') 'Reals of FFT', REAL(c(1,:))
WRITE(*,'(a/(10f10.4))') 'Imag of FFT', AIMAG(c(1,:))
!
! Backward transform
CALL fourrow(c, 1)
c = c/REAL(nx)
!
! Check
WRITE(*,'(a/(10f10.4))') 'Reals of FFT', REAL(c(1,:))
WRITE(*,'(a/(10f10.4))') 'Imag of FFT', AIMAG(c(1,:))
!
b = TRANSPOSE(c)-a
rn1 = REAL(b,8)
rn2 = AIMAG(b)
PRINT *, 'Min. max err of real', MINVAL(rn1), MAXVAL(rn1)
PRINT *, 'Min. max err of imag', MINVAL(rn2), MAXVAL(rn2)
!
! Clean up
DEALLOCATE(a,b,c, rn1,rn2)
END PROGRAM main
c4science · Help