rlft3.f90 inverse transform

[dayakaran]

Hi all

I am trying to implement rlft3 subroutine to perform a inverse fourier transform. The size of my arrays are in powers of two, and my input complex array satisfies the required symmetry properties.

However the output I obtain is not strictly a real array. I tried to fourier transform the output, and it gives the input array.

So I thought may there is a mistake with how I set up my speq array.

do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
speq(j,k)=data0(0,jc,kc)
enddo
enddo

do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
enddo
enddo

function indc(i,N)
!
implicit none
integer:: indc,i,N
!
if(i<=N/2) then
indc=i-1
else
indc=i-1-N
endif
!
return
end function indc

I ensure symmetry in the code too. Is there is a mistake with what I am doing ? kindly help me in this regard

Thanks

John

[davekw7x]

Originally Posted by dayakaran

...
I am trying to implement rlft3 subroutine to perform a inverse fourier transform.

So: You have "something" that is supposed to be the Fourier Transform of a 3-D real array. Right? (In the f90 distribution, if it's for a 2-D real array you would use rlft2. Right?)

Anyhow, where did the "something" come from? That is, where and how did you get the data that comprises the 3-D complex collection of data values that supposedly represents the FFT of a 3-D real function? (See Footnote.)

Maybe you can show us how you fill the data array that you are using with rlft3 from the complex data values.

...

Originally Posted by dayakaran

However the output I obtain is not strictly a real array.

I don't understand what you mean here. Maybe you can show array declarations and show how you call rlft3 and how you interpret the output as "not being strictly a real array."

Originally Posted by dayakaran

...
speq(j,k)=data0(0,jc,kc)
...

Regardless of the correctness of other things in your code, if data0 is a 3-D Fortran array that you are going to use with NR f90 routines, you can't have any index value of zero. Or is data0 some kind of function. Or what?

Originally Posted by dayakaran

...Is there is a mistake with what I am doing ?...

Maybe you can show some more of your work.

Regards,

Dave

Footnote: In order to understand what the heck is going on with the data array and spec and speq, some people find it instructive to start with a real data array and perform the forward FFT. (Use a non-trivial something for which is easy to predict the results.) Look at the data array and look at speq after the call to rlft2 or rlft3 and make sure you can see exactly what happened.

[dayakaran]

Hi Dave

Thanks for your detailed reply.

I ve attached my code for you to look at

here is my code

program fftcheck
!
implicit none
!
integer, parameter :: Nv=8
complex,dimension(Nv,Nv) :: speq
complex,dimension(-Nv/2:Nv/2,-Nv/2:Nv/2,-Nv/2:Nv/2) :: data0
complex,dimension(Nv/2,Nv,Nv) :: data
double precision :: ran1,f1,f2
integer :: idum,i,j,k,indc,jc,kc
!
call system_clock(idum)
!
open(unit=220,file='input_speq.dat',status='unknow n')
open(unit=221,file='outpu_real.dat',status='unknow n')
open(unit=222,file='real_inver.dat',status='unknow n')
open(unit=223,file='freq_nyqst.dat',status='unknow n')
!
do i=-Nv/2,Nv/2
do j=-Nv/2,Nv/2
do k=-Nv/2,Nv/2
f1=ran1(idum)
f2=ran1(idum)
data0(i,j,k)=cmplx(f1,f2)
enddo
enddo
enddo
!
! Symmetry due to the realness of data in real-space
!
do i=-Nv/2,Nv/2
do j=-Nv/2,Nv/2
do k=-Nv/2,Nv/2
data0(i,j,k)=conjg(data0(-i,-j,-k))
enddo
enddo
enddo
!
! Store the workspace arrays for fft
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
data(i,j,k)=data0(i-1,jc,kc)
enddo
enddo
enddo
!
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
speq(j,k)=data0(-Nv/2,jc,kc)
enddo
enddo
!
! We try to retain symmetry here
!
data(1,1,1)=cmplx(real(data(1,1,1)),0.0) ! global average field
data(1,Nv/2+1,1)=cmplx(real(data(1,Nv/2+1,1)),0.0)
data(1,1,Nv/2+1)=cmplx(real(data(1,1,Nv/2+1)),0.0)
data(1,Nv/2+1,Nv/2+1)=cmplx(real(data(1,Nv/2+1,Nv/2+1)),0.0)
speq(1,1)=cmplx(real(speq(1,1)),0.0)
speq(Nv/2+1,1)=cmplx(real(speq(Nv/2+1,1)),0.0)
speq(1,Nv/2+1)=cmplx(real(speq(1,Nv/2+1)),0.0)
speq(Nv/2+1,Nv/2+1)=cmplx(real(speq(Nv/2+1,Nv/2+1)),0.0)
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(220,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i-1,jc,kc
enddo
enddo
enddo
!
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc
enddo
enddo
!
! Inverse FFT ing the data
!
call rlft3(data,speq,Nv,Nv,Nv,-1)
!
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc
enddo
enddo
!
! FFT
!
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(221,'(2f16.8,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i-1,jc,kc
enddo
enddo
enddo
!
call rlft3(data,speq,Nv,Nv,Nv,1)
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(222,'(2f16.8,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i-1,jc,kc
enddo
enddo
enddo
!
stop
!
end program fftcheck
!
!************************************
!************************************
!
function indc(i,N)
!
implicit none
integer:: indc,i,N
!
if(i<=N/2) then
indc=i-1
else
indc=i-1-N
endif
!
return
end function indc
!
!************************************
!************************************



Thanks a lot for your help

Cheers

John

[dayakaran]

And to check if my rlft3.f90 works right, I took real data, fourier transformed it, Then inverse fourier transformed it again to check if things were fine, and it worked right

I ve attached the code here


program fftcheck
!
implicit none
!
integer, parameter :: Nv=16
complex,dimension(Nv,Nv) :: speq
complex,dimension(-Nv/2:Nv/2,-Nv/2:Nv/2,-Nv/2:Nv/2) :: data0
complex,dimension(Nv/2,Nv,Nv) :: data
double precision :: ran1,f1,f2
integer :: idum,i,j,k,indc,jc,kc,finder
!
call system_clock(idum)
!
write(*,*) 'Randomn number is ',idum
!
open(unit=220,file='input_real.dat',status='unknow n')
open(unit=221,file='outpu_speq.dat',status='unknow n')
open(unit=222,file='speq_inver.dat',status='unknow n')
open(unit=223,file='freq_nyqst.dat',status='unknow n')
!
do i=-Nv/2,Nv/2
do j=-Nv/2,Nv/2
do k=-Nv/2,Nv/2
f1=ran1(idum)
f2=0.0d0
data0(i,j,k)=cmplx(f1,f2)
enddo
enddo
enddo
!
!
! Store the workspace arrays for fft
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
data(i,j,k)=data0(i-1,jc,kc)
enddo
enddo
enddo
!
!
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(220,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i-1,jc,kc
enddo
enddo
enddo
!
call rlft3(data,speq,Nv,Nv,Nv,1)
!
do jc=1,Nv
j=finder(jc,Nv)
do kc=1,Nv
k=finder(kc,Nv)
if (data(1,jc,kc).ne.conjg(data(1,j,k))) write(*,*) 'Hello',jc,j,kc,k
enddo
enddo
!
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc
enddo
enddo
!
! Finally do the FFT
!
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(221,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i-1,jc,kc
enddo
enddo
enddo
!
call rlft3(data,speq,Nv,Nv,Nv,-1)
!
do i=1,Nv/2
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(222,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i-1,jc,kc
enddo
enddo
enddo
!
do j=1,Nv
jc=indc(j,Nv)
do k=1,Nv
kc=indc(k,Nv)
write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc
enddo
enddo
!
stop
!
end program fftcheck
!
!************************************
!************************************
function indc(i,N)
!
implicit none
integer:: indc,i,N
!
if(i<=N/2) then
indc=i-1
else
indc=i-1-N
endif
!
return
end function indc
!
!************************************
!************************************
!

Thank you once again

John

[davekw7x]

Originally Posted by dayakaran

...
I ve attached my code for you to look at...

Here's what I was trying to suggest for starters: Do the forward transform by calling filling the 8x8x8 real data matrix and calling rlft3(,,,+1). Inspect the results (the transformed 4x8x8 complex data matrix and the 8x8 complex speq matrix).

By my reckoning, here is a list of members of the transformed data matrix and the speq matrix that satisfy the symmetry requirements. (That is, I show index values of pairs that are complex conjugates.)

If the members of the matrices do not have the conjugate relationship indicated in the table below, the matrices might not represent the results of a 3-d Fourier transformation of a matrix of real numbers (so trying rlft3(,,,-1) won't give anything really meaningful in the context of what I perceive you are asking.)

Number of complex data matrix values = 256
Here is a list of data matrix conjugates:
* (1,1,1)  (1,1,1) *
  (1,1,2)  (1,1,8)
  (1,1,3)  (1,1,7)
  (1,1,4)  (1,1,6)
* (1,1,5)  (1,1,5) *
  (1,2,1)  (1,8,1)
  (1,2,2)  (1,8,8)
  (1,2,3)  (1,8,7)
  (1,2,4)  (1,8,6)
  (1,2,5)  (1,8,5)
  (1,2,6)  (1,8,4)
  (1,2,7)  (1,8,3)
  (1,2,8)  (1,8,2)
  (1,3,1)  (1,7,1)
  (1,3,2)  (1,7,8)
  (1,3,3)  (1,7,7)
  (1,3,4)  (1,7,6)
  (1,3,5)  (1,7,5)
  (1,3,6)  (1,7,4)
  (1,3,7)  (1,7,3)
  (1,3,8)  (1,7,2)
  (1,4,1)  (1,6,1)
  (1,4,2)  (1,6,8)
  (1,4,3)  (1,6,7)
  (1,4,4)  (1,6,6)
  (1,4,5)  (1,6,5)
  (1,4,6)  (1,6,4)
  (1,4,7)  (1,6,3)
  (1,4,8)  (1,6,2)
* (1,5,1)  (1,5,1) *
  (1,5,2)  (1,5,8)
  (1,5,3)  (1,5,7)
  (1,5,4)  (1,5,6)
* (1,5,5)  (1,5,5) *

Number of complex speq matrix values = 64
Here is a list of speq matrix conjugates:
* (1,1)  (1,1) *
  (1,2)  (1,8)
  (1,3)  (1,7)
  (1,4)  (1,6)
* (1,5)  (1,5) *
  (2,1)  (8,1)
  (2,2)  (8,8)
  (2,3)  (8,7)
  (2,4)  (8,6)
  (2,5)  (8,5)
  (2,6)  (8,4)
  (2,7)  (8,3)
  (2,8)  (8,2)
  (3,1)  (7,1)
  (3,2)  (7,8)
  (3,3)  (7,7)
  (3,4)  (7,6)
  (3,5)  (7,5)
  (3,6)  (7,4)
  (3,7)  (7,3)
  (3,8)  (7,2)
  (4,1)  (6,1)
  (4,2)  (6,8)
  (4,3)  (6,7)
  (4,4)  (6,6)
  (4,5)  (6,5)
  (4,6)  (6,4)
  (4,7)  (6,3)
  (4,8)  (6,2)
* (5,1)  (5,1) *
  (5,2)  (5,8)
  (5,3)  (5,7)
  (5,4)  (5,6)
* (5,5)  (5,5) *

Note that, for emphasis, I flagged the values that are real numbers. (A complex number is its own conjugate if and only if its imaginary part is equal to zero.)

Of course a more formal statement and derivation of the conditions would be appropriate, but I think you might get started by looking at some actual results. Maybe you will test a few and then try some actual derivations. Try forcing those conjugate relationships onto your randomly generated matrices and see if rlft3(,,,-1) followed by rlft3(,,,+1) show you what you want to see.

Maybe there is more to the story (and maybe there is less).

Bottom line: don't take my word for it. In other words: Your Mileage May Vary.

Regards,

Dave