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

Help
1

Hello everyone,

To reduce repeated code inside a subroutine, as in this example where we set the node coordinates of a cell, I was thinking if we could create a local anonymous subroutine that would help us do this repeated operation in a cleaner way.

Surprisingly we can.

Using contains statement, I was able to create a nested subroutine test and use it for setting the value of cc to a dummy value of -999.

I verified that this works on Intel, GNU, and AOCC Fortran 2008 compilers.

But since I never saw anyone else use this kind of nested subroutines/functions in OpenSource Fortran codes, so I have to ask if this is good practice or not.

Thanks

code

Edit: A more cleaner example of the same code …

code

2 Likes
2

I have used this in some of my codes. I do not know any reason why it should be considered bad practice, but maybe some more opinions should be listened at this point.

2 Likes
3

It is standard fortran since f90, so it should work for any modern (post-f77) compiler.

There are perhaps a couple of comments. You use implicit none within the subroutine, perhaps thinking that this will help identify if you were to mistype a variable or to inadvertently reuse the name of a temporary variable from the host routine. It doesn’t. The contained subprogram still has access to all of the host subroutine variables, or to the module variables if this occurs within a module. The burden to avoid making these kinds of mistakes is more on the programmer than the compiler in these cases, so beware.

The contained subprogram is not entirely anonymous. The host subprogram can still use the contained routine as an actual argument so that other subprograms can reference it (as a function, or with a call statement) through the dummy argument. This is in fact a nice way to expose what would otherwise be local data, while limiting the global namespace pollution. I think procedure pointers can work the same way, but with the restriction that the host subprogram is in the active call tree.

You are only allowed one level of “contains” in an external program, and two levels of “contains” within a module. There are situations where it might be nice to organize an algorithm around multiple nested contains levels, but that is not allowed.

Finally, one use for a contained subprogram is to localize the variables associated with some operation. There are now other ways to achieve that goal, namely the BLOCK and ASSOCIATE constructs. So when the contained subprogram is called from the host routine in multiple places, the contained subprogram is probably the right approach. But if it is called from just a single place in the host routine, then sometimes the BLOCK or ASSOCIATE constructs are a better fit.

6 Likes
4

There are two aspects to this: code readability (which also often means maintainability) and performance:

  1. for you and your team (who may be you yourself in various future “avatars” as you reapproach the code at different stages in your development), do such internal subprograms make it easier to understand (and maintain) your code; or, is it better to have them as MODULE subprograms?

  2. check with the processors you intend to use and how they optimize internal subprograms or not and whether overhead costs can be incurred in your “hot” loops.

Separately, look into using a Fortran 2018 compiler, or one that supports the Fortran 2018 enhancements with the IMPORT statement: c.f. WG5 document N161 - you may have meant IMPORT instead of “implicit none”!

image
image729Ă—894 76.8 KB

4 Likes
5

It is always good to create a heirachy of subroutines/functions, so that repeated actions or groups of actions can be easily defined, without duplicated code.

One aspect of using contains, that needs to be understood, is the scope of variables that are being used in the containing routine. When a sub/fun is contained, all variables that are in scope in the containing routine are available to the contained routine, unless there is a conflict in the named arguments of the contained routine.

This means that, unlike an uncontained routine, there are many other variables available to the routine that are not explicitly defined in this routine. (consider local DO indexes)
This can be contrary to the strategy of using a regular routine where the scope/impact of this simpler routine can be more clearly identified.
In some cases, it might be preferable that this simpler routine is not contained.

Nevertheless, the idea of nesting actions in a set of routines is an important approach for software development when using Fortran.

4 Likes
6

I don’t love using contains’d subprograms specifically because they don’t actually get their own scope. They have access to all the same variables available in whatever routine contains them. I find that can lead to unintended bugs, or just be annoying as I prefer to use i,j,k for loop indices everywhere, which isn’t great if the next function down a call stack now has access to the same i,j,k…

Instead, I just make a new routine in the same module, call it wherever I wanted, and don’t make it public. This way, it truly gets its own scope, and it’s still not callable by anything outside the module.

4 Likes
7

You have a point here. I use contains’d subprograms as a replacement of the deprecated inline functions. For anything any complex than that I use an independent subprogram in the module.

Nevertheless, your argument on the scope of variables is absolutely true.

1 Like
8 
   import [[ :: ] import-name-list ]
or import, only : import-name-list
or import, none
or import, all

introduced starting Fortran 2018 standard revision provides better control about variable scope in the internal subprogram.

3 Likes
9

That’s interesting I’ll have to give that a try. Do you know if it’s implemented in gfortran, ifort, or ifx?

I will say it’s still kind of a miss along the same lines as requiring implicit none, requiring interface blocks import everything again and declare implicit none again, etc. But, at least it’s possible to have an import none so that’s great.

10

The code

program main
implicit none
integer :: i, j, k
i = 2
j = 3
k = 4
call display(i)
contains
subroutine display(i)
import, only: j
integer, intent(in) :: i
print*,i,j
! line below does not compile, because k is not imported
! print*,i,j,k
end subroutine display
end program main

is compiled by ifort Version 2021.6.0 and ifx, giving output 2 3, but not by gfortran 14.0.0 20230604.

11

@Beliavsky ,

What is “does not compile” ? Possibly undeclared variable ?
What is the scope of implicit none in this example, given import, only ?

12

Activating the print*,i,j,k line, ifort says

xsub.f90(14): error #6404: This name does not have a type, and must have an explicit type.   [K]
print*,i,j,k
-----------^
compilation aborted for xsub.f90 (code 1)

The implicit none in the main program also applies to the contained subroutine. Removing the implicit none, so that the code is

program main
integer :: i, j, k
i = 2
j = 3
k = 4
call display(i)
contains
subroutine display(i)
import, only: j
integer, intent(in) :: i
print*,i,j
print*,i,j,k
end subroutine display
end program main

ifort compiles it and gives output

           2           3
           2           3           0

although of course k having value 0 in the display subroutine is not something one should count on.

13

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
14

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

15

@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
16

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

1 Like
17

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
18

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
19

@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

20

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.