bug report for CMBFAST and CAMB

Use of Healpix, camb, CLASS, cosmomc, compilers, etc.
Post Reply

Do you think this bug would cause some "effect" somewhere?

Poll ended at May 06 2006

Yes
1
100%
No
0
No votes
 
Total votes: 1

Zhiqi Huang
Posts: 12
Joined: April 22 2006
Affiliation: Sun Yat-Sen University
Contact:

bug report for CMBFAST and CAMB

Post by Zhiqi Huang » April 22 2006

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.

Antony Lewis
Posts: 1498
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: bug report for CMBFAST and CAMB

Post by Antony Lewis » April 22 2006

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}).

Zhiqi Huang
Posts: 12
Joined: April 22 2006
Affiliation: Sun Yat-Sen University
Contact:

bug report for CMBFAST and CAMB

Post by Zhiqi Huang » April 22 2006

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.

Zhiqi Huang
Posts: 12
Joined: April 22 2006
Affiliation: Sun Yat-Sen University
Contact:

bug report for CMBFAST and CAMB

Post by Zhiqi Huang » April 22 2006

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.

Zhiqi Huang
Posts: 12
Joined: April 22 2006
Affiliation: Sun Yat-Sen University
Contact:

bug report for CMBFAST and CAMB

Post by Zhiqi Huang » April 22 2006

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.

Zhiqi Huang
Posts: 12
Joined: April 22 2006
Affiliation: Sun Yat-Sen University
Contact:

bug report for CMBFAST and CAMB

Post by Zhiqi Huang » April 22 2006

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..

Antony Lewis
Posts: 1498
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: bug report for CMBFAST and CAMB

Post by Antony Lewis » April 22 2006

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.

Ben Gold
Posts: 81
Joined: September 25 2004
Affiliation: University of Minnesota
Contact:

bug report for CMBFAST and CAMB

Post by Ben Gold » April 23 2006

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:

Code: Select all

--- cmbfast4.5.1/jlgen.F        Tue Apr 27 19&#58;20&#58;32 2004
+++ cmbfast4.5.1b/jlgen.F       Tue Apr 27 19&#58;47&#58;29 2004
@@ -237,7 +237,7 @@
             endif
 
             if&#40;l .eq. 3&#41; then
-                if&#40;ax .gt. 0.001&#41; then
+                if&#40;ax .gt. 0.02&#41; then
                 jl=real&#40;&#40;cx*&#40;1.0d0-15.0d0/ax2&#41;
      1                    -sx*&#40;6.0d0-15.0d0/ax2&#41;/ax&#41;/ax&#41;
                 else
@@ -247,7 +247,7 @@
             endif
 
             if&#40;l .eq. 4&#41; then
-                if&#40;ax .gt. 0.001&#41; then
+                if&#40;ax .gt. 0.1&#41; then
         jl=real&#40;&#40;sx*&#40;1.0d0-45.0d0/&#40;ax*ax&#41;+105.0d0
      1       /&#40;ax*ax*ax*ax&#41;&#41; +cx*&#40;10.0d0-105.0d0/&#40;ax*ax&#41;&#41;/ax&#41;/ax&#41;
                 else
@@ -258,7 +258,7 @@
 
              if&#40;l .eq. 5&#41; then
 
-                if&#40;ax .gt. 0.001&#41; then
+                if&#40;ax .gt. 0.2&#41; then
         jl=real&#40;&#40;sx*&#40;15.0d0-420.0d0/&#40;ax*ax&#41;+945.0d0
      1    /&#40;ax*ax*ax*ax&#41;&#41;/ax -cx*&#40;1.0-105.0d0/&#40;ax*ax&#41;+945.0d0
      1                                    /&#40;ax*ax*ax*ax&#41;&#41;&#41;/ax&#41;

Post Reply