Coding Fortran: a sorting module

In the previous post we looked at arrays, and a program to calculate the median of a series of numbers. It used the bubblesort() function contained in a module sorting.f95. Here is what the sorting module looks like:

module sorting

contains

subroutine bubblesort(td,A,n)
   integer, intent(in) :: td
   integer, intent(in out), dimension(td) :: A
   integer, intent(in) :: n
   integer :: i, j
   logical :: swapped

   do j = n-1, 1, -1
      swapped = .false.
      do i = 1, j
         if (A(i) > A(i+1)) then
            temp = A(i)
            A(i) = A(i+1)
            A(i+1) = temp
            swapped = .true.
         end if
      end do
      if (.not. swapped) exit
   end do
end subroutine bubblesort

end module sorting

It contains one sorting algorithm at the moment, the bubblesort(), but it would be easy to add others.

 

 

Advertisements

Coding Fortran: Arrays

Having covered the basics of Fortran, we move on to the core data structure found in most languages: the array. I have covered this is some depth in a number of posts on Fortran arrays, including “PLAYING WITH ARRAYS IN FORTRAN” parts (i) and (ii). However let’s consider a sample program. This is a simple program which calculates the median of a series of numbers input by the user. Note that this program uses the module sorting, which contains a bubblesort() function (check out an upcoming post). The program has two subroutines, inputNumbers() which allows the user to input a series of integers, and calcMedian() which calculates the median of the numbers.

program median
   use sorting
   implicit none

   integer :: n=100, num
   integer, dimension(100) :: theArray
   real :: m

   call inputNumbers(n,theArray,num)
   m = calcMedian(n,theArray,num)
   write(*,*) 'The median is ', m

contains

subroutine inputNumbers(aD,A,num)

   integer, intent(in) :: aD
   integer, dimension(aD), intent(in out) :: A
   integer, intent(out) :: num
   integer :: i

   write (*,*) 'How many numbers? '
   read (*,*) num

   do i = 1, num
      read (*,*) A(i)
   end do

end subroutine inputNumbers

real function calcMedian(aD,A,num)
   integer, intent(in) :: aD
   integer, dimension(aD), intent(in out) :: A
   integer, intent(in) :: num

   call bubblesort(aD,A,num)

   if (mod(num,2) == 0) then
      calcMedian = (A(num/2+1) + A(num/2))/2.0
   else
      calcMedian = A(num/2+1)
   end if

end function calcMedian

end program median

In the main program, the array is declared in this manner:

integer, dimension(100) :: theArray

Making it an integer array with a dimension of 100. Arrays in Fortran are extremely flexible, and are capable of whole-array processing and slicing [1]. Note that manipulating elements in an array requires the use of parentheses, “(” and “)”, not square brackets as in some languages (this might take some getting use to… I still make the occasional mistake).

The empty array is then passed to the subroutine  inputNumbers(), which assigns a user specified series of numbers to it. Note the header for the subroutine:

subroutine inputNumbers(aD,A,num) 
   integer, intent(in) :: aD 
   integer, dimension(aD), intent(in out) :: A 
   integer, intent(out) :: num

There are three parameters, illustrating all three ways of passing parameters: in, out, and in out. The first parameter, aD,  indicates the size of the array, and is used when specifying the dimension of the array in the specification of the array parameter A. The array is specified as an “in out” parameter. The last parameter, num, specifies the number of integers input by the user, and stored in the array.

The main program then executes the function calcMedian(), which has a similar set of parameters. Even though calcMedian() does not alter the contents of the array A directly, the subroutine bubblesort() does, so A must be passed to  calcMedian() as “in out“.

The program can be compiled with the module sorting.f95 like this:

gfortran median.f95 sorting.f95

Check out bubblesort() in the next post!

[1] For more info on arrays, see the posts “Playing with Arrays in Fortran” part (i) and (ii).

 

 

Fortran re-engineering: Prime Number Generator (v)

The final fix-ups to the code are somewhat cosmetic. The first of these involves inlining the calculation l=maybe/iposs into the if statement. As what is trying to be calculated is the modular (i.e. remainder), it would be simpler just to use that. So the code:

l = maybe/iposs 
if (l*iposs-maybe == 0) then

Instead becomes:

if (mod(maybe,iposs) == 0) then

This means it calculates the remainder of maybe/iposs. If this is zero, then iposs is a factor of maybe, and maybe cannot be a prime number. This also eliminates the variable l. Finally, we simplify the formatting in the print 30 statement, through inline embedding. So:

