bug report for CMBFAST and CAMB
-
- Posts: 12
- Joined: April 22 2006
- Affiliation: Sun Yat-Sen University
- Contact:
bug report for CMBFAST and CAMB
The subroutine bjl, which calculates the spherical bessel functions in package CMBFAST and CAMB, has a small bug.
Use the following codes to see this bug.
!!====for CMBFAST PACKAGE (use the subroutine bjl in file jlgen.F) ==!
program main
real x,jl
integer L
L=5
x=1.E-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
x=1.001E-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
end
!!===for CAMB PACKAGE(use the subroutine bjl in file bessels.f90) ==!!
program main
doubleprecision x,jl
integer L
L=5
x=1.D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
x=1.001D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
end program
!====end of the test codes=====!
Go to my web http://www.astro.utoronto.ca/~zqhuang/a ... ademic.htm to download the corrected subroutines for spherical bessel functions.
Use the following codes to see this bug.
!!====for CMBFAST PACKAGE (use the subroutine bjl in file jlgen.F) ==!
program main
real x,jl
integer L
L=5
x=1.E-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
x=1.001E-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
end
!!===for CAMB PACKAGE(use the subroutine bjl in file bessels.f90) ==!!
program main
doubleprecision x,jl
integer L
L=5
x=1.D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
x=1.001D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
end program
!====end of the test codes=====!
Go to my web http://www.astro.utoronto.ca/~zqhuang/a ... ademic.htm to download the corrected subroutines for spherical bessel functions.
-
- Posts: 1943
- Joined: September 23 2004
- Affiliation: University of Sussex
- Contact:
Re: bug report for CMBFAST and CAMB
The output is
x= 1.000000000000000E-003 j_5(x)= 9.620009620009621E-020
x= 1.001000000000000E-003 j_5(x)= 0.000000000000000E+000
so I doubt it has any significant effect since the Bessel functions are tiny anyway; for a particular model I tried the change to the C_l was O(10^{-4}).
x= 1.000000000000000E-003 j_5(x)= 9.620009620009621E-020
x= 1.001000000000000E-003 j_5(x)= 0.000000000000000E+000
so I doubt it has any significant effect since the Bessel functions are tiny anyway; for a particular model I tried the change to the C_l was O(10^{-4}).
-
- Posts: 12
- Joined: April 22 2006
- Affiliation: Sun Yat-Sen University
- Contact:
bug report for CMBFAST and CAMB
The output of j_5(0.001001) depends on the precision of your computer. On my own laptop it's -124. On the computer in CITA it is -87. The fomula used in that subroutine for l=5, x>0.001 is an exact expression. But this does not mean that the result would always be correct. The problem is caused by substracting one huge number from another. If you are using a super computer, then it is fine. My suggestion is to try this on your own machine or somewhere else, you can see the difference.
-
- Posts: 12
- Joined: April 22 2006
- Affiliation: Sun Yat-Sen University
- Contact:
bug report for CMBFAST and CAMB
Oh, I read the CAMB code and find that it is different from the codes in CMBFAST and for x=0.001001 the output is 0.E0. Then it is fine. But if you try CMBFAST, that is really a problem.
-
- Posts: 12
- Joined: April 22 2006
- Affiliation: Sun Yat-Sen University
- Contact:
bug report for CMBFAST and CAMB
for CAMB, the code is less problematic, but still not perfect, let's try this:
program main
doubleprecision x,jl
integer L
L=5
x=1.D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
x=1.00000001D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
end
Then you will see the problem.
program main
doubleprecision x,jl
integer L
L=5
x=1.D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
x=1.00000001D-3
call BJL(L,x,jl)
print *,'x=',x,' j_5(x)=',jl
end
Then you will see the problem.
-
- Posts: 12
- Joined: April 22 2006
- Affiliation: Sun Yat-Sen University
- Contact:
bug report for CMBFAST and CAMB
Well, struggling with the difference between different computers.
Here is a way to let computer automatically find the divergent region.
program main
real*8 x,jl,y
integer L,i
L=5
y=1.d-8
do i=1,100000000
x=1.d-3+y
call BJL(L,x,jl)
y=y+1.d-8
if(dabs(jl).gt.1)then
print *,'x=',x,' j_5(x)=',jl
pause 'see the problem'
endif
enddo
end program
Then you can see the problem exists for any computer..
Here is a way to let computer automatically find the divergent region.
program main
real*8 x,jl,y
integer L,i
L=5
y=1.d-8
do i=1,100000000
x=1.d-3+y
call BJL(L,x,jl)
y=y+1.d-8
if(dabs(jl).gt.1)then
print *,'x=',x,' j_5(x)=',jl
pause 'see the problem'
endif
enddo
end program
Then you can see the problem exists for any computer..
-
- Posts: 1943
- Joined: September 23 2004
- Affiliation: University of Sussex
- Contact:
Re: bug report for CMBFAST and CAMB
Thanks for the clarification, I see what you mean. In fact I think I've seen this before, which is why CAMB never calls bjl for arguments when the result should be very small - see subroutine GenerateBessels. (I have also tested CAMB by replacing bjl by an essentially exact but slow routine, so, whilst nice to know, I don't think the fix should have any significant effect on CAMB's output). For completeness I will include it in the next update though.
bug report for CMBFAST and CAMB
I ran into and reported this bug for CMBFAST two or three years ago. I don't think anything ever came of it. As far as I can tell, it's not something you'd ever run into during the normal course of running CMBFAST since it never calls bjl with values in that range. I only ran into it because I was doing some special alternative lensing calculations which called bjl directly for many different values.
For [tex]\ell < 6[/tex], bjl calculates the full Bessel function exactly (it's a series of sines and cosines), unless the argument is small. The problem lies in that it doesn't change what it considers "small" as \ell goes from 2 to 5, and depending on how good your machine's trig functions are, you may get bad values. In any case, I made the following diff patch (to cmbfast version 4.5.1, but I doubt jlgen has changed) which fixed things on all machines I tested it on:
For [tex]\ell < 6[/tex], bjl calculates the full Bessel function exactly (it's a series of sines and cosines), unless the argument is small. The problem lies in that it doesn't change what it considers "small" as \ell goes from 2 to 5, and depending on how good your machine's trig functions are, you may get bad values. In any case, I made the following diff patch (to cmbfast version 4.5.1, but I doubt jlgen has changed) which fixed things on all machines I tested it on:
Code: Select all
--- cmbfast4.5.1/jlgen.F Tue Apr 27 19:20:32 2004
+++ cmbfast4.5.1b/jlgen.F Tue Apr 27 19:47:29 2004
@@ -237,7 +237,7 @@
endif
if(l .eq. 3) then
- if(ax .gt. 0.001) then
+ if(ax .gt. 0.02) then
jl=real((cx*(1.0d0-15.0d0/ax2)
1 -sx*(6.0d0-15.0d0/ax2)/ax)/ax)
else
@@ -247,7 +247,7 @@
endif
if(l .eq. 4) then
- if(ax .gt. 0.001) then
+ if(ax .gt. 0.1) then
jl=real((sx*(1.0d0-45.0d0/(ax*ax)+105.0d0
1 /(ax*ax*ax*ax)) +cx*(10.0d0-105.0d0/(ax*ax))/ax)/ax)
else
@@ -258,7 +258,7 @@
if(l .eq. 5) then
- if(ax .gt. 0.001) then
+ if(ax .gt. 0.2) then
jl=real((sx*(15.0d0-420.0d0/(ax*ax)+945.0d0
1 /(ax*ax*ax*ax))/ax -cx*(1.0-105.0d0/(ax*ax)+945.0d0
1 /(ax*ax*ax*ax)))/ax)