[Numpy-tickets] [NumPy] #924: problem with summing large array of float32

[NumPy]

#924: problem with summing large array of float32
------------------------+---------------------------------------------------
 Reporter:  emil        |       Owner:  somebody         
     Type:  defect      |      Status:  new              
 Priority:  high        |   Milestone:                   
Component:  numpy.core  |     Version:  1.0.1            
 Severity:  critical    |    Keywords:  sum, dot, float32
------------------------+---------------------------------------------------
 In medical imaging large arrays of single floats are often used. In
 summing over such an
 array I discovered the following problem:
 a=ones([840,2200,60],float32)
 sum(a.flatten()) yields 16777216 (4096^2)   ???

 However, sum(sum(sum(a,axis=0),axis=0),axis=0) yields 110880000 the
 correct result.
 The problem doesn't occur with float64
 dot(a.flatten() , a.flatten() ) also gives 16777216

 I don't think the problem lies only with numpy. I tried to circumvent the
 problem by
 writing my own summation in fortran (gfortran) and using f2py. The
 following fortran snippet
 gave me also incorrect results

        Real*4 sum

        sum = 0.

        do i = 1,840
              do j = 1,2200
                     do k = 1,60
                             sum = sum +1.0
                     end do
               end do
        end do

 sum ends up being 16777216 again! the incorrect result!
 on the other hand the following works:
        Real*4 sum1, sum

        sum = 0.

        do i = 1,840
              sum1 = 0.
              do j = 1,2200
                     do k = 1,60
                             sum1 = sum1 +1.0
                     end do
               end do
              sum = sum+sum1
        end do


 The same problem seems to come up in matlab.

 My computer has amd64's (opteron) and is running suse linux

-- 
Ticket URL: <http://scipy.org/scipy/numpy/ticket/924>
NumPy <http://projects.scipy.org/scipy/numpy>
The fundamental package needed for scientific computing with Python.