How to be more code efficient on fortran? I solved an exercise and I'm not satisfied with my own answer

Question

I'm studing fortran on my own for a job, I'm very new at it. I tried the following exercise and got a correct answer. However I believe there must be many more processing efficient ways to solve the same problem.

The Exercise: Develop a program in Fortran 90 that computes and provides an array with the prime numbers between 2 and 1000.

My solution:

        subroutine locate_primes(prime_array,count)
            integer(4), parameter               :: n = 1000
            integer                             :: i,j
            integer, allocatable, intent(out)   :: prime_array(:)
            integer, intent(out)                :: count
        
            j = 2
            i = 3
            count = 1
            
            do while (i<n+1) !counts how many prime numbers are there on the set
                do while (j<i)
                    if (mod(i,j) == 0 ) then
                        go to 100
                    end if
                    j = j + 1
                    if (j == i) then
                        count = count + 1
                    end if
                end do
100             i = i+1
                j = 2
            end do
            
            allocate(prime_array(count)) !allocate an array of the right size for the prime numbers
            count = 1
            prime_array(1) = 2
            j = 2
            i = 3

            do while (i<n+1) !add create an array with all the prime numbers
                do while (j<i)
                    if (mod(i,j) == 0 ) then
                        go to 101
                    end if
                    j = j + 1
                    if (j == i) then
                        count = count + 1
                        prime_array(count)=i
                    end if
                end do
101             i = i+1
                j = 2
            end do
        end subroutine

My two main problems with my own solution are that I need to check each number between 3 and 1000 twice making it processing inefficient and memory inefficient.

The idea was to create an array of rank one and size one with the first prime (number 2). Then, I would check if the number 3 is prime, and if it is, I would increase the size of the array to two and add the number 3 to the second position. I planned to repeat this process for every number up to 1000.

However, I quickly learned that this approach was not feasible because you cannot simply change the size of an array in Fortran. The only option is to deallocate it, which would result in losing the numbers already saved in it. So, here's what I tried instead:

I thought the best approach would be to have two arrays: one allocatable and one with rank one and a size of 500. Since I know that half of the numbers between 2 and 1000 can't be prime (as they are multiples of two), a size of 500 is definitely enough to store all prime numbers between 1 and 1000. I would save the prime numbers in the size 500 array and keep a count in an integer variable. After that, I could allocate the new array with the size I counted and then populate it with the numbers I had in the size 500 array. Finally, I would deallocate the size 500 array to save memory.

This was my new solution:

        subroutine locate_primes_2(prime_array_final,count)
            integer(4), parameter               :: n = 1000
            integer                             :: i,j
            integer, allocatable, intent(out)   :: prime_array_final(:)
            integer, allocatable                :: prime_array_500(:)
            integer, intent(out)                :: count
        
            allocate(prime_array_500(n/2)) !creates an temporary array to save the prime numbers
            j = 2
            i = 3
            count = 1
            prime_array_500(1) = 2
            
            do while (i<n+1) !counts how many prime numbers there are on the set and saves them
                do while (j<i)
                    if (mod(i,j) == 0 ) then
                        go to 102
                    end if
                    j = j + 1
                    if (j == i) then
                        count = count + 1
                        prime_array_500(count) = i
                    end if
                end do
102             i = i+1
                j = 2
            end do
            
            allocate(prime_array_final(count))
            prime_array_final = prime_array_500(1:count) !saves the set of the prime numbers
            deallocate(prime_array_500) !deallocate the array that ocupied unnecessary memory
            
        end subroutine

It seems to be more processing efficient since I don't need to go through every number between 3 and 1000 twice, but I still end up with a large set that occupies unnecessary memory space before deallocating it. Although I believe I have a better solution than the one I started with, I still wanted to ask the community how one should approach this kind of problem.

I understand that with such a small sample size it may not be a problem but I'm trying to learn how to program more efficiently.

I haven't studied pointers in Fortran yet, but it's my next topic. Feel free to discuss it.

How can I approach processing efficiency and memory efficiency problems in the future? Are there references I can study for this? How can I determine if a new solution is truly better than an old one? Do you have a better solution than mine?

Answer 1

Your question is more about algorithms more than about programing. As others mentioned, look around the "sieve of Eratosthenes", I won't discuss more about algorithms. Regarding programing:

