Degenerate Conic

Xojo and Memory Lane

I happened across this article on a programming language and IDE called Xojo. I never used it under this name, but I have fond memories of using it when it was called RealBASIC. I learned it shortly after it came out in the late 1990s. I found it very intuitive and it was super easy to create fairly nice looking GUI applications for the Mac. In fact, the first optimization code I ever wrote was in RealBASIC, for solving the traveling salesman problem.

In the mid-2000's when I was a grad student, I used it to build a GUI interface to JPL's BlackJack GPS receiver. This included reading data from the receiver, as well as sending commands to it. There were numerous floating windows that displayed data and plots (see screenshot at right). At the time, JPL's interface code only ran on pre-MacOS X systems, which were rapidly becoming obsolete. My RealBASIC app ran great on MacOS X. In fact, I believe JPL requested a copy when they saw it (I have no idea if it lives on somewhere at JPL).

In 2006, I traveled to Germany with a blue clamshell iBook to deliver some updated software to the TerraSAR-X satellite, which included the same type of GPS receiver. I remember during the update process, something unexpected occurred (the details escape me now). In any event, it was required to generate some binary command files to upload to the system in order to set things right. My RealBASIC interface app was already set up to do that, so I fired up my iBook and was able to generate the files in short order. I distinctly remember some surprise from the German engineers that I could do such a thing with this silly looking machine! After uploading the files, things got back on track. TerraSAR-X has been in orbit since June 2007.

I haven't used it in years, but I'm glad to know that RealBASIC is still around.

See also

• B. Keeney, Xojo: the Best Secret in the Programming Industry Part 1, October 19, 2016
• O. Montenbruck, M. Garcia-Fernandez, and J Williams, Performance comparison of semicodeless GPS receivers for LEO satellites, GPS Solut (2006) 10: 249–261.

Flang

It appears that Flang, a new open source Fortran front-end for LLVM, has appeared on GitHub recently with little fanfare. This is apparently the result of NVIDIA's previously-announced plan to open source the PGI Fortran compiler. Unfortunately, they decided to give it the same name as another earlier attempt to create a Fortran/LLVM compiler (more confusion for poor Fortran programmers). I don't really know how it compares to Gfortran or Intel (PGI appears to be lagging behind on support for the Fortran 2008 standard). Initial tests by Usenet denizens (yes Usenet still exists) indicate that maybe Flang isn't quite ready for prime time. Hopefully it will improve with time. I think it's great news to potentially have another free Fortran compiler available.

See also

• "Flang: Developing an open source Fortran front-end for LLVM", ARM HPC User Group 2016.
• Fortran + LLVM [Earlier post]

Fortran Configuration File Formats

String and file manipulation in Fortran isn't as bad as you've heard (OK, it's bad, but it's getting better). Sure, modern Fortran only provides the bare minimum of features for this sort of thing out of the box, but various libraries are starting to appear that the modern Fortran programmer can use to deal with a range of different file formats including JSON, INI, CSV, XML, etc. The Fortran user community has been late to the game on this opensource thing. But it does exist and I encourage all Fortraners to use it. You don't have to roll your own file format or write your own parser. If the existing libraries don't have what you need, then help to improve them.

The following is a look at a few different options for configuration files in your Fortran program. It is limited to text files. For binary files, I would highly recommend a standard cross-platform format (such as HDF5). For goodness' sake, do not invent your own binary file format.

Namelists

Fortran namelists are very old-school, but are still part of the standard (indeed this is the only high-level file read/write feature that is built into Fortran). Namelist have a few advantages, but they are generally outweighed by their many disadvantages. Some advantages include:

• It is super easy to use. The reading and writing of the file is automatically handled by the compiler. A namelist is declared like so:

integer :: a
real :: b
character(len=5) :: c
integer :: iunit,istat

namelist /blah/ a,b,c

!set values:
a = 1
b = 2.0
c = 'three'

and can be written by:

write(unit=iunit, nml=blah, iostat=istat)

Which produces a file like this:

&BLAH
A= 1,
B= 2.00000000 ,
C="three",
/

The file can be read by:

read(unit=iunit, nml=blah, iostat=istat)

• It automatically handles any Fortran variables you throw at it (including multidimensional arrays and derived data types).
• Variables in the file can be optional. Any variable not present in the file retains the value it had before the namelist was read.

