0
program analysis
        implicit none
        
        integer i, j, k, l
        double precision a, b, c, d, e
        integer binb(1:20),binc(1:20),bind1(1:20),bine(1:20)
        real lower(1:20),upper(1:20), next
    
        
    
        character(100) event
        
    
        upper(1)=-2.7
        lower(1)=-3.0
        binb(1:20)=0
        binc(1:20)=0
        bind1(1:20)=0
        bine(1:20)=0
        
        
        next=.3
        
    
    
            do l=2,20
                lower(l)=lower(l-1)+next
                upper(l)=upper(l-1)+next
            end do

        
            
        open(unit = 7, file="zpc_initial_momenta.dat")
            
            
            do k=1, 10
                read(7,'(A)') event
                do j=1,4000
                read(7,*) a, b, c, d, e
                
            do i=1, 20
                if(b>=lower(i) .and. b<upper(i)) then
            binb(i)=binb(i)+1
                end if
                if(c>=lower(i) .and. c<upper(i)) then 
            binc(i)=binc(i)+1
                end if
                if(d>=lower(i) .and. d<upper(i)) then
            bind1(i)=bind1(i)+1
                end if
                if(e>=lower(i) .and. e<upper(i)) then 
            bine(i)=bine(i)+1
                end if
            
            
            end do
            end do
            end do
            
            
                
                
        close(7)
        
            
        
        open(unit = 8, file="outputanalysis2.dat")
            Write(8,*) 'The bins in each column are as follows:'
            Write(8,*) 'FIRST COLUMN (MOMENTUM IN X DIRECTION)'
            write(8,*) binb(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(binb)
            Write(8,*) 'SECOND COLUMN (MOMENTUM IN Y DIRECTION)'
            write(8,*) binc(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(binc)
            Write(8,*) 'THIRD COLUMN (MOMENTUM IN Z DIRECTION)'
            write(8,*) bind1(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(bind1)
            Write(8,*) 'FOURTH COLUMN (KINETIC ENERGY)'
            write(8,*) bine(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(bine)
            
            
        close(8)
        
        end program
                

I have an issue that I'm half convinced is due to Fortran's base 2 round-off error, but I'm not entirely sure: In my if statements designed to count data for different ranges, I have an array for the upper and lower bound of each bin. I had a previous program that did not use arrays but was incredibly lengthy. In my array program, when executing the program, some bins are off by a max of 3 counts (not all of them, and some of them are off by less than that or none at all.) I did some further digging and compared it to a declared non-array variable equal to one of the bin's upper and lower bounds, and printed it out. For example: say, bin number 5, whose bin length is from (-1.8 to -1.5) for which it counts values in this range. When I used the non-array upper and lower bound to see if there was a change in the bin count, there was: it was 897, which is the correct bin count (found in the previous nonarray program.) However, when I use arrays to define each upper and lower limit, putting them into the if statement (as seen in this program) the count is 900. To see what was going on, I then had my program print out the nonarray variable equal to -1.8. When executed the program prints out: -1.79999995. When I had my program print out the array lower bound equivalent (lower(5)), it prints out -1.80000019. Can someone explain to me what is going on here, and if that is the reason the counts are off slightly for my array program? Is there a way to fix this? Thank you!

Edit: I think the best way to reproduce this is with a test program:

program help
    implicit none
    
    integer i
    real lower(1:5),next,nonarraylower2
    
    
    
    lower(1)=-2.1
    nonarraylower2=-0.9
    next=0.3
    
    do i=2,5
        lower(i)=lower(i-1)+next
    end do
    
    write(*,*) lower(5),nonarraylower2
    
    end program

If we write lower(2) and set nonarraylower2 equal to the same value that should be there according to the do loop, the values print to be the same. However, do lower(5) and set nonarraylower2 equal to what should be the same value (-.9) and have it print out both and you see there begins to show a difference. I should clarify: the data file I am pulling from has numbers small enough to where these minute changes in decimal places may matter. For example, here are some data points I am working with:

21  0.1314E+00 -0.1232E+01  0.6258E+00  0.1388E+01
    21 -0.3478E+00 -0.1268E+01 -0.3834E+00  0.1370E+01
    21  0.2138E+00  0.1995E+01 -0.9115E-01  0.2009E+01
    21  0.8936E+00 -0.5758E-01  0.7773E+00  0.1186E+01
    21 -0.5949E+00 -0.1999E+00  0.3787E+00  0.7330E+00

I hope my issue is clearer. When I originally built a program without using arrays to read each column and count how many belong to different ranges (defined by discrete bins), the count can sometimes be off by 5 maximum of my new program's output counts (it isn't consistent and sometimes I get the same counts. It seems the error is between 1 and 5 lost or extra counts.) What is going on here, and is there a way to fix it?

1 Answers1

1

I could not reproduce your error as I do not have the data, but your code should be fine unless you need high precision. I suggest you use real*8 for all float variables and parameters and assign values with the same precision: like next=0.3d0. I have uploaded a minimal code that worked for me. Hope that helps.

    program analysis
    implicit none
    
    integer :: i, j, k, l
    real*8 :: a, b, c, d, e
    integer :: binb(1:20),binc(1:20),bind1(1:20),bine(1:20)
    real*8 :: lower(1:20),upper(1:20), next

    integer :: bincheck


    upper(1)=-2.7d0
    lower(1)=-3.0d0
    binb(1:20)=0   
    
    next=0.3d0
    


        do l=2,20
            lower(l)=lower(l-1)+next
            upper(l)=upper(l-1)+next
        end do

    
        
    open(unit = 100, file="zpc_initial_momenta.dat")   
    do j=1,300
    read(100,*) b
            
    do i=1, 20
    if(b>=lower(i) .and. b<upper(i)) then
    binb(i)=binb(i)+1
    end if
        
    end do
    end do

    close(100)
      
!-------- check --------

    do i=1,20
    write(*,*) binb(i)
    enddo  
   
    open(unit = 100, file="zpc_initial_momenta.dat") 
    bincheck = 0
    do i=1,300
    read(100,*) b
    if(b>=lower(1) .and. b<upper(1)) then
    bincheck = bincheck +1
    end if
    enddo

    write(*,*) 'values in bin 1', bincheck


        
            
            
    close(100)
    
    end program

EDIT: example to check if values change when saved in arrays:

    program analysis
    implicit none
    
    integer :: i
    real*8 :: a
    real*8 :: alist(1:20)

    a=-1.79999995d0

    write(*,*) 'nonarray', a

    do i=1,20
    alist(i)=a
    write(*,*) 'array', i, alist(i)
    enddo

    end program

It gives me

 nonarray  -1.7999999499999999     
 array           1  -1.7999999499999999     
 array           2  -1.7999999499999999     
 array           3  -1.7999999499999999     
 array           4  -1.7999999499999999     
 array           5  -1.7999999499999999     
 array           6  -1.7999999499999999 
---- and so on -----

EDIT 2: modified program help

    program help
    implicit none
    
    integer i
    real*8 :: lower(1:5),next,nonarraylower2
    
    lower(1)=-2.1d0
    nonarraylower2=-0.9d0
    next=0.3d0
    
    do i=2,5
        lower(i)=lower(i-1)+next
    end do
    
    write(*,*) lower(5),nonarraylower2
    
    end program
Deb S. B.
  • 78
  • 8
  • I tried doing real*8 before my real values (columns a,b,c,d,e are all double precision numbers in the original file) but the output was the same as it was previously. Do you think it has something to do with how real type arrays store data vs. how real type nonarrays store data? – Garrett Leigh Apr 28 '23 at 04:02
  • I have just edited my post with a way to reproduce my problem. – Garrett Leigh Apr 28 '23 at 04:38
  • It is simply a matter of precision. I ran your version of ```program help``` and here is what I get: ```-0.900000036 -0.899999976``` which is not wrong because your variables do not have double precision. When I change it (EDIT 2 in my answer) I get ```-0.89999999999999991 -0.90000000000000002```. Apart from declaring variables as double precision, it is also very important to assign values that are double precision too like ```1.d0``` not ```1.0```. – Deb S. B. Apr 28 '23 at 04:55
  • So, do you think that would be what is causing the difference in bin counts across both programs? – Garrett Leigh Apr 28 '23 at 04:59
  • I would guess so. Check some of the answers around ```double precision``` like this: https://stackoverflow.com/questions/9211037/extended-double-precision – Deb S. B. Apr 28 '23 at 05:07
  • Thanks! I know what the problem is now. The *original* program without arrays was simply inaccurate, as it did not have the precision that the arrays do with real and real*8. Once I made the nonarray variable a real*8 variable and used it as the upper and lower bounds, the output was the same count of the original array function. – Garrett Leigh Apr 28 '23 at 05:12
  • 1
    `real*8` is not Fortraj. It has never been accepted into any Fortran standard, it is just an extension. Fortran has sufficient means to select precision in a standard way https://stackoverflow.com/questions/838310/fortran-90-kind-parameter AVOID the answer by John Skeet at https://stackoverflow.com/questions/10520819/what-does-real8-mean it is only updated so much because the name of the answerer even if he knows nothing about Fortran. – Vladimir F Героям слава Apr 28 '23 at 05:57