Is creating nested subroutines/functions considered good practice in Fortran?

Nested procedures offer an excellent way of creating wrappers that can be subsequently passed to other procedures as procedure arguments, much better than object-oriented style solutions, in my opinion. It’s a powerful feature of Fortran 2008. It was a source of trouble and segmentation faults with Microsoft WSL1 a long time ago and later Apple M1. Both are now resolved.

2 Likes

Could you please show an example of this? Or if you know any website that explains it, that would also be useful.

@aerosayan , you may want to keep “Modern Fortran Explained” handy as that can provide a ready handbook reference with a fair amount of prose to provide the background on the semantics and syntax of the ISO IEC language standard for Fortran which online forum comments cannot quite provide.

With that in mind, you can review this example:

module solver_m
   abstract interface
      function Ifunc( x ) result(r)
         ! Argument list
         real, intent(in) :: x
         ! Function result
         real :: r
      end function 
   end interface
contains
   subroutine eval( x, func )
      ! Argument list
      real, intent(in) :: x
      procedure(Ifunc) :: func
      print *, "In eval: f(x) = ", func(x)
   end subroutine 
end module 
module calc_m
   use solver_m, only : eval
contains
   subroutine sub( x )
      ! Argument list
      real, intent(in) :: x
      call eval( x, myfunc )   !<-- NOTE: `myfunc` as procedure argument 
   contains
      function myfunc( x ) result(r)
         ! Argument list
         real, intent(in) :: x
         ! Function result
         real :: r
         real, parameter :: a = 1.0, b = 2.0, c = 1.0
         r = a*x**2 + b*x + c
      end function 
   end subroutine 
end module
   use calc_m, only : sub
   call sub( x=99.0 )
end 
  • The response you can expect by a processor as confirmed by gfortran on Windows is as follows:
C:\temp>gfortran -ffree-form p.f -o p.exe

C:\temp>p.exe
 In eval: f(x) =    10000.0000

Now, note my strong recommendation to Fortranners will be to strive for IMPORT, NONE if they decide to pursue such internal subprograms. However, the gfortran users will have to become developers or support development to enhance the GCC/gfortran to support this Fortran 2018 feature.

..
   contains
      function myfunc( x ) result(r)
         import, none
         ! Argument list
         real, intent(in) :: x
..
  • response using Intel Fortran:
C:\temp>ifort /free /standard-semantics p.f
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.10.0 Build 20230609_000000
Copyright (C) 1985-2023 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.36.32537.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
 In eval: f(x) =  10000.00
C:\temp>ifx /free /standard-semantics p.f
Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2023.2.0 Build 20230627
Copyright (C) 1985-2023 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.36.32537.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
 In eval: f(x) =  10000.00
2 Likes

See this comment below: not quite implemented yet with gfortran.

1 Like

The concept is known as a callback. In other words you have a subroutine that takes another subroutine or function as argument.

The classic example would be a numerical solver and a parametrized function of some sort. We can take root-solving as an example:

  ! find the root f(x) = 0, in the interval [ax,bx]
  x = zeroin(ax,bx,f,tol)

The function f passed as a parameter to zeroin must match a particular interface:

interface
  real(dp) function f(x)
    import dp
    real(dp), intent(in) :: x
  end function
end interface

To make this more concrete, say you had to find the friction factor f prescribed by the Colebrook-White equation:

\frac{1}{\sqrt{f}} = -2 \log \left( \frac{\epsilon}{3.7 D_\text{h}} + \frac{2.51}{\mathit{Re} \sqrt f} \right)

Finding f can be recast as finding the root of the equation:

G(f) = \frac{1}{\sqrt{f}} +2 \log \left( \frac{\epsilon}{3.7 D_\text{h}} + \frac{2.51}{\mathit{Re} \sqrt f} \right) = 0

In Fortran you might do this as follows:

real(dp) :: f, Re, Dh, eps
real(dp), parameter :: tol = 1.0e-9_dp

Re = 10000
Dh = 0.2_dp    ! in meters
eps = 0.002_dp ! 2 mm

f = zeroin(0.001_dp, 0.01_dp,friction_equation,tol)
print *, f