Disadvantages include:

• Different compilers can and will output slightly different formats, and you have absolutely no control over that.
• If you want to generate or manipulate your input files using a different language, you'll probably have to write your own parser to do it. Hardly any other programming language will have a decent namelist parser or writer available. A notable exception is f90nml, a great Python library for reading and writing namelists.
• The variables in the file must exactly correspond (in name, type, and size) to variables in the code. This turn out to be very restrictive for a lot of reasons. For one thing, it is quite annoying for strings, since it requires you to specify a maximum string length in the variable declaration.
• It can be difficult to diagnose errors when reading the file (e.g, a syntax error, an unrecognized variable, or an array size or string length larger than the Fortran variable). The compiler will usually just throw an I/O error if anything goes wrong. However, it is possible on some compilers to retrieve the line where the error occurred by rewinding one line and then reading it with a normal READ statement.
• It is difficult to tell if a variable is actually present in the file (if one isn't there, then the corresponding variable in the code will just retain whatever value it had before the read was attempted. Thus, the only real way to check if a variable was present or not is to initialize them all to unusual values, and then check for this after the read.)

In my view, there is little reason to start using namelists at this point, since there are now better options available. Unless you have a legacy code that is using them which can't be changed, I would recommend using something else.

JSON

JSON stands for JavaScript Object Notation, and is a very popular and lightweight data-interchange format. There is an interface for practically all known programming languages, including modern Fortran. My JSON-Fortran library was built on an older Fortran 95 library, which I forked in order to be able to use some of the newer Fortran 2003/2008 features. With JSON-Fortran, the example above would be:

type(json_file) :: json
call json%add('blah.a',1)
call json%add('blah.b',2.0)
call json%add('blah.c','three')

Which produces the JSON file:

{
    "blah": {
        "a": 1,
        "b": 2.0,
        "c": "three"
    }
}

Note that f90nml can also be used to transition away from namelists by converting them to JSON using Python. To go from a namelist to JSON, all you have to do is this:

import f90nml
import json
n = f90nml.read('my_namelist.nml')
with open('my_namelist.json', 'w') as outfile:
json.dump(n, outfile, sort_keys = True, indent = 4)

Using my JSON-Fortran library, you can even convert a JSON file back into a Namelist if that sort of thing is your bag (as shown in a previous post). One disadvantage of JSON is that the standard does not allow for comments in the file. However, some parsers include support for comments (including JSON-Fortran).

INI

INI is a simple, easy-to-write format. Although not really standardized, it is based on the old MS-DOS initialization files. For modern Fortran, the FiNeR library provides a Fortran interface. The INI version of the above file would be:

[blah]
a=1
b=2.0
c=three

Here's an example of converting an INI file to a JSON file (using FiNeR and JSON-Fortran):

subroutine ini2json(ini_filename, json)

use json_module
use finer

implicit none

character(len=*),intent(in) :: ini_filename
type(json_file),intent(out) :: json

type(file_ini) :: ini
integer :: ierr
integer :: i
character(len=:),dimension(:),allocatable :: secs
character(len=:),dimension(:),allocatable :: opts

call ini%load(filename=ini_filename,error=ierr)
call ini%get_sections_list(secs)
do i=1,size(secs)
    do
        if (.not. ini%loop(trim(secs(i)),opts)) exit
        call json%add(trim(secs(i))//'.'&
                    trim(opts(1)),trim(opts(2)))
    end do
end do
call ini%free()

end subroutine ini2json

You could call this routine like so:

program main
use ini2json_module
use json_module
implicit none
type(json_file) :: json
call ini2json('test.ini',json)
call json%print_file()
call json%destroy()
end program main

And it would produce the following result for the file above:

{
    "blah": {
        "a": "1",
        "b": "2.0",
        "c": "three"
    }
}

Note that, in this example, all the variables are returned as strings (getting the numerical values as numbers is left as an exercise to the reader).

Others

• XML (Extensible Markup Language) -- I never really liked this one because it is not easily human readable or editable. There are some Fortran XML parsers out there, though.
• YAML (YAML Ain't Markup Language) -- I'm not as familiar with this one (also I hate recursive acronyms). There is a Fortran interface to it called fortran-yaml, but I haven't tried it.
• TOML (Tom's Obvious, Minimal Language) - I'm not aware of a Fortran parser for this one. This one was created by one of the founders of GitHub.
• Lua - Another interesting possibility is for your configuration file to be written in a full-up scripting language. There are a couple of Lua-Fortran interfaces out there (e.g., AOTUS).

See also

• Fortran 90: A Conversion Course for Fortran 77 Programmers, Section 6.3 NAMELIST. "The NAMELIST statement has been available as a suppliers extension to Fortran since the early days (it was available as an IBM extension to FORTRAN II in the early 60's!). It has now been included in the Fortran 90 language. However, NAMELIST is a poorly designed feature and should be avoided whenever possible."
• Parsing INI configuration files with Fortran, August 23, 2011 [shocksolution.com]
• YAML, JSON and Human Readable Configuration Files [www.marshallward.org]
• FiNeR: Fortran INI ParseR and generator [GitHub]
• JSON-Fortran: A Fortran 2008 JSON API [GitHub]
• JSON as configuration files: please don’t

GFortran 7.1

GCC (which includes GFortran) version 7.1 has just been released. Various features from Fortran 2003-2015 have been added, including:

• User-defined derived-type input/output (UDTIO).
• Partial support for derived type coarrays with allocatable and pointer components.
• Non-constant stop codes and error stop codes.
• Derived types with allocatable components of recursive type.
• Intrinsic assignment to polymorphic variables.

Partial support for OpenMP 4.5 is also included. It looks like they also added support for several non-standard extensions that are present in legacy code. GFortran still hasn't quite made it to full Fortran 2003 support, over a decade after that standard was published. But they are getting there. Thanks to the small number of volunteers who work on this. You are the real heroes!

See also

• What's new in gfortran?
• Compiler Support for Modern Fortran [Fortran Wiki]

JSON-Fortran 5.3.0

JSON-Fortran 5.3.0 is out. This release includes a few minor bug fixes and a couple of new features. For one thing, it is now possible to easily build or modify a JSON structure without using pointer variables. Here's an example:

program test

use json_module, rk=>json_rk

implicit none

type(json_file) :: f

call f%initialize()

call f%add('name', 'test')
call f%add('inputs.t', 0.0_rk)
call f%add('inputs.x', [1.0_rk,2.0_rk,3.0_rk])
call f%add('outputs.t', 10.0_rk)
call f%add('outputs.x(1)', 11.0_rk)
call f%add('outputs.x(2)', 22.0_rk)
call f%add('outputs.x(3)', 33.0_rk)

call f%print_file()

call f%destroy()

end program test

As shown, the structure can be built by specifying the variable paths. This way, objects, vectors, vector elements, etc. can all be added. The result of this example is:

{
    "name": "test",
    "inputs": {
        "t": 0.0E+0,
        "x": [
            0.1E+1,
            0.2E+1,
            0.3E+1
        ]
    },
    "outputs": {
        "t": 0.1E+2,
        "x": [
            0.11E+2,
            0.22E+2,
            0.33E+3
        ]
    }
}

So now, in addition to the previous methods of building a JSON structure (reading it from a file, reading it from a string, or building it yourself using pointer variables), there is now another way. Each method has its uses for different applications. Building it from the path strings can be very convenient and also opens up the possibility of easy conversion between a Fortran namelist and a JSON file. I'll show an example of this in a future post. In a previous post, I showed how to do the reverse (JSON to namelist).

JSON-Fortran is my most popular GitHub project (70 stars as of today). It usually makes it to GitHub's monthly "trending" Fortran projects page. I'm always interested to hear about what people are using it for. I added a page on the wiki to list some of the projects that use it, so if you're using it, feel free to add your project to the list.

Binary Search

Binary search is a search algorithm that finds the position of a target value within a sorted array. The following binary search routine is taken from NASTRAN, a venerable finite element analysis program from NASA:

    SUBROUTINE BISLOC (*,ID,ARR,LEN,KN,JLOC)
C-----
C BINARY SEARCH - LOCATE KEY WORD 'ID' IN ARRAY 'ARR', 1ST ENTRY
C IF FOUND, 'JLOC' IS THE MATCHED POSITION IN 'ARR'
C IF NOT FOUND, NON-STANDARD RETURN
C I.E.
C ID = KEY WORD TO MATCH IN ARR. MATCH AGAINST 1ST COL OF ARR
C ARR = ARRAY TO SEARCH. ARR(ROW,COL)
C LEN = LENGTH OF EACH ENTRY IN ARRAY. LEN=ROW
C KN = NUMBER OF ENTRIES IN THE ARR. KN =COL
C JLOC= POINTER RETURNED - FIRST WORD OF ENTRY. MATCHED ROW
C-----
C
    INTEGER ARR(1)
    DATA ISWTCH / 16 /
C
    JJ = LEN - 1
    IF (KN .LT. ISWTCH) GO TO 120
    KLO = 1
    KHI = KN
10  K = (KLO+KHI+1)/2
20  J = K*LEN - JJ
    IF (ID-ARR(J)) 30,90,40
30  KHI = K
    GO TO 50
40  KLO = K
50  IF (KHI-KLO -1) 100,60,10
60  IF (K .EQ. KLO) GO TO 70
    K = KLO
    GO TO 80
70  K = KHI
80  KLO = KHI
    GO TO 20
90  JLOC = J
    RETURN
100 JLOC = KHI*LEN - JJ
    J = KN *LEN - JJ
    IF (ID .GT.ARR(J)) JLOC = JLOC + LEN
110 RETURN 1
C
C SEQUENTIAL SEARCH MORE EFFICIENT
C
120 KHI = KN*LEN - JJ
    DO 130 J = 1,KHI,LEN
    IF (ARR(J)-ID) 130,90,140
130 CONTINUE
    JLOC = KHI + LEN
    GO TO 110
140 JLOC = J
    GO TO 110
    END

As you can see, it is a scene of FORTRAN 66 horrors! However, it can be refactored into a nice modern Fortran routine like so:

pure function bisloc(id,arr) result(jloc)

!! binary search of a sorted array.

implicit none

integer,intent(in) :: id
    !! key word to match in `arr`
integer,dimension(:),intent(in) :: arr
    !! array to search (it is
    !! assumed to be sorted)
integer :: jloc
    !! the first matched index in 'arr'
    !! (if not found, 0 is returned)

integer :: j,k,khi,klo,n
integer,parameter :: iswtch = 16

n = size(arr)
jloc = 0

if ( n<iswtch ) then

    ! sequential search more efficient

    do j = 1 , n
        if ( arr(j)==id ) then
            jloc = j
            return
        else if (arr(j)>id) then
            return ! error
        end if
    end do
    return ! error

else

    klo = 1
    khi = n
    k = (klo+khi+1)/2
    do
        j = k
        if ( id<arr(j) ) then
            khi = k
        else if ( id==arr(j) ) then
            jloc = j
            return
        else
            klo = k
        end if
        if ( khi-klo<1 ) then
            return ! error
        else if ( khi-klo==1 ) then
            if ( k==klo ) then
                k = khi
            else
                k = klo
            end if
            klo = khi
        else
            k = (klo+khi+1)/2
        end if
    end do

end if

end function bisloc

I did change a couple of things:

• The original routine was using an obsolete alternate return (note the * argument) for an error condition. The new one just returns 0 if there is an error.
• The original routine was passing in the array as a vector and treating it like a matrix (it's either searching a row or column of this matrix, I admit I had a hard time telling). The new version assumes the input is a rank-1 array. You can pass in whatever matrix slices you want (say a column: a(:,1) or a row: a(1,:)).
• I also made it a pure function, because that's how I roll.

Note that NASTRAN was released under the NASA Open Source Agreement version 1.3, a somewhat unusual open source license. It isn't considered a "free" license by the FSF since it includes a statement requiring changes to be your "original creation".

Unfortunately, the trend at NASA (and elsewhere) is to declare old Fortran code like this to be totally unmaintainable. So, instead of refactoring them into modern Fortran standards, they are thrown away and rewritten from scratch using a more "modern" language like C++ (whose roots go back to 1972!).

See also

• NASTRAN-95 at NASA's GitHub.

Latest Library Updates

I just released some updates to my two most popular Fortran libraries on GitHub: JSON-Fortran and bspline-fortran. Coincidently, both are now at v5.2.0. Details of the updates are:

JSON-Fortran

There are several new features in this release. The biggest update is that now the code can parse JSON files that include comments (it just ignores them). You can even specify the character that identifies a comment. For example, if using # as a comment character, the following file can now be parsed just fine:

{
    "t": 0.12345,  # this is the time
    "x": 123.7173  # this is the state
}

Technically, comments are not part of the JSON standard. Douglas Crockford, the creator of JSON, had his reasons [1] for not including them, which I admit I don't understand (something about parsing directives and interoperability?) I mainly use JSON for configuration files, where it is nice to have comments. Crockford's suggestion for this use case is to pipe your commented JSON files through something called JSMin before parsing, a solution which seems somewhat ridiculous for Fortran users. So, never fear, now we can have comments in our JSON files and continue not using JavaScript for anything.

Another big change is the addition of support for the RFC 6901 "JSON Pointer" path specification [2]. This can be used to retrieve data from a JSON structure using its path. Formerly, JSON-Fortran used a simple path specification syntax, which broke down if the keys contained special characters such as ( or ). The new way works for all keys.

Bspline-Fortran

A minor update to Bspline-Fortran is the addition of an extrapolation mode. Formerly, if an interpolation was requested outside the bounds of the data, an error was returned. Now, the user has the option of enabling extrapolation.

See also

1. D. Crockford, "Comments in JSON", Apr 30, 2012 [Google Plus]
2. JavaScript Object Notation (JSON) Pointer, RFC 6901, April 2013 [IETF]

Fortran CSV Module

The lowly CSV (comma-separated values) file never gets any respect. Sure, it isn't really standardized, but it is a very useful text file format for columns of data, and is frequently encountered in science/engineering fields. Of course, there is no good modern Fortran library for reading and writing CSV files (see [1-2] for some older Fortran 90 code for writing only). Searching the internet (e.g., [3-4]) you will find various ad hoc suggestions for handling these files. It's amazing how little Fortran users expect their language to support non-numerical algorithms. No Fortran programmer would ever suggest that someone roll their own matrix inverse routine, but they don't see anything wrong with having to write their own data file parser (if you told a Python user they had to write their own CSV parser, you'd be laughed out of the place).

So, I created fortran-csv-module (now available on GitHub), a modern Fortran library for reading and writing CSV files. Everything is handled by an object-oriented csv_file class. Here is an example for writing a file:

program csv_write_test

use csv_module
use iso_fortran_env, only: wp => real64

implicit none

type(csv_file) :: f
logical :: status_ok

! open the file
call f%open('test.csv',n_cols=4,status_ok=status_ok)

! add header
call f%add(['x','y','z','t'])
call f%next_row()

! add some data:
call f%add([1.0_wp,2.0_wp,3.0_wp],real_fmt='(F5.3)')
call f%add(.true.)
call f%next_row()
call f%add([4.0_wp,5.0_wp,6.0_wp],real_fmt='(F5.3)')
call f%add(.false.)
call f%next_row()

! finished
call f%close(status_ok)

end program csv_write_test

Which produces the following file:

 x,y,z,t 1.000,2.000,3.000,T 4.000,5.000,6.000,F

Real, integer, logical, or character data can be added as scalars, vectors, and matrices.

When reading a CSV file, the data is stored internally in the class as allocatable character strings, which can be retrieved as real, integer, logical or character vectors as necessary. For example, to get the x, y, z, and t vectors from the previously-generated file:

program csv_read_test

use csv_module
use iso_fortran_env, only: wp => real64

implicit none

type(csv_file) :: f
character(len=30),dimension(:),allocatable :: header
real(wp),dimension(:),allocatable :: x,y,z
logical,dimension(:),allocatable :: t
logical :: status_ok
integer,dimension(:),allocatable :: itypes

! read the file
call f%read('test.csv',header_row=1,status_ok=status_ok)

! get the header and type info
call f%get_header(header,status_ok)
call f%variable_types(itypes,status_ok)

! get some data
call f%get(1,x,status_ok)
call f%get(2,y,status_ok)
call f%get(3,z,status_ok)
call f%get(4,t,status_ok)

! destroy the file
call f%destroy()

end program csv_read_test

Various options are user-selectable for specifying the format (e.g., changing the quote or delimiter characters). You can choose to enclose strings (or all fields) in quotes or not. The library works pretty well, and there are probably additional improvements that could be made. For one thing, it doesn't properly handle the case of a string that contains the delimiter character (I'll eventually fix this). If anybody has any other improvements, fork it and send me a pull request. The license is BSD, so you can use it for whatever you want.

See also

1. J. Burkardt, CSV_IO -- Read and Write Comma Separated Values (CSV) Files. [note: it doesn't actually contain any CSV reading capability]
2. A. Markus, FLIBS -- Includes a csv_file module for writing CSV files.
3. How to read a general csv file [Intel Fortran Forum]
4. Read data from a .csv file in fortran [Stack Overflow]

Multidimensional Linear Interpolation (Part 2)

In an earlier post, I mentioned that we needed an object-oriented modern Fortran library for multidimensional linear interpolation. Well, here it is. I call it finterp, and it is available on GitHub. It can be used for 1D-6D interpolation/extrapolation of data on a regular grid (i.e., not scattered). It has a similar object-oriented interface as my bspline-fortran library, with initialize, evaluate, and destroy methods like so:

type(linear_interp_3d) :: interp
call interp%initialize(xvec,yvec,zvec,fmat,iflag)
call interp%evaluate(x,y,z,f)
call interp%destroy()

For example, the low-level 3D interpolation evaluate routine is:

pure subroutine interp_3d(me,x,y,z,fxyz)

implicit none

class(linear_interp_3d),intent(inout) :: me
real(wp),intent(in) :: x
real(wp),intent(in) :: y
real(wp),intent(in) :: z
real(wp),intent(out) :: fxyz !! Interpolated ( f(x,y,z) )

integer,dimension(2) :: ix, iy, iz
real(wp) :: p1, p2, p3
real(wp) :: q1, q2, q3
integer :: mflag
real(wp) :: fx11, fx21, fx12, fx22, fxy1, fxy2

call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag)
call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag)
call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag)

