I'm working on a modern Fortran version of the simulated annealing optimization algorithm, using this code as a starting point. There is a Fortran 90 version here, but this seems to be mostly just a straightforward conversion of the original code to free-form source formatting. I have in mind a more thorough modernization and the addition of some new features that I need.
Refactoring old Fortran code is always fun. Consider the following snippet:
C Check termination criteria.
QUIT = .FALSE.
FSTAR(1) = F
IF ((FOPT - FSTAR(1)) .LE. EPS) QUIT = .TRUE.
DO 410, I = 1, NEPS
IF (ABS(F - FSTAR(I)) .GT. EPS) QUIT = .FALSE.
410 CONTINUE
This can be greatly simplified to:
! check termination criteria.
fstar(1) = f
quit = ((fopt-f)<=eps) .and. (all(abs(f-fstar)<=eps)
Note that we are using some of the vector features of modern Fortran to remove the loop. Consider also this function in the original code:
FUNCTION EXPREP(RDUM)
C This function replaces exp to avoid under- and overflows and is
C designed for IBM 370 type machines. It may be necessary to modify
C it for other machines. Note that the maximum and minimum values of
C EXPREP are such that they has no effect on the algorithm.
DOUBLE PRECISION RDUM, EXPREP
IF (RDUM .GT. 174.) THEN
EXPREP = 3.69D+75
ELSE IF (RDUM .LT. -180.) THEN
EXPREP = 0.0
ELSE
EXPREP = EXP(RDUM)
END IF
RETURN
END
It is somewhat disturbing that the comments mention IBM 370 machines. This routine is unchanged in the f90 version. However, these numeric values are no longer correct for modern hardware. Using modern Fortran, we can write a totally portable version of this routine like so:
pure function exprep(x) result(f)
use, intrinsic :: ieee_exceptions
implicit none
real(dp), intent(in) :: x
real(dp) :: f
logical,dimension(2) :: flags
type(ieee_flag_type),parameter,dimension(2) :: out_of_range = [ieee_overflow,ieee_underflow]
call ieee_set_halting_mode(out_of_range,.false.)
f = exp(x)
call ieee_get_flag(out_of_range,flags)
if (any(flags)) then
call ieee_set_flag(out_of_range,.false.)
if (flags(1)) then
f = huge(1.0_dp)
else
f = 0.0_dp
end if
end if
end function exprep
This new version is entirely portable and works for any real kind. It contains no magic numbers and uses the intrinsic IEEE exceptions module of modern Fortran to protect against overflow and underflow, and the huge() intrinsic to return the largest real(dp) number.
So, stay tuned...
See also
- S. Kirkpatrick, C. D. Gelatt Jr., M. P. Vecchi, "Optimization by Simulated Annealing", Science 13 May 1983, Vol. 220, Issue 4598, pp. 671-680
- W. L. Goffe, SIMANN: A Global Optimization Algorithm using Simulated Annealing, Studies in Nonlinear Dynamics & Econometrics, De Gruyter, vol. 1(3), pages 1-9, October 1996.
- Original simann.f source code [Netlib]
It amuses me somewhat to see the push to get people to stop using Python 2. Python 2 was sort of replaced with Python 3 in 2008. However, Python 3 broke backward compatibility, and Python 2 continued to be supported and even developed (some of the Python 3 features have even been backported to Python 2 which is somewhat bizarre). However, the official Python maintainers have declared January 1st, 2020 to be the end of life for Python 2. However, many Linux distos (as well as MacOS) still include Python 2 as the default.
The situation with Fortran is significantly different. You can still get (even purchase) a Fortran compiler that can compile obsolescent FORTRAN 77 written 40 years ago. A lot of people haven't gotten the memo that Fortran has moved on. Not a day goes by where there isn't a StackOverflow question from somebody about some godawful FORTRAN 77 code that they can't figure out how to get working. So, I wonder how successful this effort will be with Python. Will there still be Python 2 holdouts 30 years from now? Of course, the Python community is significantly different from the Fortran community (to whatever extent Fortran could be said to even have a community). The Python implementation is open source so people are free in theory to just fork it and continue to update it. But, as the official source dries up I suspect eventually people will just move on.
There is basically no one centralized place for Fortran users to learn about Fortran, download Fortran, work on or even comment on changes to Fortran, or anything else really. Backward compatibility is actually one of the major strengths of Fortran, but there just isn't anyone to tell you to move on. Many of the major Fortran libraries freely available on the internet are still, in 2018, written in Fortran 77 (see NetLib, where Fortran codes go to die). There is also SOFA and SPICELIB (two libraries still being developed in Fortran 77 for some reason). There is LAPACK, probably the most visible FORTRAN 77 library. It turns out, you can't even link LAPACK and SPICELIB in the same program anymore, because both now have a routine called DPSTRF! If these libraries weren't using this depreciated source form and were actually using Fortran 90 modules (or god forbid, Fortran 2008 submodules) we wouldn't have this problem. Nobody should be publishing a Fortran library at this point that is just a pile of subroutines. And yet, they do.
Sure, there is a Fortran standards committee. But they write the standard and don't do much else. There isn't any official website that serves as any kind of central location to learn about the language or download any compilers. There is no nonprofit "Fortran Foundation". There isn't even a Fortran logo. What we basically have is a few compiler vendors such as NAG, Intel, PGI, as well as the open source gfortran developers, each with their own websites and their own schedules. Even the C++ language (another ISO standard language) has a GitHub site and even an actual homepage. The Fortran committee website is extremely underwhelming, and is basically just a list of links to plain text files that have meeting minutes and some other documents. The odds of a new user even finding this site are pretty slim (and of course, there really isn't anything useful here anyway). The Fortran standards process seems very opaque (I guess it involves infrequent meetings and typing up plain text files).
In many ways, the standards committee has failed the user base. The slowness of the language development process, the refusal to adopt modern practices of collaboration, and the refusal to address major shortcomings of the language has allowed other languages to overtake Fortran in the fields it was once dominant. In fields such as machine learning, one of the major computational activities of modern times, there is no significant presence of Fortran. The committee has given us incremental, somewhat half-baked features that don't really solve real-world problems in any kind of satisfactory manner (like ISO_VARYING_STRING, parametrized derived types, user-defined IO, floating point exceptions). While shortcomings from decades ago are left unaddressed (lack of a useful mechanism for generic programming, lack of access to the file system, no intrinsic high-level classes for dynamic structures such as stacks and queues, no useful varying string class, a non-standard module file format that is a constant source of annoyance for both users and library developers, even the lack of a standardized file extension for Fortran, which leads to nonsense like .f90, .f03, .f08, .f18? being used). It's incredible but true that the Fortran standard doesn't actually define a file extension to be used for Fortran. One of the respondents to a recent survey on Fortran language development had this to say:
The standard committee is too inbred with compiler developers who only see the language from the inside out and lacking in users who know what features they need for their particular application space.
Amazing libraries in the scientific/technical/numerical domain are being written in other programming languages. DifferentialEquations.jl is the type of high-quality math-oriented library that no one is writing anymore in Fortran (indeed it makes use of features of Julia that aren't even possible in Fortran). This article about using Julia as a differentiable language is the sort of thinking we desperately need in the Fortran world. Not more excuses about why fundamental changes can't be made to Fortran for fear of breaking somebody's 60 year old spaghetti code.
See also
- Ben James, Stop Using Python 2: What you Need to Know About Python 3, August 15, 2018 [hackaday.com]
- https://pythonclock.org
- PEP 373 -- Python 2.7 Release Schedule
- M. Innes, et. al., Building a Language and Compiler for Machine Learning, Dec 3, 2018 [julialang.org]
- Some ideas for Fortran, from a newbies perspective, July 22, 2019 [comp.lang.fortran]
I came across this old NASA document from 1963, where a cross product subroutine is defined in Fortran (based on some of the other routines given, it looks like they were using FORTRAN II):
A cross product function is always a necessary procedure in any orbital mechanics simulation (for example when computing the angular momentum vector \(\mathbf{h} = \mathbf{r} \times \mathbf{v}\)). The interesting thing is that this subroutine will still compile today with any Fortran compiler. Of course, there are a few features used here that are deemed obsolescent in modern Fortran (fixed-form source, implicit typing, and the DIMENSION statement). Also it is using single precision reals (which no one uses anymore in this field). A modernized version would look something like this:
subroutine cross(a,b,c)
implicit none
real(wp),dimension(3),intent(in) :: a,b
real(wp),dimension(3),intent(out) :: c
c(1) = a(2)*b(3)-a(3)*b(2)
c(2) = a(3)*b(1)-a(1)*b(3)
c(3) = a(1)*b(2)-a(2)*b(1)
end subroutine cross
But, basically, except for declaring the real kind, it does exactly the same thing as the one from 1963. The modern routine would presumably be put in a module (in which the WP kind parameter is accessible), for example, a vector utilities module. While a subroutine is perfectly acceptable for this case, my preference would be to use a function like so:
pure function cross(a,b) result(c)
implicit none
real(wp),dimension(3),intent(in) :: a,b
real(wp),dimension(3) :: c
c(1) = a(2)*b(3)-a(3)*b(2)
c(2) = a(3)*b(1)-a(1)*b(3)
c(3) = a(1)*b(2)-a(2)*b(1)
end function cross
This one is a function that return a \(\mathrm{3} \times \mathrm{1}\) vector, and is explicitly declared to be PURE (which can allow for more aggressive code optimizations from the compiler).
This document, which is the programmer's manual for an interplanetary error propagation program, contains a lot of other Fortran gems. It also has some awesome old school flow charts:
See also
There are various ways to interface Fortran and Python. However, many of the online examples on this subject seem to be content with a 40 year old variant of the language, and do not tend to focus on modern Fortran. One particular thing that is hard to figure out is how to send arbitrary or dynamic data from Python to Fortran (and vice versa)? This could be variable-sized arrays, variable-sized strings (or even a variable number of variable sized strings), or just any random data structure that contains different types of variables that aren't necessarily known at compile time. It is not easy to find out how to do this, but one option presented here is to use JSON as an intermediate format.
Sending a dict from Python to Fortran
Say we want to communicate an arbitrary Python dictionary structure to Fortran. We can do this fairly easily by dumping the dict to a JSON string, sending the string to Fortran, and then parsing the string in Fortran using JSON-Fortran. Then we have access to all the data in the dictionary in the Fortran code. An example of this is presented here.
First we need a function python_str_to_fortran() to convert a Python dict object into a c_char_p string that can be sent via the C interface to a Fortran routine:
from ctypes import *
import json
def python_str_to_fortran(s):
return cast(create_string_buffer(s.encode()),c_char_p)
def python_dict_to_fortran(d):
return python_str_to_fortran(json.dumps(d))
Note that we are using Python's ctypes module here to create a c_char_p string that can be passed to Fortran. I admit, this bit of code is the result of some trial and error to figure out what would work (the documentation is less than helpful and of course doesn't mention Fortran at all). Ctypes requires that the Fortran code be compiled as a shared library (a .DLL on Windows, or .so on Linux). We can then define the interfaces to the routines in the DLL we want to call from Python. As long as we get the types of the arguments right (which can be tricky for strings as will be shown) it works great. Consider this example:
# create a dict:
a = {'Generated in Python': True,
'scalar': 1,
'vector': [1,2,3],
'string': 'hello'}
cp = python_dict_to_fortran(a) # creates the c_char_p
dll = CDLL ('lib/json.so') # load the shared library
# define the interface to one of the Fortran routines:
send_json_to_fortran = dll.send_json_to_fortran
send_json_to_fortran.argtypes = [POINTER(c_char_p)]
send_json_to_fortran.restype = None
send_json_to_fortran(byref(cp)) # send the string to fortran
Here we are creating a dictionary, converting it into a string, loading our Fortran DLL, defining the interface to the Fortran function we want to call (send_json_to_fortran()), and then calling it. On the Fortran side, we first define a routine to convert the input string into a JSON structure (using JSON-Fortran):
subroutine c_ptr_to_json(cp,json)
type(c_ptr),intent(in) :: cp
!! a `c_char_p` from python containing a JSON string.
type(json_file),intent(inout) :: json
!! the JSON data structure
character(len=:),allocatable :: fstr
!! string containing the JSON data
if (c_associated(cp)) then
fstr = c_ptr_to_f_string(cp)
call json%load_from_string(fstr) ! parse JSON
deallocate(fstr)
end if
end subroutine c_ptr_to_json
Note that we are using the following routine to convert the Python c_char_p string into a normal Fortran string:
function c_ptr_to_f_string(cp) result(fstr)
type(c_ptr),intent(in) :: cp !! a `c_char_p` from python
character(len=:),allocatable :: fstr
!! the corresponding fortran string
integer :: ilen !! string length
ilen = strlen(cp) ! C library function
block
! convert the C string to a Fortran string
character(kind=c_char,len=ilen+1),pointer :: s
call c_f_pointer(cp,s)
fstr = s(1:ilen)
nullify(s)
end block
end function c_ptr_to_f_string
The example function we are calling from Python is defined below. It is important to note that Fortran is expecting the cp argument to be passed by reference. That is why we use byref(cp) when we call it from Python.
subroutine send_json_to_fortran(cp) &
bind (c, name='send_json_to_fortran')
type(c_ptr),intent(in) :: cp
!! a `c_char_p` from python containing a JSON string.
type(json_file) :: json
call c_ptr_to_json(cp,json)
! do something with the data
! (in this case, just print it)
call json%print_file()
call json%destroy() ! free memory
end subroutine send_json_to_fortran
This prints the JSON structure (from Fortran) like so:
{
"Generated in Python": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello"
}
Note that or send_json_to_fortran() routine has the BIND(c) attribute to make it callable from C (or in this case, Python). Once we have the JSON data on the Fortran side, we can also access the variables, query what's in the structure, get the values, etc., like we can any JSON data. If you want to modify the data on the Fortran side, keep reading.
Sending a JSON structure from Fortran to Python
The inverse operation (sending arbitrary data using JSON from Fortran back to Python) is also possible, but a little bit more complicated. First, in Fortran, we will define a container type that will hold a deferred-length (allocatable) string:
type,public :: container
character(len=:),allocatable :: str
end type container
The reason we need this is so that we can point to it with a pointer, which will be come apparent shortly. We then define a Fortran function to convert a JSON structure into a string, and then allocate a pointer to a container to hold the string:
subroutine json_to_c_ptr(json,cp,destroy)
implicit none
type(json_file),intent(inout) :: json !! JSON data
type(c_ptr) :: cp
!! a pointer to a container
!! containing the JSON data as a string
logical,intent(in) :: destroy
!! to also destroy the JSON file
!! (must be destroyed on the fortran
!! side somewhere to prevent memory leak)
character(len=:),allocatable :: str
type(container),pointer :: c
call json%print_to_string(str)
if (destroy) call json%destroy()
allocate(c)
c%str = str
cp = c_loc(c)
end subroutine json_to_c_ptr
We also will need a routine that is callable from Python to get the length of a string within a container:
function get_string_length(cp) result(ilen) &
bind (c, name='get_string_length')
type(c_ptr),intent(in) :: cp !! pointer to a container
integer(c_int) :: ilen !! the length of the string
type(container),pointer :: c
ilen = 0
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) ilen = len(c%str)
end if
end function get_string_length
Now, we also need another routine callable from Python to populate a c_char_p string buffer:
subroutine populate_char_string(cp,string) &
bind(c,name='populate_char_string')
type(c_ptr),intent(in) :: cp !! pointer to a container
type(c_ptr),intent(inout) :: string
!! a preallocated string buffer that
!! the string will copied into
type(container),pointer :: c
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) then
call f_string_to_c_ptr(c%str, string)
end if
end if
end subroutine populate_char_string
Finally, we also need a way to destroy the string that we are creating in Fortran. This has to be done in Fortran, since the container pointer was allocated in Fortran. Python will garbage collect its own variables, but it doesn't know how to garbage collect a Fortran pointer, so we'll have to call this routine from Python when we no longer need the variable.
subroutine destroy_string(cp) &
bind (c, name='destroy_string')
type(c_ptr),intent(inout) :: cp
!! pointer to a container
type(container),pointer :: c
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) deallocate(c%str)
deallocate(c)
end if
cp = c_null_ptr
end subroutine destroy_string
Now, say we have the following Fortran routine that generates some JSON data and sends it to Python:
subroutine get_json_from_fortran(cp) &
bind (c, name='test_send_json_to_python')
type(c_ptr) :: cp
!! pointer to a container containing a json string
type(json_file) :: json
! sample data:
call json%add('Generated in Fortran', .true.)
call json%add('scalar', 1)
call json%add('vector', [1,2,3])
call json%add('string', 'hello')
call json%add('string array', ['1','2','3'])
! convert it to a c_ptr (and destroy JSON structure)
call json_to_c_ptr(json,cp,destroy=.true.)
end subroutine get_json_from_fortran
On the Python side, we have:
# declare some more interfaces to the Fortran routines:
get_string_length = dll.get_string_length
get_string_length.argtypes = [POINTER(c_void_p)]
get_string_length.restype = c_int
populate_char_string = dll.populate_char_string
populate_char_string.argtypes = [POINTER(c_void_p),POINTER(c_char_p)]
populate_char_string.restype = None
destroy_string = dll.destroy_string
destroy_string.argtypes = [POINTER(c_void_p)]
destroy_string.restype = None
def fortran_str_to_python_dict(cp,destroy=True):
# get the length of the string:
length = c_int()
length = get_string_length(byref(cp))
# preallocate a string buffer of the correct size to hold it:
s = c_char_p()
s = cast(create_string_buffer(b' '.ljust(length)),c_char_p)
# convert it to a normal python string:
populate_char_string(byref(cp),byref(s))
string = s.value.decode()
if (destroy):
# destroy it on the Fortran side:
destroy_string(byref(cp))
return json.loads(string)
The fortran_str_to_python_dict() routine is the main routine that converts the string that comes from Fortran into a Python dictionary. First we get the string length, then allocate a string buffer to hold it, we populate the buffer, and then parse the JSON. Note that we also have the option to destroy the Fortran string (which if you recall was allocated as a pointer in json_to_c_ptr() so it needs to be deallocated to prevent a memory leak).
Now we can use these routines from Python like so:
get_json_from_fortran = dll.get_json_from_fortran
get_json_from_fortran.argtypes = [POINTER(c_void_p)]
get_json_from_fortran.restype = None
cp = c_void_p()
get_json_from_fortran(byref(cp))
d = fortran_str_to_python_dict(cp,destroy=True)
print(json.dumps(d, indent=2))
Which prints the result:
{
"Generated in Fortran": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello",
"string array": [
"1",
"2",
"3"
]
}
From Python to Fortran and Back
Finally, another use case would be to generate some data in Python, pass it to Fortran where it is modified, and then return it back to Python. In this case, we need to use the container approach, with a few additional routines. On the Python side:
c_ptr_to_container_c_ptr = dll.c_ptr_to_container_c_ptr
c_ptr_to_container_c_ptr.argtypes = [POINTER(c_char_p),POINTER(c_void_p)]
c_ptr_to_container_c_ptr.restype = None
def python_str_to_container(s):
cp = python_str_to_fortran(s)
ccp = c_void_p()
c_ptr_to_container_c_ptr(byref(cp),byref(ccp))
return ccp
def python_dict_to_container(d):
return python_str_to_container(json.dumps(d))
And we also need a routine on the Fortran side to convert from the string buffer format to the container format:
subroutine c_ptr_to_container_c_ptr(cp,ccp) &
bind (c, name='c_ptr_to_container_c_ptr')
type(c_ptr),intent(in) :: cp
!! a `c_char_p` from python
type(c_ptr),intent(out) :: ccp
!! pointer to a container that contains the string
character(len=:),allocatable :: str
!! fortran version of the string
type(container),pointer :: c
!! container to hold the string
if (c_associated(cp)) then
! get the fortran string:
str = c_ptr_to_f_string(cp)
! return a pointer to the container:
allocate(c)
c%str = str
ccp = c_loc(c)
else
ccp = c_null_ptr
end if
end subroutine c_ptr_to_container_c_ptr
Here is an example use case:
subroutine modify_json_in_fortran(cp) &
bind (c, name='modify_json_in_fortran')
type(c_ptr),intent(inout) :: cp !! a pointer to a container
type(json_file) :: json
type(container),pointer :: c
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) then
! parse JSON:
call json%load_from_string(c%str)
!do something with the data:
call json%print_file()
call json%add('Added in Fortran', [9,10])
! convert it to a c_ptr (and destroy JSON structure)
call json_to_c_ptr(json,cp,destroy=.true.)
end if
end if
end subroutine modify_json_in_fortran
Which we can call from Python like so:
modify_json_in_fortran = dll.modify_json_in_fortran
modify_json_in_fortran.argtypes = [POINTER(c_void_p)]
modify_json_in_fortran.restype = None
cp = python_dict_to_container(a)
modify_json_in_fortran(byref(cp))
d = fortran_str_to_python_dict(cp,destroy=True)
print('')
print('modified by Fortran:')
print(json.dumps(d, indent=2))
Which prints:
{
"Generated in Python": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello"
}
modified by Fortran:
{
"Generated in Python": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello",
"Added in Fortran": [
9,
10
]
}
Summary
Once all the interface routines are created, this technique is fairly easy to use and seems to work great. It is probably fine for a lot of use cases where extreme performance is not required. For example, it may be a good enough technique for exchanging data from a Python GUI and a core application written in Fortran. However, for very large data sets, this may not be very efficient since it does involve converting everything to and from strings.
See also
Updated on 7 June 2018: the code is now available on GitHub.
