Interoperability : Calling Fortran from C (ii)

Here is a second example, using some pass-by-reference parameters in C. This time, let’s calculate the hypotenuse of a triangle, using a function in Fortran to perform the actual calculation. First, let’s consider the Fortran side of the code. This time, we have just built a single subroutine (file = pythagoras.f95).

subroutine pythagoras(a, b, c) bind(C)

   use iso_c_binding
   implicit none
   real (C_FLOAT), intent(in) :: a, b
   real (C_FLOAT), intent(out) :: c
 
   c = sqrt(a*a + b*b)
  
end subroutine pythagoras

From the perspective of writing the Fortran code, there is very little difference from a normal subroutine, except for the specification of the bind,  use of the iso_c_binding module, and specification of the subroutine parameters as C types equivalents.

From the perspective of the C program (hypotenuse.c), the subroutine pythagoras is cited as being external, and when called, both pass-by-value and pass-by-reference parameters are tagged with &. This is because although there is a difference in C, in Fortran, all parameters are “pass-by-reference”.

#include <stdio.h>

extern void pythagoras(float *a, float *b, float *c);

int main(void)
{
   float hyp, sideA, sideB;
   sideA = 3;
   sideB = 4;
   pythagoras(&sideA, &sideB, &hyp);
   printf("The hypotenuse of triangle with sides");
   printf(" %.2f and %.2f is %.2f\n", sideA, sideB, hyp);

return 0;
}

The programs can then be compiled in the following manner:

gcc -c hypotenuse.c
gfortran pythagoras.f95 hypotenuse.o

When the program (a.out) is executed, here is the output:

The hypotenuse of triangle with sides 3.00 and 4.00 is 5.00

 

 

 

Advertisements

Coding Fortran: strings

As mentioned in a previous post, strings in Fortran are different to character arrays. To illustrate some of the features of strings, we will look at a program which reads in a strings from the keyboard, and strips all the blanks from the string. Instead of having a dimension, like arrays, they are declared using a length, in the case of the example program below len=30.

program stripBlanks

   character (len=30) :: str, strip
   integer :: strL, i, ipos

   write (*,*) 'Enter a sentence: '
   read (*,'(a)') str
   strip = ' '

   strL = len_trim(str)
   ipos = 0
   do i = 1,strL
      if (str(i:i) /= ' ') then
         ipos = ipos + 1
         strip(ipos:ipos) = str(i:i)
      end if
   end do

   write (*,'(a)') str
   write (*,'(a)') strip

end program

Notice that when reading and writing strings, we use the format string ‘(a)’. After the input has been obtained from the user, the length of the string is calculated. This is done using the function len_trim(), which calculates the length of the string minus any trailing blanks. This means the length calculated will not be 30, as it would using the function len(). Also notice the use of a whole array operation, setting the string strip (which will hold the sentence stripped of blanks) to blanks.

Then a loop is used to traverse the string str, and extract any element of the string which is not a blank, assigning it to the string strip. Note that unlike normal arrays, which use one index, because this is a string it uses substrings, and therefore a single index has to be expressed using the substring ipos:ipos.

Here is a sample run of the program:

Enter a sentence:
the cat in the hat came back
the cat in the hat came back
thecatinthehatcameback

Strings also have a series of built-in functions that normal character arrays don’t. Strings in Fortran are actually treated as “collections of characters”, and character arrays are just that – arrays of characters, i.e. they are treated in the same way as arrays of integers or reals.

Interoperability : Calling Fortran from C (i)

In the world of integrated programming languages, sometimes you might want to call a subroutine written in Fortran from C. There are a multitude of Fortran scientific subroutines out there, and there is no inherent need to reinvent the wheel.

Here is a sample Fortran module containing a subroutine called fortransub() :

module fortmodule
use iso_c_binding
implicit none

integer(C_INT), bind(C), dimension(5) :: numbers

contains

subroutine fortransub() bind(C)

   print *, "Hello from Fortran!"
   numbers(1) = 1
   numbers(2) = 2
   numbers(3) = 3
   numbers(4) = 4
   numbers(5) = 5

end subroutine

end module

The Fortran piece of code uses the Fortran ISO C Binding, which became part of the Fortran 2003 language standard. The ISO C Binding, when used in the declaration of a Fortran subroutine or function, causes the Fortran compiler to use the C calling conventions so that that subprogram can be directly called from C.