q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1)))
q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1)))
q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1)))
p1 = one-q1
p2 = one-q2
p3 = one-q3

fx11 = p1*me%f(ix(1),iy(1),iz(1)) + q1*me%f(ix(2),iy(1),iz(1))
fx21 = p1*me%f(ix(1),iy(2),iz(1)) + q1*me%f(ix(2),iy(2),iz(1))
fx12 = p1*me%f(ix(1),iy(1),iz(2)) + q1*me%f(ix(2),iy(1),iz(2))
fx22 = p1*me%f(ix(1),iy(2),iz(2)) + q1*me%f(ix(2),iy(2),iz(2))
fxy1 = p2*( fx11 ) + q2*( fx21 )
fxy2 = p2*( fx12 ) + q2*( fx22 )

fxyz = p3*( fxy1 ) + q3*( fxy2 )

end subroutine interp_3d

The finterp library is released under a permissive BSD-3 license. As far as I can tell, it is unique among publicly-available Fortran codes in providing linear interpolation for up to 6D data sets. There used to be a Fortran 77 library called REGRIDPACK for regridding 1D-4D data sets (it had the option to use linear or cubic interpolation independently for each dimension). Written by John C. Adams at the National Center for Atmospheric Research in 1999, it used to be hosted at the UCAR website, but the link is now dead. You can still find it on the internet though, as part of other projects (for example here, where it is being used via a Python wrapper). But the licensing is unclear to me.

