my Fortran sieve slows dramatically when array is larger than 800,000,000

Question

I'm a physicist, and I've been doing work with Fortran a lot recently. Originally i used Java extensively for recreation because it was the first language I learned, but I've abandoned it for Fortran and C++. I have an amateur passion for prime numbers and so I created a prime number sieve. I was able to find all the prime numbers up to 2^31 in 15 seconds. This is the maximum array size for Java so that was the end of that. I ported the code carefully( I mean CAREFULLY, I was so upset my code was slow and that I couldn't find an error that I ported my Fortran code back into Java to verify that it wasn't my fault and then ported it back to Fortran, erasing each iteration!). The problem is that around 800,000,000 Fortran will grind to a halt. Up to that point its beating Java, but after that it's insanely slow. I spent a few hours graphing it and fitting a curve. It slows exponentially and could take hundreds of years to solve up to Javas level. I've asked many people to no avail. Why is this happening to me?!?! Are there any wise Fortran coders who can help me? I'm running a Macbook Pro late 2013 i5. My code is below.

program sieve
integer(4),allocatable:: prime(:)
integer(4)::a,max,b,primeCount
write(*,*)"Welcome to the slow prime number sieve!"
write(*,*)"--------------------------------------------"
write(*,*)"Up to what numbers do you need to find primes for?"
write(*,*)"Enter a number below 2^(32-1)"
read*, max

primeCount=0
allocate(prime(max))
prime(1)=1

    do a=2,int(sqrt(real(max))) !the main loop
        if(prime(a)==0)then !if the number is marked as prime procede
            do b=2*a,max,a  !eliminate all the numbers that are multiples of the number
                if(prime(b)==0)then !but only spend time eliminating if the number is still marked prime
                    prime(b)=1
                end if
            end do
        end if
    end do

do a=1,max
    if(prime(a)==0)then
        primeCount=primeCount+1
    end if
end do

print*, primeCount

end program

EDIT: My machine has 4 Gigs of ram installed

Answer 1

I can see several things you could do to speed the code up although none of them seems likely to explain the sharp drop in performance you experience. The most likely culprit seems to be the RAM constraint as Alexander Vogt suggests.

The first thing you should do is to change prime from integer to logical array. This reduces the memory requirements and also speeds up the evaluations of if (prime(a)==0).

The relevant parts of the code would then look as follows

logical(1),allocatable:: prime(:)

primeCount=0
allocate(prime(max))
prime = .false.
prime(1)=.true.

    do a=2,int(sqrt(real(max))) !the main loop
        if(.not. prime(a))then !if the number is marked as prime procede
            do b=2*a,max,a  !eliminate all the numbers that are multiples of the number
                if(.not. prime(b))then !but only spend time eliminating if the number is still marked prime
                    prime(b)=.true.
                end if
            end do
        end if
    end do

do a=1,max
    if(.not. prime(a))then
        primeCount=primeCount+1
    end if
end do

I don't do any Java programming but in Matlab if you declare prime(1:max)=0 and then only switch values bewteen 0 and 1 I would think that Matlab treats the array as a logical array. Java might be doing the same. This could explain why your Java code does not suffer from the performance degradation (assuming the RAM constraint is really the issue).

EDIT: I have done a few experiments.

On my machine (Debian Linux 64bit, i3, 16GB RAM, ifort 14 with default flags) it took 22 seconds for max=800 million (8E8). It took 60 seconds for max=2E9. This is not hours as reported in the question. Also in each case the prime array happened to be initialized to zero.

Moreover using integer(1) makes the programme run 33% faster then using integer(4). With logical(1) it runs less then 5% faster then with integer(1). This behaviour is probably due to the better use of cash as each element of prime occupies less space in memory therefore the processor can do a larger number of iterations with the data currently in cash going through the loops much faster.

The conclusions I would draw from this would be that the culprit was the lack of RAM as noted by Alexander Vogt and that there is a fair chance the author's experience was not affected by the omission to initialize the prime array (although that definitely should not have happened) as HighPerformanceMark pointed out. Further I suspect that Java declared prime as a logical array and that is why the problem did not occur there. (Although 2^31 in 15 seconds in Java? The Fortran code used here comes nowhere close to this. Was really the same code compared?)

Answer 2

This is more of an extended version of my earlier comment than an answer. I think that your mistake in not setting all elements of prime to 0 is affecting the run time of the code, though perversely my initial assessment is that it is making it faster.

Your original code fails to set to 0 all the elements of primes. It is highly likely that when these statements are first encountered:

do a=2,int(sqrt(real(max))) !the main loop
    if(prime(a)==0)then !if the number is marked as prime procede

the value of prime(2) is not 0 so your algorithm fails to mark 2 and its multiples as prime. And it gets worse ! Unless the elements of prime are initialised to 0 then prime(a) cannot ever be guaranteed to be 0 and your program could run to completion without ever marking a number as prime. I'd expect this mistake to make your code faster since it will not enter the inner loop other than by chance.

Faster perhaps but so broken it's not worth measuring its performance.