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]
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.
The upcoming Fortran standard formerly known as Fortran 2015 has a new name: Fortran 2018. It was decided to change it in order to match the expected year of publication. This makes sense. The previous standard (Fortran 2008) was published in 2010.
Waiting for an updated Fortran standard is an exercise in Zen-like patience. Almost a decade after Fortran 2008, we'll get a fairly minor update to the core language. And it will be years after that before it's fully supported by any compiler that most users will have available (gfortran still doesn't have a bug-free implementation of all of Fortran 2008 or even 2003). Fortran was essentially reborn in the Fortran 2003 standard, which was an amazing update that brought Fortran into the modern world. It's a terrific programming language for scientific and technical computing. However, the limitations are all too clear and have been for a long time:
- We need better facilities for generic programming. It's impossible to do some things without having to duplicate code or use "tricks" like preprocessing or include files (see JSON-Fortran for examples).
- We need some kind of exception handling. Fortran does have a floating-point exception handling feature, but honestly, it's somewhat half-baked.
- We need a better implementation of strings. Allocatable strings (introduced in Fortran 2003) are great, but not enough, since they can't be used in all instances where strings are needed.
- We need the language to be generally less verbose. I'm tired to having to type multiple nested SELECT TYPE statements to do something that is a one liner in Python (sure I know Fortran will never be as succinct as Python, but some of the verbosity required for object-oriented Fortran is just perverse).
- We need any number of new features to make it easier to extend the language with third-party libraries (so we don't have to wait two decades for a feature we want).
- We also need the language development process to embrace a more open collaborative model using modern tools (Usenet is not the future). I guess the recent survey was unprecedented, but it's not enough.
Fortran is a programming language that needs a better PR department. Legacy Fortran codes are being rewritten in C++, Python, or even Julia. NASA frequently throws massive Fortran 77 libraries with decades of heritage (e.g., DPTRAJ/ODP, GTDS, SPICELIB) into the trash in order to rewrite it all from the ground up in C++, without ever considering Fortran 2003+ (or maybe not realizing it exists?). The information about modern Fortran on the internet is spotty at best, outdated, or downright wrong (what is the deal with REAL*8?). In popular consciousness Fortran is mostly a punchline (usually something to do with punchcards and your granddad). A language like Python (which was never designed for technical computing) is now seen by many as a superior solution for technical computing. Matlab refers to itself as "the only top programming language dedicated to mathematical and technical computing"! The Julia website lists somewhat misleading benchmarks than implies that C, Julia, and even Lua are faster than Fortran.
Now, get off my lawn, you kids!
See also
Python's matplotlib.pyplot is a very nice collection of functions that provide an easy Matlab-like interface for data plotting. It can be used to generate quite professional looking plots. There is a lot of information on the internet about calling Fortran from Python, but what if all you want to do is generate some plots from your (modern) Fortran program? With this in mind, I've created a very simple Fortran wrapper for Pyplot, allowing you to make pretty good plots by writing only Fortran code. Consider the following example:
program test
use,intrinsic :: iso_fortran_env, only: wp => real64
use pyplot_module
implicit none
real(wp),dimension(100) :: x,sx,cx,tx
type(pyplot) :: plt
integer :: i
!generate some data:
x = [(real(i,wp), i=0,size(x)-1)]/5.0_wp
sx = sin(x)
cx = cos(x)
tx = sx * cx
!plot it:
call plt%initialize(grid=.true.,xlabel='angle (rad)',&
title='python_plot test',legend=.true.)
call plt%add_plot(x,sx,label='\$\sin (x)\$',&
linestyle='b-o',markersize=5,linewidth=2)
call plt%add_plot(x,cx,label='\$\cos (x)\$',&
linestyle='r-o',markersize=5,linewidth=2)
call plt%add_plot(x,tx,label='\$\sin (x) \cos (x)\$',&
linestyle='g-o',markersize=2,linewidth=1)
call plt%savefig('test.png')
end program test
The main user interface is the pyplot class, which has methods such as initialize, add_plot, and savefig. This example produces the following plot:
All the module does is generate a Python script, and then runs it and saves the result. The Python call is completely transparent to the Fortran user.
I posted the project to GitHub (pyplot-fortran). Maybe someone else will find it useful, or help to expand it.
See also
- matplotlib
- F2PY Users Guide and Reference Manual --Fortran to Python interface generator
There is a lot of confusion and misinformation about the Fortran programming language on the internet, and a general ignorance about it among programmers. Most younger programmers who use languages invented five minutes ago probably have never seen it, and may only be dimly aware of it as some obsolete language that nobody uses anymore.
You may be surprised to learn that Fortran is a modern, object-oriented, general-purpose programming language. Fortran is not a programming language created by computer scientists to write operating systems, nor does it include every single programming concept anyone's ever heard of. It is a language designed for computational efficiency and efficient array manipulation, and has a clear uncluttered syntax that can be read and understood by non-experts with only a little effort and training. It is particularly suited for numerical and scientific programming by non-expert programmers such as engineers and scientists. Learning all of Fortran is considerably easier than learning all of C++ (as an example). Sure, you can't do template metaprogramming, but few engineer/scientist types would ever want to do that anyway (besides, it's bad for you and could make you go blind).
By "Fortran", I mean modern Fortran (i.e., Fortran 2003 and 2008, which is the latest standard). Yes, the roots of Fortran go way back to the 1950s. Sure, early Fortran programs were written on punched cards. So what? Latin was once scratched into wax tablets, but that isn't really relevant to modern Italian and French speakers. In fact, the Fortran language has evolved considerably since it was first standardized in 1966. It generally has followed a cycle where a major update is followed by a minor update (1977=minor, 1990=major, 1995=minor, 2003=major, 2008=minor). It has been said that the 2003 update was as big an update to Fortran 95 as C++ was to C! Ten years later, the GNU Fortran compiler is still not fully F2003 compliment (the Intel compiler only recently became so).
People who attempt to enter the world of Fortran programming are easily corrupted and discouraged by misinformation. Even a lot of old-school Fortran users are unaware of the later standards. This is too bad, because modern Fortran is actually quite a respectable programming language for a lot of technical applications. This article is a pretty good overview of Fortran for C/C++ programmers. However, it is outdated, since it is confined to Fortran 95. Most of the limitations it mentions (no procedure pointers, clunky character strings, lack of an intent attribute for pointer dummy arguments, the nonstandardness of the ; character) have been rectified in subsequent standards.
The fact is the internet is not really the best source of information for modern Fortran. One day, maybe, there will be vibrant community of Fortran users on the internet, extensive online documentation, open source projects, and all your questions will simply be a web search away (cf., Python). But for now, you'll probably have to buy some books. If a book has the numbers 77, 90, or 95 in the title, don't open it, it will only confuse you. This is not to say that there aren't friendly Fortran folk on the internet who will also help you out. Two of the best places to go with questions are the Intel Fortran forum and the comp.lang.fortran newsgroup (yes, apparently, Usenet still exists).
References
- N. Maclaren, Why (and Why Not) to Use Fortran: Instead of C++, Matlab, Python etc., University of Cambridge Computing Service, June 2012.
- L. Phillips, Scientific computing’s future: Can any coding language top a 1950s behemoth?, May 7, 2014 [arstechnica.com]
- Fortran Wiki -- an open venue for discussing all aspects of the Fortran programming language and scientific computing.
- A. Koenig, C Traps and Pitfalls, AT&T Bell Laboratories.
- Why C and C++ are Awful Programming Languages
I just started a new modern Fortran software library called bspline-fortran, which is for multidimensional (multivariate) b-spline interpolation of data defined on a regular grid. It is available on GitHub, and released under a permissive BSD-style license.
It seems impossible to find code for higher than 3D spline interpolation on the internet. Einspline only has 1D-3D, as do the NIST Core Math Library DBSPLIN and DTENSBS routines. Python's SciPy stops at 2D (Bivariate splines). At first glance, the SciPy routine map_coordinates seems to do what I want, but not really since there doesn't seem to be any way to specify the abscissa vectors (x, y, z, etc.). Note: maybe there is a way, but looking at the less-than-great documentation, the minimally-commented C code, and then reading several confusing Stack Overflow posts didn't help me much. Apparently, I'm not the only one confused by this routine.
After much searching, I finally just took the 2D and 3D NIST routines (written in 1982), refactored them significantly so I could better understand what they were doing, and then extended the algorithm into higher dimensions (which is actually quite easy to do for b-splines). The original FORTRAN 77 code was somewhat hard to follow, since it was stuffing matrices into a giant WORK vector, which required a lot of obscure bookkeeping of indices. For example, I replaced this bit of code:
IZM1 = IZ - 1
KCOLY = LEFTY - KY + 1
DO 60 K=1,KZ
I = (K-1)*KY + 1
J = IZM1 + K
WORK(J) = DBVALU(TY(KCOLY),WORK(I),KY,KY,IDY,YVAL,INBV1,WORK(IW))
60 CONTINUE
with this:
kcoly = lefty - ky + 1
do k=1,kz
temp2(k) = dbvalu(ty(kcoly:),temp1(:,k),ky,ky,idy,yval,inbv1,work)
end do
which makes it a lot easier to deal with and add the extra dimensions as necessary.
The new library contains routines for 2D, 3D, 4D, 5D, and 6D interpolation. It seems like someone else may find them useful, I don't know why they don't seem to be available anywhere else. Eventually, I'll add object-oriented wrappers for the core routines, to make them easier to use.
I just came across this article from 2005 written by Jef Raskin (probably most famous for having initiated the Macintosh project at Apple) on source code documentation [1]. I agree with much of what he recommends. I have never believed that code can be "self-documenting" without comments, and generally will start writing the comments before I start the actual code. I also like some aspects of Knuth's "literate programming" concept [2].
I do find automatic documentation generators useful, however. There aren't a lot of solutions for modern Fortran projects though. I've not tried to use Doxygen, but am told that support for modern Fortran features is somewhat lacking [3]. ROBODoc is pretty good, although it does not perform any code analysis, and only extracts specific comment blocks that you indicate in the source. A brand new tool called ford is very promising, since it was created with modern Fortran in mind.
I also have an idea that the best code (including comments) should be understandable to a reasonably intelligent person who is not familiar with the programming language. Some languages (Fortran, Python) make this easier than others (C++, Haskell). There is a good article here, discussing the "beauty" of some C++ code (although I disagree with his views on comments). I was struck by the discussion on how to write functions so that more information is provided to the reader of the code [4]. The example given is this:
int idSurface::Split( const idPlane &plane, const float epsilon,
idSurface **front, idSurface **back,
int *frontOnPlaneEdges, int *backOnPlaneEdges ) const;
Of course, if you are not familiar with the C++ *, **, and & characters and the const declaration, there is no information you can derive from this. It is complete gibberish. Modern Fortran, on the other hand, provides for the declaration of an intent attribute to function and subroutine arguments. The equivalent Fortran version would look something like this:
integer function split(me,plane,eps,front,back,&
frontOnPlaneEdges,backOnPlaneEdges)
class(idSurface),intent(in) :: me
type(idPlane), intent(in) :: plane
real, intent(in) :: eps
type(idSurface), intent(out) :: front
type(idSurface), intent(out) :: back
integer, intent(out) :: frontOnPlaneEdges
integer, intent(out) :: backOnPlaneEdges
Here, it should be more obvious what is happening (for example, that plane is an input and front is an output). The use of words rather than symbols definitely improves readability for the non-expert (although if taken to the extreme, you would end up with COBOL).
References
- J. Raskin, "Comments are More Important than Code", Queue, Volume 3 Issue 2, March 2005, Pages 64-65.
- D. E. Knuth, "Literate Programming", The Computer Journal (1984) 27 (2): 97-111.
- F03/08 supporting Documentation tools, comp.lang.fortran.
- S. McGrath, "The Exceptional Beauty of Doom 3's Source Code", kotaku.com, Jan. 14, 2013.