The dearth of easy-to-find, easy-to-use, high-quality open source modern Fortran libraries is a big problem for the long-term future of the language. Fortran users don't seem to be as interested in promoting their language as users of some of the newer programming languages are. There is a group of us working to change this, but we've got a long way to go. For example, Julia, a programming language only four years old, already has a ton of libraries, some of which are just wrappers to Fortran 77 code that no one's ever bothered to modernize. They even have a yearly conference. An in depth article about this state of affairs is a post for another day.

References

• G. Hoffmann, Multidimensional Linear Interpolation

Numerical Differentiation

I present the initial release of a new modern Fortran library for computing Jacobian matrices using numerical differentiation. It is called NumDiff and is available on GitHub. The Jacobian is the matrix of partial derivatives of a set of \(m\) functions with respect to \(n\) variables:

Typically, each variable \(x_i\) is perturbed by a value \(h_i\) using forward, backward, or central differences to compute the Jacobian one column at a time (\(\partial \mathbf{f} / \partial x_i\)). Higher-order methods are also possible [1]. The following finite difference formulas are currently available in the library:

• Two points:
  • \((f(x+h)-f(x)) / h\)
  • \((f(x)-f(x-h)) / h\)
• Three points:
  • \((f(x+h)-f(x-h)) / (2h)\)
  • \((-3f(x)+4f(x+h)-f(x+2h)) / (2h)\)
  • \((f(x-2h)-4f(x-h)+3f(x)) / (2h)\)