• Don't use goto's when there are simple alternatives. In your exemple use the exit statement instead to exit from the loop. There is also the cycle statement, which skips the rest of the code in the loop and starts the next iteration.
• do while <expression> is deprecated (*) Instead use a infinite loop with an exit test:

do
    if (.not.<expression>) exit
    ...
end do

• Your first do while loop could be replaced by a classical loop since the upper bound does not change during iterations: do i = 1, n
• 1D allocatable arrays can be dynamically extended, for instance to append an element x to an existing array A(:) you can write A = [A, x]. This is the "allocate on assignment" feature. Beware, there are some gotchas here, but it's convenient when you don't know the final size beforehand (although it is also generally not efficient performance wise).
• Fortran pointers are quite different from the C pointers, they are strongly typed and fully describe the pointed objects. Whenever possible, allocatables should be used instead of pointers, but still, there are cases where pointers are useful.

Last remark: Stack Overflow is maybe not the best place when learning a language, because extended discussions are sometimes needed. Therefore I raise your attention to this forum, which is a kind of "unofficial official" forum of the Fortran community, and also to the associated site for some resources.

(*) edit: it's not, my mistake, it's just the recommendation from the authors of a well known reference book... btw there's nothing like "deprecated features" in the official standard, but rather "obsolescent features" (generally features that are meant to be deleted in future revisions).

Answer 2

Well, FWIW, here's a sieve implementation.

A few Fortran things to note:

You can allocate arrays to exactly the required size at run time: allocate( primes( num ) )

You can use count as a useful array function to avoid writing loops.

You can use array-slicing notation is_prime( i * i : max_size : 2 * i ) = .false. to avoid writing loops.

You can (since a Fortran standard that I can't remember) return an allocatable array from a function.

Some books:

For learning: Chapman, S.J., 2017, Fortran For Scientists and Engineers (4th Ed.), McGraw-Hill (but you might need a new mortgage for the 4th edition; the 3rd edition is cheaper)

For reference: Metcalf, M., Reid, J. and Cohen, M., 2018, Modern Fortran Explained, OUP (there's a 2023 version in the pipeline)

program test
   implicit none
   integer, allocatable :: p(:)

   p = sieve( 1000 )
   print "( 'Number of primes = ', i0 )", size( p )
   print "( *( i0, 1x ) )", p

contains
   function sieve( max_size ) result( primes )
      integer, allocatable :: primes(:)
      integer, intent(in) :: max_size
      logical is_prime( max_size )
      integer i, imx, p, num

      ! Set up sieve and remove 1 and all even numbers > 2
      is_prime = .true.;   is_prime(1) = .false.;   is_prime( 4 : max_size : 2 ) = .false.
      ! Knock out odd composite numbers
      imx = sqrt( max_size + 0.0 )
      do i = 3, imx, 2
         if ( is_prime(i) ) is_prime( i * i : max_size : 2 * i ) = .false.
      end do

      ! Count primes (the remaining holes in the sieve)
      num = count( is_prime )

      ! Allocate a prime-containing array of exactly the right size
      allocate( primes( num ) )
      p = 1
      do i = 1, max_size
         if ( is_prime(i) ) then
            primes(p) = i
            p = p + 1
         end if
      end do

   end function sieve

end program test

Output:

Number of primes = 168
2 3 5 7 11 13 17 19 ... lots of numbers ... 971 977 983 991 997 

It's possible to shorten the final collecting-up of primes with the pack() function. Well, short, but a little obscure. Note that this does an automatic allocation of the primes(:) array.

program test
   implicit none
   integer, allocatable :: p(:)

   p = sieve( 1000 )
   print "( 'Number of primes = ', i0 )", size( p )
   print "( *( i0, 1x ) )", p

contains
   function sieve( max_size ) result( primes )
      integer, allocatable :: primes(:)
      integer, intent(in) :: max_size
      logical is_prime( max_size )
      integer i, imx

      ! Set up sieve and remove 1 and all even numbers > 2
      is_prime = .true.;   is_prime(1) = .false.;   is_prime( 4 : max_size : 2 ) = .false.
      ! Knock out odd composite numbers
      imx = sqrt( max_size + 0.0 )
      do i = 3, imx, 2
         if ( is_prime(i) ) is_prime( i * i : max_size : 2 * i ) = .false.
      end do

      primes = pack( [(i,i=1,max_size)], is_prime )

   end function sieve

end program test