Fortran - fftw - loosing accuracy in cosine transform

Author

fftw - loosing accuracy in cosine transform

Kamaraju Kusumanchi

2004-09-23, 3:56 am

Hi all,
Sorry if this is a wrong place for posting fftw related question.
But on the www.fftw.org, I could not find any email address other than
fftw@fftw.org which I think is meant to leave feedback rahter than
asking questions. In the past I got very good response on fft related
issues from this news group, so I am posting it here.

I am using fftw 3.0.1, absoft 8.0, fedora core 1. An example probably
explains my problem better. The output and sample code is given at the end.

My questions
(1) why I am getting only half the accuracy when doing a cosine
transform as opposed to the fouerier transform?
For example in fourier transform the round off error is O(1e-18) whereas
in cosine transform it is O(1e-9).
(2) Is it something that can be avoided or is it inherent in the cosine
transform algorithm?
Better still, can I improve the code to reduce the round-off errors in
the cosine transform?

output of the fourier transform
-4.92871531348332D-017 -6.22399809642460D-018
1.00000000000000 -1.03683365493055D-016
2.93123004154575D-017 1.25603286447767D-017
7.51060464024562D-018 7.65378971138986D-018
1.31628920003318D-017 1.45926836152971D-017
-8.53809210832335D-018 2.77146745197110D-017
1.54345126076430D-017 2.74725077800306D-018
1.81322035841033D-017 7.65378971138986D-018
6.22399809642460D-018 7.65378971138986D-018
1.84314369322536D-018 -1.11783997592600D-018
-1.26933184863306D-019 2.74725077800306D-018
-7.51060464024562D-018 7.65378971138986D-018
-7.14895807482630D-019 1.45926836152971D-017
-5.98208548668877D-017 -1.47589455874086D-017
-1.40047209926778D-017 1.25603286447767D-017
-1.81322035841033D-017 7.65378971138986D-018
output of the cosine/chebyshev transform
0.499999985996426
0.999999976990354
0.500000025327541
1.54879742012126D-008
-5.75090005171598D-009
3.17029677751287D-009
-2.08134160024177D-009
1.51237007364410D-009
-1.17638936092401D-009
9.62117874234755D-010
-8.18547139920208D-010
7.19508711332980D-010
-6.50485120223127D-010
6.02994678319307D-010
-5.71958809144919D-010
5.54385367040194D-010
-2.74346247976232D-010

program end_of_spectrum
use nrtype, WPC => DPC, WP => DP
implicit none
include "fftw3.f"

integer, parameter :: N=16
integer :: j, npo = N+1 ! npo = n plus one
integer(I8B) :: plan_fou, plan_cheb
complex(WPC) :: sample(0:N-1), output(0:N-1)
real(WP) :: cheb_sample(0:N), cheb_output(0:N)
! sample array contains input to the FFTW program
! The output of fftw is stored in output array
! Initialize sample to a simple function such that the actual fourier
! coeffficients will be zero after some inital modes.
! Then by observing the end values of output we know what the
round-off errors
! from the fftw program are
!
! cheb_sample and cheb_output contain the input and output to the
chebyshev
! transform.

call dfftw_plan_dft_1d( plan_fou, N, sample, output, FFTW_FORWARD,
FFTW_PATIENT)
call dfftw_plan_r2r_1d( plan_cheb, npo, cheb_sample, cheb_output, &
FFTW_REDFT00, FFTW_PATIENT )

! initialize sample
do j=0,N-1
sample(j) = cmplx(cos(2.0 * PI_D * j/N), sin(2.0 * PI_D * j /N))
end do

call dfftw_execute(plan_fou)
output = output/N

write(*,*) 'output of the fourier transform'
do j=0,N-1
write(*,*) real(output(j)), aimag(output(j))
end do

do j=0,N
! initialize the sample to x^2 + x
cheb_sample(j) = cos(j*PI/N) * (cos(j*PI/N) + 1.0)
end do

call dfftw_execute(plan_cheb)
cheb_output(1:N-1) = cheb_output(1:N-1)/N
cheb_output(0) = cheb_output(0)/(2.0 * N)
cheb_output(N) = cheb_output(N)/(2.0 * N)

write(*,*) 'output of the cosine/chebyshev transform'
do j=0,N
write(*,*) cheb_output(j)
end do

! conclusions:
! for 1D fourier dft using double precision, the ratio of the largest
to the
! smallest coefficients is of the order of 1e-16
! for 1D cosine transform the ratio of the largest to the smallest is
of the
! order of 1e-8
end program end_of_spectrum

James Van Buskirk

2004-09-23, 3:56 am

"Kamaraju Kusumanchi" <kk288@cornell.edu> wrote in message
news:cithr4$2p5$1@news01.cit.cornell.edu...

Afraid I can't help you much because I don't grok any of the
FFTW stuff, but I will make a couple of observations:

> ! initialize sample
> do j=0,N-1
> sample(j) = cmplx(cos(2.0 * PI_D * j/N), sin(2.0 * PI_D * j /N))
> end do

In the above I assume that sample is COMPLEX(8) and PI_D is REAL(8).
Did you know that if you don't supply a KIND optional argument to
the CMPLX intrinsic, the KIND of the result defaults to KIND(1.0)?
Surely not what you want here and you weren't punished with this
Fortran processor, but in the future supply a KIND = WPC keyword
actual argument at the end.

> do j=0,N
> ! initialize the sample to x^2 + x
> cheb_sample(j) = cos(j*PI/N) * (cos(j*PI/N) + 1.0)
> end do

