More on functions and subroutines

Yutaka Masuda

February 2020

Back to index.html.

More on functions and subroutines

Basic usage

Fortran: a procedural language

Fortran is categorized as procedural language, which is one of the programming paradigms. A procedural language has procedures, and each procedure has a group of instructions that is structured with flow-control statements (if and do). The programmer describes a program to call a sequence of procedures. Some procedures and related variables are packed into a module. Fortran shows its capabilities through procedural programming.

A procedure is a program unit that performs a particular task. Fortran has two types of procedures: function and subroutine. A function receives the arbitrary number of arguments (input values) and returns one value. A subroutine receives the arbitrary number of arguments (input values) and returns the arbitrary number of output, possibly zero output.

There are many reasons why it is better to use procedures. One of the advantages is better readability that makes the programmers quickly understood the code. For example, the following is an updated program that I have shown in the introduction to compare to R.

program compute_mean_and_sd
   use data_manipulation
   implicit none
   real :: val,mean
   real,allocatable :: data(:)

   n = get_number_of_records("data.txt")
   allocate(data(n))
   call read_reacords("data.txt",n,data)

   print *,"n    = ",n
   print *,"mean = ",get_mean(data)
   print *,"sd   = ",get_sd(data)
end program compute_mean_and_sd

We assume that all the procedures are in a custom module, data_manipulation. This program should be more transparent for anybody than the original program.

Another advantage is reusability. Once you write a procedure, you do not have to write the same instructions anymore. You can give the (external) procedure to another programmer who can use it in various programs. You may have many procedures, and a module becomes a collection of procedures. The module efficiently exports its procedures to other programmers with less effort.

The other one is the ease of maintenance and tests. If a procedure is aimed to have a single task, a modification of the procedure is straightforward. A short procedure is also helpful to test. Although you have to write a sophisticated program, each component should be relatively simple, and the element can form a procedure or a module.

Definition of function and subroutine

Function

A function is defined by the following form. You should write implicit none if it is an external function i.e., defined outside a program unit.

function name(arguments...) result(variable)
   {implicit none}
   argument_type,attribute :: argument_list
   result_type :: variable

   {main body}
end function name

A simple example is as follows.

function sumup(n) result(val)
   integer,intent(in) :: n
   integer :: val,i
   val = 0.0
   do i=1,n
      val = val + i
   end do
end function sumup

A function receives the input as variables called arguments (n above). The argument variable (known as dummy variable or dummy argument) is a reference to the original data. An argument has a type with attributions, other characteristics to clarify the role of the argument. A typical attribute is intent() which defines whether the argument is supposed to be changed or not.

It is the programmer’s responsibility to assign a value to the return variable. Without any assignment to the variable, the return value is undefined, and it causes trouble.

Subroutine

A subroutine is defined by the following form. As I mentioned, you should put implicit none if it is external.

subroutine name(arguments...)
   {implicit none}
   argument_type,attribute :: argument_list

   {main body}
end subroutine name

A simple example is as follows.

subroutine sumup(n,val)
   integer,intent(in) :: n
   integer,intent(out) :: val
   integer :: i
   val = 0.0
   do i=1,n
      val = val + i
   end do
end subroutine sumup

The argument can have the same attributes as used in a function.

Type of procedures and scope of variables

External procedure

The procedure is external if it is defined as an independent program unit (i.e., no “parent” procedure).

Internal procedures

The procedure is internal if it is nested within another program unit under contains.

Arguments

Optional arguments

It is useful to have optional arguments in function or subroutines. For example, the following subroutine shows a vector, optionally with a text message.

subroutine print_vector(v,message)
   real,intent(in) :: v(:)
   character(len=*),intent(in),optional :: message

   if(present(message)) then
      print "(a)",trim(message)
   end if
   print *,v
end subroutine print_vector

You can call this subroutine in two ways.

   real :: v(3)
   ...

   ! usage 1: only values
   call print_vector(v)

   ! usage 2: a message with values
   call print_vector(v,"values")

   ! another form
   call print_vector(v,message="values")

Using this feature, you can unify several subprograms into a single form. You can define more than 1 optional arguments. Here are the basic rules to use optional arguments.

When you access the optional argument, but the argument is not given, you hit a weird error. Some compilation options can detect this error.