contains

  real(dp) function friction_equation(f) result(y)
    import, only: Re, Dh, eps  ! if supported by your Fortran compiler
    real(dp), intent(in) :: f
    y = 1.0_dp/sqrt(f) + 2*log(eps/(3.7_dp*Dh) + 2.51_dp/(Re*sqrt(f))) 
  end function

This could be made part of general subroutine taking Re, D_h and \epsilon as inputs, making sure they are valid and finally returning the friction factor f as output.

2 Likes

For reference (not advocating it now that IMPORT is defined in the
standard) note that a BLOCK/ENDBLOCK in the containing procedure as
well as a nearly null containing procedure have been used to simulate
the behavior IMPORT now provides.

I remember that the default behavior where everything declared in
the containing procedure is in scope in the contained procedure was
surprising to me when I first tried contained procedures. I remember
looking for something to turn it off, and thought “IMPLICIT NONE”
explicitly in the contained procedure might turn off the “importing”;
so I missed something like IMPORT immediately.

The two work-arounds for providing IMPORT-like functionality can be
confusing if the intent is not made clear with a comment, which is
a red flag there is missing functionality in the language, now resolved
via IMPORT.

If a compiler that supported IMPORT had a flag to warn if a contained
procedure did not have an explicit IMPORT statement I would use it.
Once IMPORT is widely available I expect my practice will be to always
have an IMPORT statement in a contained procedure for clarity (like IMPLICIT NONE, some extra typing will be required to get what in my opinion is the default expected behavior in new code).

Artificial examples of the two ways to provide what IMPORT
now provides, just for demonstration purposes:

! nothing to inherit because everything is in contained procedures
program main
!implicit none
call top()
contains

subroutine top()
  integer,parameter :: bad=-huge(0)
  integer :: i=bad,j=bad,k=bad
  call printme()
end subroutine top

subroutine printme() ! intentionally undefined i,j,k 
   write(*,*)i,j,k
end subroutine printme

end program main
! exclude what you do not want in scope to contained
! procedures by declaring the excluded variables in a BLOCK
program main
! implicit none
localize: block
  integer,parameter :: bad=-huge(0)
  integer :: i=bad,j=bad,k=bad
  call printme()
endblock localize
contains

subroutine printme()
   ! i,j,k not defined
   write(*,*)i,j,k
end subroutine printme

end program main

urbanjs@mercury:~/github/fun/try_import/app$

2 Likes

@urbanjost
I find your block example confusing.
printme is contained by program main (not block), so as i,j,k are local to block, are they available to printme ?
I still use a simpler structure where scope implications are easier to understand

Proving a negative (that i,j,k are not defined) would have been better illustrated with an allocatable variable perhaps, as doing “normal” things like an implicit none make the code
not compile and so on. It is not totally intuitive but the scope of variables declared in the block
are not seen as host variables of the broken but are only defined in the confines of the block/endblock so they are not seen by the contained procedures. In this case the compiler errors might make it clear; so if you add the IMPLICIT NONE to the top of the program you will see the compiler complains that I,J,K in PRINTME() are not defined.

When a compiler does not support IMPORT I sometimes just put a block and endblock in the procedure hosting the contained procedures to make sure I get only complaints about the variables I wanted and then remove them, with commented out IMPORT statements to remind me how I actually want it done :> The BLOCK/ENDBLOCK syntax is my favorite for not exporting anything I do not want exported, as anything above the BLOCK statement is visible to the contained procedures and the statements can have a label like LOCALIZE that gives some idea as to their purpose; but even BLOCK is not implemented in all processors so the near-empty hosting procedure syntax is currently more portable.

Over time I’ve come to believe that callbacks, procedure arguments and procedure pointers are under-utilized by Fortran practitioners. The majority of Fortran codes remain written in a procedural style.

If we think about which codes have lasted long – ODEPACK, MINPACK, COBYLA, SLSQP, ARPACK … – and are reused by hundreds of programmers including those from other programming languages, what many of these program have in common is they accept a procedure parameter or provide a reverse communication interface. Both of these are a form of dependency injection. The reasons these programs are long-lasting is because they help achieve separation of concerns. A numerical expert writes an algorithm; the users just need to provide their problem and receive immediate value from such programs.