print 30, maybe 
30 format(1h,10x,i7) 

becomes:

print '(1h,10x,i7)', maybe

Finally, we can remove the stop, and maybe change a couple of the variables to something a little more meaningful: mayber to maybeSQ, and iposs to poss.

Here is the final program:

! prime number generator
program png
   implicit none
   integer :: maybe, maybeSQ, poss
   logical :: isPrime

   print *, "Prime numbers from 3 to 10,000,000"

   do maybe = 3,9999997,2
      isPrime = .true.
      maybeSQ = sqrt(float(maybe))
      do poss = 2,maybeSQ
         if (mod(maybe,poss) == 0) then
            isPrime = .false.
            exit
         end if
      end do
      if (isPrime .eqv. .true.) then
         print '(1h,10x,i7)', maybe
      end if
   end do
end program png

This program makes much more sense, and we have effectively eliminated all “jump” labels. The program as it stands is pretty much re-engineered.

 

Fortran re-engineering: Prime Number Generator (iv)

Looking at the code as it stands (below), you will notice that the lower arithmetic if statement calculates the difference between the current value of maybe, and 9999997, meaning that while maybe has the values 3→9999997, it jumps up to label 10, and starts the process again with the next value of maybe. When its value hits 9999999, the program terminates. This is the classic type of legacy code which indicates a global loop.

!     prime number generator
      program png
         implicit none
         integer :: maybe, mayber, iposs, l

      print *, "Prime numbers from 3 to 10,000,000"

      maybe = 1
   10 maybe = maybe + 2
      mayber = sqrt(float(maybe))
      do iposs = 2,mayber
         l = maybe/iposs
         if (l*iposs-maybe == 0) then
            goto 40
         end if
      end do
      print 30, maybe
   30 format(1h,10x,i7)
   40 if (maybe-9999997) 10,10,50
   50 stop
      end program png

Converting this code, just requires the addition of an appropriate do loop, and the removal of the arithmetic if. Consequently, we are also able to dispose of the code at label 10, and the line of code before that (and the label 50, which is now superfluous). The do loop now starts at 3, ends at 9999997, and increments by 2. The arithmetic if at label 40 is also gone, but notice we have had to add a cycle keyword (similar to continue in C), to deal with the goto 40 from within the loop. Fortran won’t allow for empty statements next to labels, and the goto 40 actually jumps out of the inner loop, and bypasses the print 30, which activates if a number is a prime. If we used exit to exit the loop, it would still print the number, even if it was not a prime.

!     prime number generator
      program png
         implicit none
         integer :: maybe, mayber, iposs, l

      print *, "Prime numbers from 3 to 10,000,000"

      do maybe = 3,9999997,2
         mayber = sqrt(float(maybe))
         do iposs = 2,mayber
            l = maybe/iposs
            if (l*iposs-maybe == 0) then
               goto 40
            end if
         end do
         print 30, maybe
   30    format(1h,10x,i7)
   40    cycle
      end do
      stop
      end program png

This is actually a good case for using goto to jump to a place in a program, bypassing code.

However there is a better solution, which does not use a goto, instead it uses a logical variable to control flow in the program. In this case, we use a variable isPrime, which is toggled to true at the start of testing each value from 3 to 9999997. If within the inner loop, the number being tested is found to be divisible by some factor, isPrime is set to false, and the inner loop exits. Printing only occurs if the inner loop runs its course, and isPrime is still true (i.e. the number *is* a prime).

!     prime number generator
      program png
         implicit none
         integer :: maybe, mayber, iposs, l
         logical :: isPrime

      print *, "Prime numbers from 3 to 10,000,000"

      do maybe = 3,9999997,2
         isPrime = .true.
         mayber = sqrt(float(maybe))
         do iposs = 2,mayber
            l = maybe/iposs
            if (l*iposs-maybe == 0) then
               isPrime = .false.
               exit
            end if
         end do
         if (isPrime .eqv. .true.) then
            print 30, maybe
            30 format(1h,10x,i7)
         end if
      end do
      stop
      end program png

This program makes much more sense, and we have effectively eliminated all “jump” labels. The program as it stands is pretty much re-engineered. In the final post in this series, we will look at tidying up some loose ends.

 

Coding Fortran: functions