• Four points:
  • \((-2f(x-h)-3f(x)+6f(x+h)-f(x+2h)) / (6h)\)
  • \((f(x-2h)-6f(x-h)+3f(x)+2f(x+h)) / (6h)\)
  • \((-11f(x)+18f(x+h)-9f(x+2h)+2f(x+3h)) / (6h)\)
  • \((-2f(x-3h)+9f(x-2h)-18f(x-h)+11f(x)) / (6h)\)
• Five points:
  • \((f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)) / (12h)\)
  • \((-3f(x-h)-10f(x)+18f(x+h)-6f(x+2h)+f(x+3h)) / (12h)\)
  • \((-f(x-3h)+6f(x-2h)-18f(x-h)+10f(x)+3f(x+h)) / (12h)\)
  • \((-25f(x)+48f(x+h)-36f(x+2h)+16f(x+3h)-3f(x+4h)) / (12h)\)
  • \((3f(x-4h)-16f(x-3h)+36f(x-2h)-48f(x-h)+25f(x)) / (12h)\)

The basic features of NumDiff are listed here:

• A variety of finite difference methods are available (and it is easy to add new ones).
• If you want, you can specify a different finite difference method to use to compute each column of the Jacobian.
• You can also specify the number of points of the methods, and a suitable one of that order will be selected on-the-fly so as not to violate the variable bounds.
• I also included an alternate method using Neville's process which computes each element of the Jacobian individually [2]. It takes a very large number of function evaluations but produces the most accurate answer.
• A hash table based caching system is implemented to cache function evaluations to avoid unnecessary function calls if they are not necessary.
• It supports very large sparse systems by compactly storing and computing the Jacobian matrix using the sparsity pattern. Optimization codes such as SNOPT and Ipopt can use this form.
• It can also return the dense (\(m \times n\)) matrix representation of the Jacobian if that sort of thing is your bag (for example, the older SLSQP requires this form).
• It can also return the \(J*v\) product, where \(J\) is the full (\(m \times n\)) Jacobian matrix and v is an (\(n \times 1\)) input vector. This is used for Krylov type algorithms.
• The sparsity pattern can be supplied by the user or computed by the library.
• The sparsity pattern can also be partitioned so as compute multiple columns of the Jacobian simultaneously so as to reduce the number of function calls [3].
• It is written in object-oriented Fortran 2008. All user interaction is through a NumDiff class.
• It is open source with a BSD-3 license.

I haven't yet really tried to fine-tune the code, so I make no claim that it is the most optimal it could be. I'm using various modern Fortran vector and matrix intrinsic routines such as PACK, COUNT, etc. Likely there is room for efficiency improvements. I'd also like to add some parallelization, either using OpenMP or Coarray Fortran. OpenMP seems to have some issues with some modern Fortran constructs so that might be tricky. I keep meaning to do something real with coarrays, so this could be my chance.

So, there you go internet. If anybody else finds it useful, let me know.

References

1. G. Engeln-Müllges, F. Uhlig, Numerical Algorithms with Fortran, Springer-Verlag Berlin Heidelberg, 1996.
2. J. Oliver, "An algorithm for numerical differentiation of a function of one real variable", Journal of Computational and Applied Mathematics 6 (2) (1980) 145–160. [A Fortran 77 implementation of this algorithm by David Kahaner was formerly available from NIST, but the link seems to be dead. My modern Fortran version is available here.]
3. T. F. Coleman, B. S. Garbow, J. J. Moré, "Algorithm 618: FORTRAN subroutines for estimating sparse Jacobian Matrices", ACM Transactions on Mathematical Software (TOMS), Volume 10 Issue 3, Sept. 1984.