Why do you persecute your presumed REAL(8) cheb_sample with
the presumed REAL(4) PI instead of PI_D? I betcha that if
you changed PI (resulting in invoking the single precision
flavor of the COS generic intrinsic) to PI_D (now invoking the
double precision flavor of COS) you would be the happier.
Also please insert the KIND optional argument to the CMPLX
intrinsic as noted above just to make me happy...

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Kamaraju Kusumanchi

2004-09-23, 3:58 pm

James Van Buskirk wrote:
> "Kamaraju Kusumanchi" <kk288@cornell.edu> wrote in message
> news:cithr4$2p5$1@news01.cit.cornell.edu...
>
> Afraid I can't help you much because I don't grok any of the
> FFTW stuff, but I will make a couple of observations:
>
>
>
>
> In the above I assume that sample is COMPLEX(8) and PI_D is REAL(8).
> Did you know that if you don't supply a KIND optional argument to
> the CMPLX intrinsic, the KIND of the result defaults to KIND(1.0)?
> Surely not what you want here and you weren't punished with this
> Fortran processor, but in the future supply a KIND = WPC keyword
> actual argument at the end.
>
>

Thanks for the subtle point. Though it is not affecting here, I am sure
it will be of significant impact in the future when I run the code on,
say a different compiler. Anyway, now I included kind argument in my
cmplx statements.

>
>
> Why do you persecute your presumed REAL(8) cheb_sample with
> the presumed REAL(4) PI instead of PI_D? I betcha that if
> you changed PI (resulting in invoking the single precision
> flavor of the COS generic intrinsic) to PI_D (now invoking the
> double precision flavor of COS) you would be the happier.
> Also please insert the KIND optional argument to the CMPLX
> intrinsic as noted above just to make me happy...
>

That did the trick. If I use PI_D instead of PI, I am getting the same
accuracy as that of the usual fourier transform. Now I feel stupid for
not being able to find this out in 2 days.

In case any one is interested, the results now look like
output of the fourier transform
-4.92871531348332D-017 -6.22399809642460D-018
1.00000000000000 -1.03683365493055D-016
2.93123004154575D-017 1.25603286447767D-017
7.51060464024562D-018 7.65378971138986D-018
1.31628920003318D-017 1.45926836152971D-017
-8.53809210832335D-018 2.77146745197110D-017
1.54345126076430D-017 2.74725077800306D-018
1.81322035841033D-017 7.65378971138986D-018
6.22399809642460D-018 7.65378971138986D-018
1.84314369322536D-018 -1.11783997592600D-018
-1.26933184863306D-019 2.74725077800306D-018
-7.51060464024562D-018 7.65378971138986D-018
-7.14895807482630D-019 1.45926836152971D-017
-5.98208548668877D-017 -1.47589455874086D-017
-1.40047209926778D-017 1.25603286447767D-017
-1.81322035841033D-017 7.65378971138986D-018
output of the cosine/chebyshev transform
0.500000000000000
1.00000000000000
0.500000000000000
-3.61606816949529D-017
-2.77555756156289D-017
-2.61405850032218D-017
1.47722546001150D-018
2.49773075486348D-017
2.77555756156289D-017
3.05338436826230D-017
-2.41641559192707D-017
-1.92061708609845D-017
-2.77555756156289D-017
-2.95148649033565D-017
-1.62223750058144D-017
-4.04949511423336D-017
1.38777878078145D-017

I am assuming there is no way to decrease these round off errors below
O(1e-17) using the double precision. Is that correct?

raju

Leif Harcke

2004-09-23, 3:58 pm

On Thu, 23 Sep 2004 12:00:34 -0400, Kamaraju Kusumanchi wrote:
> I am assuming there is no way to decrease these round off errors below
> O(1e-17) using the double precision. Is that correct?

< M A T L A B >
Copyright 1984-2004 The MathWorks, Inc.
Version 7.0.0.19901 (R14)
May 06, 2004

To get started, type one of these: helpwin, helpdesk, or demo.
For product information, visit www.mathworks.com.
[color=darkred]

ans =

1.1921e-07
[color=darkred]

ans =

2.2204e-16
[color=darkred]
EPS Spacing of floating point numbers.
D = EPS(X), is the positive distance from ABS(X) to the next larger in
magnitude floating point number of the same precision as X.
X may be either double precision or single precision.
For all X, EPS(X) = EPS(-X) = EPS(ABS(X)).

EPS, with no arguments, is the distance from 1.0 to the next larger double
precision number, that is EPS = 2^(-52).

EPS('double') is the same as EPS, or EPS(1.0).
EPS('single') is the same as EPS(single(1.0)), or single(2^-23).

[...]

James Van Buskirk

2004-09-23, 3:58 pm

"Kamaraju Kusumanchi" <kk288@cornell.edu> wrote in message
news:ciurv2$i0a$1@news01.cit.cornell.edu...

> That did the trick. If I use PI_D instead of PI, I am getting the same
> accuracy as that of the usual fourier transform. Now I feel stupid for
> not being able to find this out in 2 days.

When you have a symbol lying around called PI that isn't actually
the value of pi that you want (you wanted PI_D) it can create errors
that are absolutely invisible like this. Imagine what would have
happened if this weren't a test code where the consequences of
unitended use of the single precision version was so obvious. I
think somehow you should reconsider making the undesirable symbol PI
available -- maybe in your USE statement you could have an ONLY
clause that renames PI => PI_D.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Sponsored Links