Hi,

I've got s small problem with **ZFFT2DX** in **ACML 4.2.0**.

I wrote a sample program - all it does is it takes **2D** array of **real** data, then it takes forward Fourier Transform, multiplies it by **1/scale** and then takes backward Fourier Transform. In ideal case I would get original array. But, in real world I get tiny imaginary parts, which are close to 0. Is there any way to get rid of this tiny non-zero elements in the final result?

The reason I'm asking is that in my basic program I combine two real quantaties into complex number and perform the procedure described above (in order to take one FFT instead of two). And currently it means that one quantaty influences another.

I use ** ACML 4.2.0** in conjunction with

*(lates version) OR*

**ifort****. Same results for both compilers.**

*gfortran***Compiling lines:**

*gfortran -m64 -static acml_fft_check.f -L/opt/acml4.2.0/gfortran64/lib -lacml*

*ifort -O3 -static acml_fft_check.f -L/opt/acml4.2.0/ifort64/lib/ -lacml*

**The code of sample program:**

program compare_ffts

parameter (N=5)

implicit real*8 (a-h,o-z)

complex*16 COMM(N*N+6*N+200), m1(N,N)

complex*16 CMPT2(N,N), diff(N,N), CHECK(N,N)

real*8 cor1, sc, scale

logical ltrans, inpl

CMPT2=dcmplx(0.0d0,0.0d0)

sc=1.0d0

ltrans=.false.

inpl=.false.

call zfft2dx(100,sc,ltrans,inpl,N,N,CMPT1,1,N,m1,1,N,COMM,INFO)

do j=1,N

do i=1,N

CMPT2(i,j)=((dble(j)-1.0d0)*dble(N)+dble(i))*1.0d0

CHECK(i,j)=CMPT2(i,j)

end do

end do

sc=1.0d0

ltrans=.true.

inpl=.true.

CALL zfft2dx(-1,sc,ltrans,inpl,N,N,CMPT2,1,N,m1,1,N,COMM,INFO)

scale=N**2

do j=1,N

do i=1,N

CMPT2(i,j)=CMPT2(i,j)/scale

end do

end do

CALL zfft2dx(1,sc,ltrans,inpl,N,N,CMPT2,1,N,m1,1,N,COMM,INFO)

write(*,*) 'CMPT_ACML AFTER'

diff=CMPT2-CHECK

write(*,*) 'mean_real= ', sum(dble(diff))/scale

write(*,*) 'mean_imag= ', sum(dimag(diff))/scale

end

I should also mention that there is no problem if **N<=4.**

Also, these non-zero imaginary values are less if **N = 4, 16, 64, 256.** However further increase to **N=512, 1024** doesn't lead to decrease in the error.

In my basic program I use **N=512**.

**GCC/GFORTRAN:** version 4.3.2

**IFORT:** Version 11.0

The non-zero imaginary parts in the results that you observe are just due to rounding error. When the forward transform is computed, the results you get are not exact, they have been rounded to double precision - it's just that the rounding error is hard to see because the results are non-zero. And it's not just the final results that got rounded, but many intermediate results as well.

Then when you perform the backward FFT, the results get rounded again. This time the rounding errors stick out like a sore thumb because you know that you are expecting to get exact zeros in the solution. The magnitude of the non-zero parts - around 10^(-16) - is consistent with what you'd expect from rounding errors, because double precision floating-point arithmetic is accurate to around 16 decimal places.

You can easily observe the effects of rounding error in much simpler computations than an FFT - try running this program, which in exact arithmetic would produce a result of zero:

program t

double precision a, b, c

a = 7.0d0

b = 11.0d0

c = (1/a)*(1/b)*77.0d0 - 1.0d0

write(*,*) 'c = ', c

end

Compiling with gfortran or ifort, I get this result:

c = -1.11022302462515654E-016