The function fortransub() merely prints out a string, and assigns values to the integer array numbers, which is global to the module, and hence can be accessed externally. The C program that calls fortransub() looks like this:

#include <stdio.h>

extern void fortransub();
extern int numbers[5];

int main(void)
{
   int i;
   printf("Hello from C!\n");
   fortransub();
   for (i=0; i<5; i=i+1)
      printf("%d ", numbers[i]);
   printf("\n");

return 0;
}

Both the subroutine fortransub() and the integer array numbers are declared as being external to the program. The function fortransub() can then be called in the program, and the values it modifies in the array numbers are printed out in the C program.

The programs can then be compiled in the following manner:

gcc -c ccode.c
gfortran ccode.o fortrancode.f95

When the program (a.out) is executed, here is the output:

Hello from C!
 Hello from Fortran!
1 2 3 4 5

Here is a great resource for C-Fortran interoperability.

 

 

Partial re-engineering and frameworks (entombing)

Not all code re-engineering involves taking an old piece of code and re-writing the whole program. As Fortran is backwards compatible, it is often not too difficult to get the old program running, assuming it isn’t already. An easier form of re-engineering involves parts of programs, typically subprograms, which exist, but have no framework, i.e. a program, supporting them.

Take for example, this old Fortran subroutine designed to calculate “Site index curves for Douglas-fir in New Mexico“, from 1976.

Given just this subroutine, there is one of two options, both of which require the development of a framework around the code. The first option involves merely entombing the old code within a modern framework, say F95. This would first require the construction of the framing program, and then entombing the subroutine DFSITE within it. This would involve replacing features which have been deprecated, because the code is now a hybird, i.e. not completely Fortran 95. In most cases this is elements like comment delimiters, and in the case of the code above the use of the label 1 as a line continuation symbol – replaced with & on the end of the previous line. You also can’t use “implicit none”, unless you want to explicitly declare all the variables, but then you are re-engineering the subroutine. Here is a sample of the subroutine entombed in a program:

The second option is performing a re-engineering of the subroutine before placing it in the framework. A re-engineered piece of code might look something like this:

subroutine dfsite(ht,age,site)
   intrinsic none
   real, intent(in) :: ht, age
   real, intent(out) :: site
 
   integer :: i1, i2
   real :: a1, a2, htsia
   real, dimension(25) :: htsi
   htsi = (/4.5, 15.0, 25.5, 31.8, 36.4, 40.3, 43.7, 46.7, &
            49.5, 52.0, 54.3, 56.4, 58.4, 60.2, 61.8, 63.2, &
            64.4, 65.4, 66.2, 66.8, 67.2, 67.5, 67.7, 67.8, 67.8/)

   tem = age + 0.001
   if (tem > 230.001) then
      tem = 230.001
   end if
   i1 = (int(tem)/10) + 1
   i2 = i1 + 1
   a1 = real(i1-1) * 10.0
   a2 = a1 + 10.0
   htsia = (htsi(i1)*(tem-a2)-htsi(i2)*(tem-a1))/(a1-a2)
   site = 54.2 * ht / htsia
   return
end

 

 

Arrays and memory in Fortran

As far as programming languages goes, Fortran is only marginally slower than C, and yet when it comes to arrays way more flexible. From the perspective of large arrays, Fortran does not force the programmer to have to rely on the use of dynamic arrays (which can be helpful for those that despise pointers). The behaviour of Fortran’s automatic arrays depends on the type of compiler. Using Intel Fortran means they are stored on the stack (by default), whereas GNU Fortran stores them on the heap. Using the heap to store arrays is extremely useful, as the programmer does not have to worry so much about creating an array that might cause a stack overflow (as it might if the array is stored on a stack).

Consider the following Fortran program:

program arraymem
   implicit none
   call meanArray(8000000)

contains

subroutine meanArray(n)
   integer, intent(in) :: n
   integer, dimension(n) :: a, b
   a = 12
   b = 255
   a = (a + b) / 2
   print *, sum(a) 
end subroutine meanArray

end program arraymem