Rob Pike, a famous C programmer, has captured this essence of function pointers (procedure arguments in Fortran) in his Notes on Programming in C (pg. 5),

Another result of the tyranny of Pascal is that beginners don’t use function pointers. (You can’t have
function-valued variables in Pascal.) Using function pointers to encode complexity has some interesting
properties.

Some of the complexity is passed to the routine pointed to. The routine must obey some standard
protocol — it’s one of a set of routines invoked identically — but beyond that, what it does is its business
alone. The complexity is distributed.

There is this idea of a protocol, in that all functions used similarly must behave similarly. This
makes for easy documentation, testing, growth and even making the program run distributed over a network — the protocol can be encoded as remote procedure calls.

I argue that clear use of function pointers is the heart of object-oriented programming. Given a set of
operations you want to perform on data, and a set of data types you want to respond to those operations, the
easiest way to put the program together is with a group of function pointers for each type. This, in a nutshell, defines class and method. The O-O languages give you more of course — prettier syntax, derived
types and so on — but conceptually they provide little extra.

Combining data-driven programs with function pointers leads to an astonishingly expressive way of
working, a way that, in my experience, has often led to pleasant surprises. Even without a special O-O
language, you can get 90% of the benefit for no extra work and be more in control of the result. I cannot
recommend an implementation style more highly. All the programs I have organized this way have survived comfortably after much development — far better than with less disciplined approaches. Maybe
that’s it: the discipline it forces pays off handsomely in the long run.

When I think of scientific computing, here are some of the callbacks which immediately come to my mind:

  • scalar functions (root solving, integration, fitting)
  • multivariate functions (systems of nonlinear equations, systems of ODEs, integration)
  • matrices (typically a Jacobian of some sort, or a sparsity pattern)
  • operators (a transformation returning a vector, e.g. a matrix-vector multiply)

Through composition, you can arrive at quite complex cases. For example you may need to solve a boundary value problem (BVP). For this you may use the shooting method where you combine an ODE solver and a root-finder. Afterward you’d like to fit the BVP to experimental data using a least squares procedure. In this case you will have a hierarchy of callbacks working in sync.

A quite advanced example of callback usage is the Unified Framework for Finite Element Assembly (also known as UFC), which is an interface between the problem-specific and general purpose components of a finite element program used in the FEniCS project. Although written in C and C++, it allows one to provide a set of callback functions to assemble the bilinear and linear forms originating from the weak form of the PDE. The general parts including assembly and solving the system of equations is handled in the C++ framework.

