Some years ago, I needed a Fortran routine to do linear interpolation of a multidimensional (up to 6D) data set. Not wanting to reinvent the wheel, I ended up using a routine called FINT from CERNLIB (documentation here). It is written in Fortran 66, and has some hard-coded limits on the number of dimensions that can be used, but these are not fundamental to the algorithm. With a few minor changes, the routine can be used for interpolation of any number of dimensions. I present below a modern Fortran version that removes the limitations. Note that the original code is GPL licensed, thus so is this modification.
function fint(narg,arg,nent,ent,table)
use iso_fortran_env, only: wp => real64
implicit none
real(wp) :: fint
integer,intent(in) :: narg
integer,dimension(narg),intent(in) :: nent
real(wp),dimension(narg),intent(in) :: arg
real(wp),dimension(:),intent(in) :: ent
real(wp),dimension(:),intent(in) :: table
integer,dimension(2**narg) :: index
real(wp),dimension(2**narg) :: weight
logical :: mflag,rflag
real(wp) :: eta,h,x
integer :: i,ishift,istep,k,knots,lgfile,&
lmax,lmin,loca,locb,locc,ndim
!some error checks:
if (size(ent)/=sum(nent)) &
error stop 'size of ent is incorrect.'
if (size(table)/=product(nent)) &
error stop 'size of table is incorrect.'
if (narg<1) &
error stop 'invalid value of narg.'
lmax = 0
istep = 1
knots = 1
index(1) = 1
weight(1) = 1.0_wp
main: do i = 1 , narg
x = arg(i)
ndim = nent(i)
loca = lmax
lmin = lmax + 1
lmax = lmax + ndim
if ( ndim>2 ) then
locb = lmax + 1
do
locc = (loca+locb)/2
if ( x<ent(locc) ) then
locb = locc
else if ( x==ent(locc) ) then
ishift = (locc-lmin)*istep
do k = 1 , knots
index(k) = index(k) + ishift
end do
istep = istep*ndim
cycle main
else
loca = locc
end if
if ( locb-loca<=1 ) exit
end do
loca = min(max(loca,lmin),lmax-1)
ishift = (loca-lmin)*istep
eta = (x-ent(loca))/(ent(loca+1)-ent(loca))
else
if ( ndim==1 ) cycle main
h = x - ent(lmin)
if ( h==0.0_wp ) then
istep = istep*ndim
cycle main
end if
ishift = istep
if ( x-ent(lmin+1)==0.0_wp ) then
do k = 1 , knots
index(k) = index(k) + ishift
end do
istep = istep*ndim
cycle main
end if
ishift = 0
eta = h/(ent(lmin+1)-ent(lmin))
end if
do k = 1 , knots
index(k) = index(k) + ishift
index(k+knots) = index(k) + istep
weight(k+knots) = weight(k)*eta
weight(k) = weight(k) - weight(k+knots)
end do
knots = 2*knots
istep = istep*ndim
end do main
fint = 0.0_wp
do k = 1 , knots
i = index(k)
fint = fint + weight(k)*table(i)
end do
end function fint
As I recall, I wrote wrappers to this function in order to be able to input the data in a more reasonable format. For example, for a five dimensional data set it is more natural to input the five rank-1 abscissa vectors and the rank-5 ordinate matrix. The disadvantage of this routine is that you have to stuff the ordinate matrix into a vector (as well as stuffing all the abscissa vectors into a single vector). This could also be an advantage in some cases, since Fortran arrays have a maximum rank. For Fortran 90, it was 7, while Fortran 2008 expanded this to 15 (the current Intel Fortran Compiler allows up to 31 as a non-standard extension). If you ever, for example, needed to interpolate 32-dimensional data (assuming you had enough RAM), you could do it with FINT without having to declare any rank-32 matrices. In any event, an example 5D wrapper for FINT is given here:
function fint_5d_wrapper(x,y,z,q,r,&
xvec,yvec,zvec,qvec,&
rvec,fmat) result(f)
use iso_fortran_env, only: wp => real64
implicit none
integer,parameter :: narg = 5
real(wp) :: f
real(wp),intent(in) :: x,y,z,q,r
real(wp),dimension(:),intent(in) :: xvec
real(wp),dimension(:),intent(in) :: yvec
real(wp),dimension(:),intent(in) :: zvec
real(wp),dimension(:),intent(in) :: qvec
real(wp),dimension(:),intent(in) :: rvec
real(wp),dimension(:,:,:,:,:),intent(in) :: fmat
integer,dimension(narg) :: nent
nent(1) = size(xvec)
nent(2) = size(yvec)
nent(3) = size(zvec)
nent(4) = size(qvec)
nent(5) = size(rvec)
f = fint(narg,[x,y,z,q,r],nent,&
[xvec,yvec,zvec,qvec,rvec],&
reshape(fmat,[product(nent)]))
end function fint_5d_wrapper
This works well enough, but note that it isn't very computationally efficient since we are creating temporary copies of the potentially large data arrays (including having to reshape a rank-5 matrix) every time we do an interpolation. While the logic of FINT is a bit hard to fathom, the actual equations are fairly simple (see here, for example). What we really need now is an object-oriented modern Fortran library for multidimensional linear interpolation. Stay tuned...
References
DDEABM is variable step size, variable order Adams-Bashforth-Moulton PECE solver for integrating a system of first order ordinary differential equations [1-2]. It is a public-domain code originally developed in the 1970s-1980s, written in FORTRAN 77, and is available from Netlib (as part of the SLATEC Common Mathematical Library). DDEABM is primarily designed to solve non-stiff and mildly-stiff differential equations when derivative evaluations are expensive, high accuracy results are needed or answers at many specific points are required. It is based on the earlier ODE/STEP/INTRP codes described in [3-4].
DDEABM is a great code, but like many of the greats of Fortran, it seems to have been frozen in amber for 30 years. There are a couple of translations into other programming languages out there such as IDL and Matlab. It is also the ancestor of the Matlab ode113 solver (indeed, they were both written by the same person [5]). But, it looks like poor Fortran users have been satisfied with using it as is, in all its FORTRAN 77 glory.
So, I've taken the code and significantly refactored it to bring it up to date to modern standards (Fortran 2003/2008). This is more than just a conversion from fixed to free-form source. The updated version is now object-oriented and thread-safe, and also has a new event finding capability (there is a version of this code that had root-finding capability, but it seems to be based on an earlier version of the code, and also has some limitations such as a specified maximum number of equations). The new event finding feature incorporates the well-known ZEROIN algorithm [6-7] for finding a root on a bracketed interval. Everything is wrapped up in an easy-to-use class, and it also supports the exporting of intermediate integration points.
The new code is available on GitHub and is released under a permissive BSD-style license. It is hoped that it will be useful. There are some other great ODE codes that could use the same treatment (e.g. DLSODE/DVODE from ODEPACK, DIVA from MATH77, and DOP853 from Ernst Hairer).
References
- L. F. Shampine, H. A. Watts, "DEPAC - Design of a user oriented package of ode solvers", Report SAND79-2374, Sandia Laboratories, 1979.
- H. A. Watts, "A smoother interpolant for DE/STEP, INTRP and DEABM: II", Report SAND84-0293, Sandia Laboratories, 1984.
- L. F. Shampine, M. K. Gordon, "Solving ordinary differential equations with ODE, STEP, and INTRP", Report SLA-73-1060, Sandia Laboratories, 1973.
- L. F. Shampine, M. K. Gordon, "Computer solution of ordinary differential equations, the initial value problem", W. H. Freeman and Company, 1975.
- L. F. Shampine and M. W. Reichelt, "The MATLAB ODE Suite" [MathWorks].
- R. P. Brent, "An algorithm with guaranteed convergence for finding a zero of a function", The Computer Journal, Vol 14, No. 4., 1971.
- R. P. Brent, "Algorithms for minimization without derivatives", Prentice-Hall, Inc., 1973.
Good news everyone! The US government just announced that it has reached an agreement with NVIDIA to produce an open source Fortran front-end for the LLVM compiler infrastructure. It will be based on the existing commercial Portland Group compiler (NVIDIA purchased the Portland Group a couple of years ago). Source code for the Fortran front-end is expected to be available in late 2016. From the announcement:
The project is being spearheaded by the Lawrence Livermore, Sandia and Los Alamos national laboratories in response to the need for a robust open-source Fortran solution to complement and support the burgeoning use of LLVM and the CLANG C++ compiler in the HPC community. Large HPC applications, such as those developed by the NNSA Laboratories, are often built on mixed-language modules, and require a common compiler infrastructure that supports both C/C++ and Fortran. Fortran also remains widely used in the broader scientific computing community, supporting simulation science to advance national security, medicine, energy, climate and basic science missions.
Does Fortran support for LLVM mean that we'll eventually be able to have Fortran code running on our iPads? Only time will tell...
See also
My bspline-fortran multidimensional interpolation library is now at version 4.0. The documentation can be found here. Since I first mentioned it here, I've made many updates to this library, including:
- Added object-oriented wrappers to the core routines. The user now has the choice to use the older subroutine interface or the new object-oriented interface.
- Added a set of 1D routines (suggested by a user on GitHub). The library now works for 1D-6D data sets.
- Everything is now thread-safe for your multithreaded pleasure.
The new object-oriented interface makes it pretty easy to use from modern Fortran. The classes have only three methods (initialize, evaluate, and destroy). For example, a 3D case would look like this:
type(bspline_3d) :: s
call s%initialize(x,y,z,fcn,kx,ky,kz,iflag)
call s%evaluate(xval,yval,zval,idx,idy,idz,f,iflag)
call s%destroy()
The core routines of bspline-fortran were written in the early 1980s. Good Fortran code can live on for decades after it is written. While this is great, it also means that there is a lot of old Fortran code out there that is written in a coding style that no modern programmer should accept. This can be OK for well documented library routines that the user never needs to change (see, for example SPICE or LAPACK). However, refactoring old code using modern concepts can provide many advantages, as is demonstrated in this case.
See also
- DBSPLIN and DTENSBS from the NIST Core Math Library.
- Carl de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978.
Fortran's TRANSFER function (introduced in Fortran 95) is a somewhat strange function that can be used to do some interesting things not otherwise possible in earlier versions of the language. In fact, it can be used to achieve a sort of poor-man's polymorphism. It can also still be useful in modern Fortran.
All this function does is transfers the bitwise representation of one variable into another variable. This is different from casting a variable of one type to another. For example:
program transfer_example
use iso_fortran_env
integer(int64) :: i
real(real64) :: d
i = 1
d = i !d is 1.0000000000000000
d = transfer(i,d) !d is 4.9406564584124654E-324
end program transfer_example
The following example shows how to use TRANSFER to create a hash function that can operate on any variable or derived type. The dummy argument of the hash_anything function is an unlimited polymorphic variable (class(*)). The bitwise representation of this variable is then transferred into a character string big enough to hold it (making use of the Fortran 2008 STORAGE_SIZE intrinsic), which can then be hashed using any hash function that takes a character string as an input (a basic DJB method is used here).
module hash_module
!integer kind to use for the hash:
use iso_fortran_env, only: ip => INT32
implicit none
contains
function hash_anything(v) result(hash)
!! Hash anything, using unlimited
!! polymorphic variable and transfer()
class(*),intent(in) :: v
integer(ip) :: hash
!!a default character string
!!with the same size as v:
character(len=ceiling(dble(storage_size(v))/&
dble(storage_size(' ')))) :: str
str = transfer(v,str) !transfer bits of v into
!string of the same size
hash = djb_hash(str) !hash it
end function hash_anything
function djb_hash(str) result(hash)
!! Character string hashing with
!! the DJB algorithm. From [1].
character(len=*),intent(in) :: str
integer(ip) :: hash
integer :: i
hash = 5381_ip
do i=1,len(str)
hash = (ishft(hash,5_ip) + hash) + &
ichar(str(i:i))
end do
end function djb_hash
end module hash_module
See also
- Fortran hashing algorithm, July 6, 2013 [Fortran Dev]
- Hash table example [Fortran Wiki]
Hey, there's an article on the internet that mentions Coarray Fortran (introduced in the Fortran 2008 standard), and how it is the right tool for a particular massive computational job. It doesn't mention punchcards, the 1950s, or express any surprise that people are still using Fortran in modern times. Progress?
Excerpt:
Applications such as ECMWF's IFS model use OpenMP (a standard for shared memory parallelization) for computations on the nodes and MPI for communications across the nodes. Generally, this arrangement works well for many applications, but for high-resolution weather modeling using the spectral transform method, it can bog down the application, wasting precious runtime—it simply does not easily allow for any overlap of computation and communication. The use of coarrays (also called PGAS, for Partitioned Global Address Space) in Fortran, however, allows a simple syntax for data to be communicated between tasks at the same time computations are being performed.
See also
Intel just announced the availability of version 16.0 of the Intel Fortran Compiler (part of Intel Parallel Studio XE 2016). New features include:
- Submodules (Fortran 2008)
IMPURE ELEMENTAL (Fortran 2008)EXIT from BLOCK (Fortran 2008)- Full "Further Interoperability with C" implementation from TS29113 (Fortran 2015)
- Improved coarray performance on Linux and Windows
- On Windows, support for Visual Studio 2015
The new "C Descriptor" feature from Fortran 2015 looks pretty awesome, and finally plugs some of the holes in the C interoperability feature first introduced in Fortran 2003. It now allows C code to interact with Fortran pointers, allocatable variables, assumed shape variables, and CHARACTER(*) strings.
I used to be excited about submodules, which have the potential to reduce compilation cascades for very large projects. However, since Intel added the multi-processor compile option on Windows a few versions ago, that has been less bothersome for me. I'll probably try it out anyway, since it also might be useful for organizational purposes (splitting up very large modules into multiple submodules).
See also
- Intel Fortran Compiler 16.0 Release Notes
- What’s New in Intel Fortran 16.0
The 4.2.0 release of the JSON-Fortran library is now available on GitHub. This version has a few new features and a couple of minor bug fixes. The source code documentation is also now produced by FORD, which is a great new tool for modern Fortran documentation. It has been out for less than a year and is already far better than other tools that have been around for years. See an example of the output here.
It looks like JPL just opensourced their MATH77 Fortran library. MATH77 is a library of Fortran 77 subroutines implementing various numerical algorithms. It was developed over decades at JPL, and contains some very high-quality and time-tested code. The code is released under a BSD-type license. There is also a C version for people who love semicolons.
See also
- "MATH77/mathc90 Libraries Available at netlib", July 17, 2015. [comp.lang.fortran]