Allocatable arguments

An allocatable array can be allocated within a program unit that declares the array. Also, you can allocate an array in a function or a subroutine. In this case, the array argument should have the allocatable attribute.

Here is a code to allocate and initialize an array within a subroutine. The array is a temporary vector, and its size is determined with some formula in the subroutine. A caution is that you should care about a case where the array is already allocated.

subroutine allocate_temporary_vector(v,m,n)
   real,intent(inout),allocatable :: v(:)
   integer,intent(in) :: m,n

   if(allocated(v)) then
      deallocate(v)
   end if
   if(max(m,n)<1) then
      allocate(v(2))
   else
      allocate(v(max(m,n)))
   end if
   v = 0.0
end subroutine allocate_temporary_vector

Longevity of local variables

Longevity of variable: basics

As seen in the previous chapter, a variable has a scope. A variable is valid within a program unit that has defined the variable. If the program init is done, the value of the variable disappears.

It is valid for allocated arrays. The allocated array survives only within a program unit where the array is defined. If the procedure that defines the allocated arrays is over, the program automatically deallocates the allocated arrays.

Do not be confused with the allocatable argument. It just receives an allocatable array defined outside the unit, so the caller defines the longevity.

Saved variables

You can define a variable that holds value after the program unit is over. With this feature, a procedure keeps the previous value and reuse it in the next call. The save option makes a local variable permanent.

function normal_pdf(x) result(f)
   real,intent(in) :: x
   real :: f
   logical,save :: first=.true.
   real,save :: w

   ! calculate w only for the first time
   if(first) then
      w = 1.0/sqrt(2.0*3.1415926)
      first = .false.
   end if

   ! pdf of standard normal distribution
   f = w*exp(-0.5*x*x)
end function normal_pdf

This function uses 2 saved variables; first indicates whether the call is the first time or not, and w holds a constant (\(1/\sqrt{2\pi}\)). For the first call, the function calculates w and first becomes .false. so that w is never recalculated in the subsequent calls. The declaration of saved variables can be simplified by assigning a value. See a piece of code.

   real,intent(in) :: x
   real :: f

   ! same as the "save" attribute
   logical :: first=.true.
   real :: w=0.0

That is, when you assign a value to a local variable, the variable must have the save attribute. It may cause unexpected results, and it is tough to detect such a problem. For this reason, you should not assign any values to local variables in the declaration block. If you need the saved variable, you should use the save attribute explicitly.

You may wonder if the save statement is needed to write a program. In the above case, you can pre-calculate w somewhere else and use it so that no save statement is needed. Or, you can pass w as an additional argument (it is ugly but may be useful for another case). In my opinion, you can use save only if you think it is useful because saved variables may be tricky, and it creates side effects in the program. A case where save is useful, is like a function that needs many constants and temporary variables, and it is awkward to pass them through the arguments.

Recursion

A procedure can call itself with different arguments. Mathematically, it is called a recursive formula. Fortran supports recursions.

A simple recursive formula is factorial, \(n!\). Given \(n\), the factorial \(f(n)=n!\) is defined by recursion.

\[ \begin{aligned} f(1)&=1\quad\text{for}\quad n=1\\ f(n)&=n\times f(n-1)\quad\text{for}\quad n\geq 2\\ \end{aligned} \]

In Fortran, this function is written with the recursive keyword.

! It works when n>=1.
recursive function factorial(n) result(f)
   integer,intent(in) :: n
   integer :: f
   if(n==1) then
      f = 1
   else
      f = n * factorial(n-1)
   end if
end function factorial

The recursive function needs a condition that stops the recursion (\(n=1\) in the above code). It means, for many recursive formulas, you may use loops instead of recursions to do the same computations. You can see that the above code is easily rewritten using a loop.

function factorial(n) result(f)
   integer,intent(in) :: n
   integer :: f,i
   f = 1
   do i=1,n
      f = f * i
   end do
end function factorial

So, when should we use the recursive call? In my opinion, you do not have to use it. The recursion is useful in pedigree manipulation in quantitative genetics, e.g., the tabular method. But it may not be evident for the first view, although the recursive formula can be more straightforward than the loop in coding. Also, the function is called too many times, and it may fail.

Back to index.html.