In other fields of computing callbacks are used for all sorts of things:

  • comparison functions for sorting (see qsort in C)
  • drawing and GUI callback functions
  • audio waveforms (see Open Sound Control - #2 by ivanpribec)
  • teardown functions (see atexit in C)
  • games (switching between different behaviors, e.g. spells or weapons)
3 Likes

@FortranFan your example works with LFortran also:

$ lfortran x.f90 
In eval: f(x) =  1.00000000e+04

I must say I am pleasantly surprised that it actually worked! We are in alpha, meaning I expect things to break. :slight_smile:

@ivanpribec yes, I use nested functions for callbacks often, I think it’s a clean and fast way to do it.

4 Likes

I’m developing an Equations of State library to calculate thermodynamic properties. One of the key principles is that everything can be resumed as a single function of residual Helmholtz energy

A_r(z, V, T)

and with that function and its derivatives, one can calculate most of the thermodynamic and phase equilibria properties. One of the main complexities of this is that there are many different expressions for this function (a common joke is that each research group has its own expression for Helmholtz energy). This is the only thing that is required to define a model, I’ve seen a lot of codes (mostly the F77 codes I’ve inherited from my supervisor) that use some flag and a selector to use the correct function based on the desired model. This is very error-prone since every time you want to add a model you have to edit the source code and be careful to find all the occurrences of this (also it is pretty dirty having to edit the main source just for this)

After learning about procedure pointers it clicked right away, it’s a pretty simple way to add extensibility without going straight into derived types. I’ve started developing (still a WIP) a new framework to ease the process and I’m very happy with the results so far. No more need to edit old sources and feel dirty, now every new model can be a single module that contains its parameters, its hidden procedures for any inner workings, and a setup procedure that sets the needed procedure pointer.

I love easy extensibility and couldn’t find any other library (in any language) that could have a process so simple as this one, even better with how simple it is to write vectorial math in Fortran. A little example of how this is done:

here is where the procedure is held and it’s setter routine:

I’ve ended up writing a lot more than I thought but I wanted to express how procedure pointers can be really helpful with so little effort to implement!

2 Likes

I like pure functions, modern replacement for statement functions to do this if they fit your needs . They should be in lined and have minimal impact on performance.

2 Likes

Hi all,

I just want to chime in on the performance point. Some implementations use trampolines to achieve calling Fortran internal subprograms that use variables of the host subprogram via host association. This extra “indirection” does have performance impact, but it is probably less than the consequences of poor optimizations resulting from the call via a procedure pointer itself (comparing to a direct call that might be inlined, enabling more optimizations like vectorization of the loop containing the call). So I think it is not the best approach to use procedure pointers in performance critical spots.

3 Likes

@szakharin Yes. The trampolines are also problematic that they require executable stack (at least in the implementations I’ve seen), which is often not allowed. We implemented the nested functions (closures) by creating a module and sharing the needed variables via pointers. I think this approach is cleaner and should be faster, but we have not done any thorough benchmarking yet.

2 Likes

Quite interesting! As far as I understand, writing import, all is redundant, in the sense that the deafult bahavior is that the contained function/subroutine has access to all variables defined in the containing program. Is that correct?

By default, the internal subprogram indeed “imports” all the objects from the containing scope.

What IMPORT, ALL provides is documentation of that situation - the readers of your program (who may be you yourself in a future avatar such as some period of time later where the state of memory is different!) may find that clearer to understand / beneficial .

1 Like

This is interesting. Is it working with recursion?

Example:

module other
  integer, parameter :: CONST_N = 10
contains
  subroutine foo(fptr)
    integer :: fptr
    print *, fptr()
  end subroutine foo
end module other

recursive subroutine host(local)
  use other
  integer, value :: local
  call foo(callee)
  return
contains
  function callee()
    integer :: callee
    callee = local
    if (local .le. CONST_N) then
       call host(local + 1)
    endif
  end function callee
end subroutine host

As others already pointed out, the contained procedure is fully aware of everything defined in the main program. This abuses the data hiding principle (often called encapsulation or information hiding.) I would strongly recommend to avoid abusing data hiding unless it’s just a trivial case and you know what you are doing. For that reason I only use contains in trivial cases, where I just want to quickly add an auxiliary procedure (which typically isn’t long and doesn’t need many variables anyway.)

Other than trivial cases, I think it’s best to avoid contained functions, unless your compiler supports import, none and import, only, which is the solution to the data hiding abuse. But since I mainly use gfortran and that feature is not yet supported, the best strategy is to either implement a module or write the auxiliary procedure in the same file, but after the end program statement, like this:

program main
implicit none
...
call foo(...)
...
end program
!-----------------------------------------------------------------------------
subroutine foo(...)
! This subroutine is completely independent and unaware of whatever is defined
! in the calling program.
...
end subroutine foo

Yes this is a bit old fashion, but it works: the auxiliary function does respect the data hiding principle, and it’s still right there, easy to see and maintain. The drawback is when the main program uses several other modules and you want the auxiliary function to use them as well. If you write an “independent” procedure in the same file as above, you would have to rewrite the use statements in that procedure.

Writing a module, declare it private and have the auxiliary procedure(s) there is, without any doubt, the way to go for larger programs. In simple cases, however, when you just want a quick auxiliary procedure, implementing a module just for that seems overkill. Although import, only and import, none will solve the issue when it is widely available in compilers, I would argue that this also has a potential drawback, especially for beginners. One could have loads of auxiliary procedures then, contained in the main program, and still independent (or partially independent, as appropriate.) This shouts “create a module instead”, but it will be possible without a module, which I consider bad programming style.

Not in our current implementation, but in general yes: you have to maintain a stack in the module, and each recursion increases the stack position. I think the NAG compiler does exactly that. We’ll implement that later.

And to answer the initial question, it is neither a good nor a bad practice in itself, all depends on the context.

1 Like