Compiling it with gfortran will store all 8 million elements of the arrays a and b in the heap. You can also declare them as dynamic using the keyword allocatable, but it isn’t required (eat that C!) Now Fortran’s use in industrial applications  is heavily contingent on the use of arrays, so they have to work. Automatically using the heap for array storage adds a layer of transparency for the programmer, and for novice programmers that’s a good thing because it reduces their cognitive load and allows them to concentrate on the problem, rather then get bogged down in syntax.

 

 

Coding Fortran: Allocatable (dynamic) arrays

The last couple of posts covered what we would classically call static arrays, now let’s use the same task, i.e. calculating the median of a set of numbers, using allocatable, or dynamic arrays in Fortran. Here is the program (with minimal modifications really):

program median
   use sorting
   implicit none

   integer :: num
   integer, allocatable, dimension(:) :: theArray
   real :: m

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

contains

subroutine inputNumbers(A,num)
   implicit none
   integer, allocatable, dimension(:), intent(out) :: A
   integer, intent(out) :: num
   integer :: i

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

   allocate(A(num))

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

end subroutine inputNumbers

real function calcMedian(A,num)

   integer, allocatable, dimension(:), intent(in out) :: A
   integer, intent(in) :: num
   integer :: aD
   aD = 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, allocatable, dimension(:) :: theArray

It uses the keyword allocatable, and sets the dimensions to :, implying that they will be set later. Now when the array is passed to the subroutine inputNumbers(), no dimension is required. Note the header for the subroutine:

subroutine inputNumbers(A,num)  
   implicit none  
   integer, allocatable, dimension(:), intent(out) :: A  
   integer, intent(out) :: num

There are two parameters. The first parameter, A, specifies the array being passed, while the second parameter num specifies the size of the array being passed back from the subroutine.

Inside the body of inputNumbers() is a single statement, occurring after the user has specified the number if integers to input:

allocate(A(num))

This allocates the array in dynamic memory, creating enough storage for num integers. Other than that, it is treated in the same way any array is treated.

Character arrays vs strings in Fortran

In many languages the concept of a string and a character array are similar. Not so in Fortran. Here they are distinct entities.

A character array in Fortran is defined as follows:

character, dimension(12) :: chrA

A character string on the other hand has a length, the number of characters the string contains. For example:

character (len=12) :: str

or in legacy Fortran:

character str*12

If the length is omitted, it is assumed to be 1. Strings and character arrays are not the same.

Here is how they differ:

A string declaration defines the maximal length of the string. Assignment of a shorter string results in the automatic padding of the rest of the string with blanks. For example:

str = ‘FORTRAN’

|--------------- Physical Length ---------------|
+---+---+---+---+---+---+---+---+---+---+---+---+
| F | O | R | T | R | A | N |   |   |   |   |   |
+---+---+---+---+---+---+---+---+---+---+---+---+
|------ Logical Length -----|---- Blank tail ---|

So using the declarations above, str will read in a string a string up to 12 characters in length, whereas chrA will read in 12 characters individually until the array is full:

read(*,*) str

read(*,*) chrA

 

 

 

Coding Fortran: modules

Modules are somewhat similar to libraries in C (i.e. the .h, .c file combination). They are easy to create because they only require one file. Modules are useful for:

  • Defining global variables.
  • Subroutine/function interface information is generated to aid in checking that proper arguments are passed.
  • New data types can created.

We will integrate the subroutine version of the factorial calculations described previously into a module. Here is what the module (factmodule.f95) looks like:

module factmodule

contains

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

end module factmodule

LEGEND: module wrapper, subroutine

It can be accessed from a calling program (factcallModule.f95) in the following manner:

program factorial
   use factmodule
   implicit none
   integer :: i, num

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

   call fact(num)

end program factorial

And is compiled in the following manner:

gfortran factcallModule.f95 factmodule.f95 -o main

Alternatively, it is possible to create an object file, which will result in both .mod, and .o files:

gfortran -c factmodule.f95

Note that you should use filenames with all lowercase for modules, otherwise factModule.f95 would result in the creation of factModule.o and factmodule.mod, possibly resulting in problems trying to find a file.

This can then be compiled with the main program in the following manner:

gfortran factcallModule.f95 factmodule.o

Alternatively, you can just compile them like this:

gfortran factcallModule.f95 factmodule.f95

 

 

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, temp
   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.

 

 

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).