Now let’s convert the program to deal with Fortran functions, which are subprograms that return a value. In this case we will have the function fact(), only calculate a factorial value.  Functions return a single value. Here is what the program looks like with the addition of a function.

program factorial
   implicit none
   integer :: num

   write (*,*) "n! ? "
   read(*,*) num

   write(*,'(I3,"! = ", I20)') num, fact(num)

contains

integer (kind=16) function fact(n)

   implicit none
   integer, intent(in) :: n
   integer :: i
   integer (kind=16) :: f

   f = 1
   do i = 2, n
      f = f * i
   end do

   fact = f

end function fact

end program factorial

LEGEND: function wrapper, function parameter, function code, function call

The function is situated within the program using the keyword contains.  It has one parameter, n, which is declared on the line after the function declaration intro (with intent in). The remainder of the code is as it would appear in the previous program. At the end of the function, the value calculated in the variable f is assigned to the name of the function, i.e. fact = f. If done this way, no explicit return statement is needed. The function is called from the main program in the same way it is called in most languages.

Functions in this manner are different to subroutines, as they are declared internally. It is also possible to declare then externally, either external to the main program in the same file or in a different file. For example, they can be declared external in the same file in the following manner:

program factorial
   implicit none
   integer :: num
   integer (kind=16), external :: fact

   write (*,*) "n! ? "
   read(*,*) num

   write(*,'(I3,"! = ", I20)') num, fact(num)

end program factorial

integer (kind=16) function fact(n)

   implicit none
   integer, intent(in) :: n
   integer :: i
   integer (kind=16) :: f

   f = 1
   do i = 2, n
      f = f * i
   end do
   fact = f

end function fact

Notice that the function is declared, using appropriate return data types in the main program. Moving the function to an external file is then as simple as making a file fact.f95 containing only the function fact(). It is compiled in the following manner:

gfortran factFunction3.f95 fact.f95

 

 

 

 

Fortran re-engineering: Prime Number Generator (iii)

The next step involves choosing a path. In this case we will opt to modify the easy things first, and the easiest component is the first if statement. This arithmetic if statement represents a form of 3-way decision statement, depending on whether the expression in the parentheses evaluates to negative, 0, or positive. The basic form of the if statement in this case is:

if (l*iposs-maybe) 20,40,50

The value of l represents whether or not a number (maybe) is divisible by iposs. If maybe=7, and iposs=3, then l=2 (so not divisible without a remainder). So the expression in the if statement basically branches to label 20 if the l*iposs-maybe is negative, i.e. there is a remainder, and label 40 if the value is zero, i.e. divisible. The label 50 is superfluous (because a third possibility had to be provided).

This can be converted to a contemporary if statement by modifying the expression, so that rather than perform a calculation, it performs a logical comparison. For example, if (l*iposs-maybe == 0) is true, then the program branches to label 40, otherwise do nothing, in which case, it drops into the next iteration of the inner loop to try the next value of iposs. Here is what the new if statement looks like:

!     prime number generator
      program png
         implicit none
         integer :: maybe, mayber, iposs, l

      print *, "Prime numbers from 3 to 10,000,000"

      maybe = 1
   10 maybe = maybe + 2
      mayber = sqrt(float(maybe))
      do 20 iposs = 2,mayber
      l = maybe/iposs
      if (l*iposs-maybe == 0) then
         goto 40
      end if
   20 continue
      print 30, maybe
   30 format(1h,10x,i7)
   40 if (maybe-9999997) 10,10,50
   50 stop
      end program png

Now we can convert the do-continue loop to a contemporary “do-end do” loop, eliminating the label 20 in the process. This is achieved by removed the label 20 from the start of the loop, and replacing the 20 continue at the end of the loop with end do. We also perform a little housekeeping indentation.

!     prime number generator
      program png
         implicit none
         integer :: maybe, mayber, iposs, l

      print *, "Prime numbers from 3 to 10,000,000"

      maybe = 1
   10 maybe = maybe + 2
      mayber = sqrt(float(maybe))
      do iposs = 2,mayber
         l = maybe/iposs
         if (l*iposs-maybe == 0) then
            goto 40
         end if
      end do
      print 30, maybe
   30 format(1h,10x,i7)
   40 if (maybe-9999997) 10,10,50
   50 stop
      end program png

The program now has a loop that checks to see if a potential prime number is divisible by numbers less than its square root. If a number is divisible, then it jumps to label 40, and tests the next potential prime number. The program is starting to shed some of its legacy features, and become more structured. We have added an implicit goto statement, and there is the remaining arithmetic if to contemplate.

The next step is to turn the remaining arithmetic if statement into a global loop.

 

Fortran re-engineering: Prime Number Generator (ii)

The first thing to do after analyzing the code is to perform some housekeeping. First thing? – convert the code to lowercase, and convert the comment delimiters from C to !. Now the program will compile as a “contemporary” Fortran 95 program using a .f95 extension (even though it is legacy code, it compiles due to Fortran’s wonderful backwards compatibility).  Now we can start the process of re-engineering it.

!     prime number generator
      program png

      print 5
    5 format(36h1prime numbers from 3 to 10,000,000 ,///)
      maybe = 1
   10 maybe = maybe + 2
      mayber = sqrt(float(maybe))
      do 20 iposs = 2,mayber
      l=maybe/iposs
      if (l*iposs-maybe) 20,40,50
   20 continue
      print 30, maybe
   30 format(1h,10x,i7)
   40 if (maybe-9999997) 10,10,50
   50 stop
      end program png

Now, we will remove the implicit variables. This is achieved by adding the clause “implicit none“, after which the variables used in the program will have to be explicitily declared. With implicit variables, any variables used starting with i, j, k, l, m, or n are automatically declared as integers. The format statement associated with label 5 is integrated directly into the preceding print statement. Making changes liking this is more along the lines of refactoring, making the code cleaner. Here is the program after these changes.

! prime number generator 
      program png
      implicit none
      integer :: maybe, mayber, iposs, l

      print *, "Prime numbers from 3 to 10,000,000"

      maybe = 1
   10 maybe = maybe + 2
      mayber = sqrt(float(maybe))
      do 20 iposs = 2,mayber
      l = maybe/iposs
      if (l*iposs-maybe) 20,40,50
   20 continue
      print 30, maybe
   30 format(1h,10x,i7)
   40 if (maybe-9999997) 10,10,50
   50 stop
      end program png

At this stage one has to decide what path to take when re-engineering the code. This is a small piece of code, and yet the unstructured components of it make it more challenging. Imagine an unstructured program with 1,000 LOC? 10,000 LOC?

Fortran re-engineering: Prime Number Generator (i)

This is an example of re-engineering a small Fortran program. This Fortran program, written in some version of Fortran pre-F90, computes all prime numbers from 3 to 10,000,000. The program calculates a prime number by taking a number and successively dividing it by all numbers less than the square root of the number. Here is what the legacy program looks like:

C     PRIME NUMBER GENERATOR
      PRINT 5
    5 FORMAT(36H1PRIME NUMBERS FROM 3 TO 10,000,000 ,///)
      MAYBE = 1
   10 MAYBE = MAYBE + 2
      MAYBER = SQRT(FLOAT(MAYBE))
      DO 20 IPOSS = 2,MAYBER
      L=MAYBE/IPOSS
      IF (L*IPOSS-MAYBE) 20,40,50
   20 CONTINUE
      PRINT 30, MAYBE
   30 FORMAT(1H,10X,I7)
   40 IF (MAYBE-9999997) 10,10,50
   50 STOP
      END

This may seem like a small program, at 15 LOC, however it is sometimes best to start small. The first task at hand, apart from converting it to lowercase involves read the code, and analyzing the syntax. Now, because Fortran compilers like gfortran are backwards-compatible, this code actually compiles and runs.

While the code does not contain any explicit goto statements, it does contain two arithmetic if statements, which are essentially early forms of if statements which offer three alternatives. Kind-of de facto goto statements. Now let’s look at the code if we print it out and do some thinking about how the code flows. I like to do this by annotating the printed code. If possible I use different colours to denote different labels. In this case I am trying to get an idea of how the program flows.

The first thing to do is find out what the variables do.

  • MAYBE signifies the odd values starting at 2, and ending at 9,999,997.
  • MAYBER is the square root of MAYBE.
  • IPOSS represents the prime dividers for each number, 2→MAYBER
  • L is the value calculated when MAYBE is divided by IPOSS. Used in the IF statement that follows to determine if IPOSS divides into MAYBE exactly.

In annotating the code, I have used square boxes, and a highlighter to mark the occurrence of the labels. Looking first at the last IF statement, we notice that it is essentially a two-way branch, forming an outer loop. The inner IF statement controls the DO loop. Analyzing the program may allow you to create a flow-chart of some sort. Here is a flowchart which represents what happens (with I=MAYBE, and K=IPOSS).

 

Coding Fortran: subroutines

In the previous example we looked at loops in Fortran. Now we will convert the Factorial program to deal with subprograms. In Fortran there are essentially two types of subprograms: functions, which return a single value, and subroutines, which are called and can have any number of inputs and outputs. We’ll deal with subroutines first.

Here is the program sans subprograms:

program factorial

   integer :: i, num
   integer (kind=16) :: fact

   write (*,*) "n! ? "
   read(*,*) num

   fact = 1
   do i = 2, num
      fact = fact * i
      write(*,'(I3,"! = ", I20)') i, fact
   end do

end program factorial

Converting the actual code that performs the calculation into a subroutine is relatively simple. Here is what the program looks like with the addition of a subroutine.

program factorial
   implicit none
   integer :: i, num

   write (*,*) "n! ? "
   read(*,*) num
  
   call fact(num)

end program factorial

subroutine fact(n)
   integer, intent(in) :: n
   integer :: i
   integer (kind=16) :: f

   f = 1
   do i = 2, n
      f = f * i
      write(*,'(I3,"! = ", I20)') i, f
   end do

end subroutine fact

LEGEND: subroutine wrapper, subroutine parameter, subroutine code, subroutine call

The subroutine is situated at the end of the actual program that calls it (the subroutine can also be included inside the main program or in an external file). It has one parameter, n, which is declared on the line after the subroutine declaration intro.  It is declare with intent in, which means that it is really a “pass-by-value”, variable with a value that is coming in to the subroutine. Other options are out, and inout. The remainder of the code is as it would appear in the previous program.

The subroutine is called using the keyword call (which also differentiates it from a function).

Note that if we declare the subroutine fact() inside the main program, it looks like this:

program factorial
   implicit none
   integer :: i, num

   write (*,*) "n! ? "
   read(*,*) num

   call fact()

contains

subroutine fact()
   integer :: i
   integer (kind=16) :: f

   f = 1
   do i = 2, num
      f = f * i
      write(*,'(I3,"! = ", I20)') i, f
   end do

end subroutine fact

end program factorial

Note, used in this manner the variable num from the main program can be used directly in the subroutine fact(), because it is accessible via the scope of the program. I wouldn’t necessarily promote this way of doing things, as it can lead to side-effects if a variable is in accidentally modified in the wrong place.

 

 

 

 

Coding Fortran: loops (or rather the do loop)

With “Hello World!” done, let’s move onto something slightly more complex: the factorial. This will require using a datatype, some variables, a loop, and formatted output. Here is the program:

program factorial

   integer :: i, num
   integer (kind=16) :: fact

   write (*,*) "n! ? "
   read(*,*) num

   fact = 1
   do i = 2, num
      fact = fact * i
      write(*,'(I3,"! = ", I20)') i, fact
   end do

end program factorial

You can see that the program has three integer variables: i, num, and fact. Note that for the variable fact, it has an attribute kind = 16. which implies it has a size of 16-bytes. The next portion of the program deals with prompting the user for the number of factorials to calculate, which is stored in num.

Now the core part of the program, which deals with calculating each of the factorials. This is done using a do loop, the Fortran equivalent of the for loop. The variable i is used to control the number of iterations of the loop, in this case from 2 → num. Notice that the for loop uses a control statement terminator, in the form of end do. This encapsulates the loop as a closed block, regardless of the number of lines of code inside the loop. The loop isn’t that different from loops in other languages, apart from the fact that it is relatively simple in structure.

The last piece of interesting code is the write statement. Note that in the first write statement, the formatting sequence is (*,*), implying standard output, no formatting. The write statement in the for loop uses standard output (the first *), but replaces the second * with a formatting statement:

'(I3,"! = ", I20)'

where the following formatting:

I3     → an integer with a width of 3, right-justified
"! = " → a string, printed as written
I20    → an integer with a width of 20, right-justified

Here is what is output from the program.

 n! ?
20
  2! =                   2
  3! =                   6
  4! =                  24
  5! =                 120
  6! =                 720
  7! =                5040
  8! =               40320
  9! =              362880
 10! =             3628800
 11! =            39916800
 12! =           479001600
 13! =          6227020800
 14! =         87178291200
 15! =       1307674368000
 16! =      20922789888000
 17! =     355687428096000
 18! =    6402373705728000
 19! =  121645100408832000
 20! = 2432902008176640000