Nov 06, 2015
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.
Oct 17, 2015
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]
Aug 27, 2015
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
Aug 25, 2015
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
Aug 05, 2015
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.
Jul 17, 2015
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]
Jun 30, 2015
Good news for starving students: Intel just announced that their compilers (including the Intel Fortran compiler) are now available under a free, non-commercial license for qualified students on Linux, MacOS X and Windows.
Jun 21, 2015
SLSQP [1-2] is a sequential least squares constrained optimization algorithm written by Dieter Kraft in the 1980s. Today, it forms part of the Python pyOpt package for solving nonlinear constrained optimization problems. It was written in FORTRAN 77, and is filled with such ancient horrors like arithmetic IF statements, GOTO statements (even computed and assigned GOTO statements), statement functions, etc. It also assumes the implicit saving of variables in the LINMIN subroutine (which is totally non-standard, but was done by some early compilers), and in some instances assumes that variables are initialized with a zero value (also non-standard, but also done by some compilers). It has an old-fashioned "reverse communication" style interface, with generous use of "WORK" arrays passed around using assumed-size dummy arguments. The pyOpt users don't seem to mind all this, since they are calling it from a Python wrapper, and I guess it works fine (for now).
I don't know if anyone has ever tried to refactor this into modern Fortran with a good object-oriented interface. It seems like it would be a good idea. I've been toying with it for a while. It is always an interesting process to refactor old Fortran code. The fixed-form to free-form conversion can usually be done fairly easily with an automated script, but untangling the spaghetti can be a bit harder. Consider the following code snipet from the subroutine HFTI:
C DETERMINE PSEUDORANK
DO 90 j=1,ldiag
90 IF(ABS(a(j,j)).LE.tau) GOTO 100
k=ldiag
GOTO 110
100 k=j-1
110 kp1=k+1
A modern version that performs the exact same calculation is:
!determine pseudorank:
do j=1,ldiag
if (abs(a(j,j))<=tau) exit
end do
k=j-1
kp1=j
As you can see, no line numbers or GOTO statements are required. The new version even eliminates one assignment and one addition operation. Note that the refactored version depends on the fact that the index of the DO loop, if it completes the loop (i.e., if the exit statement is not triggered), will have a value of ldiag+1. This behavior has been part of the Fortran standard for some time (perhaps it was not when this routine was written?).
There is also a C translation of SLSQP available, which fixes some of the issues in the original code. Their version of the above code snippet is a pretty straightforward translation of the FORTRAN 77 version, except with all the ugliness of C on full display (a[j + j * a_dim1], seriously?):
/* DETERMINE PSEUDORANK */
i__2 = ldiag;
for (j = 1; j <= i__2; ++j) {
/* L90: */
if ((d__1 = a[j + j * a_dim1], fabs(d__1)) <= *tau) {
goto L100;
}
}
k = ldiag;
goto L110;
L100:
k = j - 1;
L110:
kp1 = k + 1;
See Also
- Dieter Kraft, "A Software Package for Sequential Quadratic Programming", DFVLR-FB 88-28, 1988.
- Dieter Kraft, "Algorithm 733: TOMP–Fortran Modules for Optimal Control Calculations," ACM Transactions on Mathematical Software, Vol. 20, No. 3, p. 262-281 (1994).
- John Burkardt, LAWSON Least Squares Routines, 21 October 2008. [Fortran 90 versions of the least squares routines used in SLSQP].
May 17, 2015
Topological sorting can be used to determine the order in which a collection of interdependent tasks must be performed. For example, when building a complex modern Fortran application, there can be many modules with complex interdependencies (via use association). A Fortran building system (like FoBiS, Foray, or SCons) will perform a topological sort to determine the correct compilation order. In graph theory, the collection of tasks are vertices of a directed graph. If there are no circular dependencies, then it is a directed acyclic graph (DAG).
An example modern Fortran implementation of a topological sorting algorithm is given here. The only public class is dag, which is used to define the graph and perform the sort. The toposort method performs a "depth-first" traversal of the graph using the recursive subroutine dfs. Each vertex is only visited once, and so the algorithm runs in linear time. The routine also triggers an error message for circular dependencies.
module toposort_module
!! Topological sorting, using a recursive
!! depth-first method. The vertices are
!! integer values from (1, 2, ..., nvertices)
implicit none
private
type :: vertex
!! a vertex of a directed acyclic graph (DAG)
integer,dimension(:),allocatable :: edges
integer :: ivertex = 0 !vertex number
logical :: checking = .false.
logical :: marked = .false.
contains
procedure :: set_edges
end type vertex
type,public :: dag
!! a directed acyclic graph (DAG)
type(vertex),dimension(:),allocatable :: vertices
contains
procedure :: set_vertices => dag_set_vertices
procedure :: set_edges => dag_set_edges
procedure :: toposort => dag_toposort
end type dag
contains
subroutine set_edges(me,edges)
!! specify the edge indices for this vertex
implicit none
class(vertex),intent(inout) :: me
integer,dimension(:),intent(in) :: edges
allocate(me%edges(size(edges)))
me%edges = edges
end subroutine set_edges
subroutine dag_set_vertices(me,nvertices)
!! set the number of vertices in the dag
implicit none
class(dag),intent(inout) :: me
integer,intent(in) :: nvertices !! number of vertices
integer :: i
allocate(me%vertices(nvertices))
me%vertices%ivertex = [(i,i=1,nvertices)]
end subroutine dag_set_vertices
subroutine dag_set_edges(me,ivertex,edges)
!! set the edges for a vertex in a dag
implicit none
class(dag),intent(inout) :: me
integer,intent(in) :: ivertex !! vertex number
integer,dimension(:),intent(in) :: edges
call me%vertices(ivertex)%set_edges(edges)
end subroutine dag_set_edges
subroutine dag_toposort(me,order)
!! main toposort routine
implicit none
class(dag),intent(inout) :: me
integer,dimension(:),allocatable,intent(out) :: order
integer :: i,n,iorder
n = size(me%vertices)
allocate(order(n))
iorder = 0 ! index in order array
do i=1,n
if (.not. me%vertices(i)%marked) &
call dfs(me%vertices(i))
end do
contains
recursive subroutine dfs(v)
!! depth-first graph traversal
type(vertex),intent(inout) :: v
integer :: j
if (v%checking) then
error stop 'Error: circular dependency.'
else
if (.not. v%marked) then
v%checking = .true.
if (allocated(v%edges)) then
do j=1,size(v%edges)
call dfs(me%vertices(v%edges(j)))
end do
end if
v%checking = .false.
v%marked = .true.
iorder = iorder + 1
order(iorder) = v%ivertex
end if
end if
end subroutine dfs
end subroutine dag_toposort
end module toposort_module
An example use of this module is given below. Here, we define a graph with five vertices: task 2 depends on 1, task 3 depends on both 1 and 5, and task 4 depends on 5.
program main
use toposort_module
implicit none
type(dag) :: d
integer,dimension(:),allocatable :: order
call d%set_vertices(5)
call d%set_edges(2,[1]) ! 2 depends on 1
call d%set_edges(3,[5,1]) ! 3 depends on 5 and 1
call d%set_edges(4,[5]) ! 4 depends on 5
call d%toposort(order)
write(*,*) order
end program main
The result is:
Which is the evaluation order. As far as I can find, the above code is the only modern Fortran topological sorting implementation available on the internet. There is a FORTRAN 77 subroutine here, and a Fortran 90-ish one here (however, neither of them check for circular dependencies).
See also
- Topological Sorting [Everything Under The Sun], June 26, 2013.
- Topological sorting [Wikipedia]
- tsort -- UNIX command for performing topological sorting. [Note that this gives the result in the reverse order from the code listed above.]
May 16, 2015
Part of the confusion about modern Fortran is what file extension to use. This confusion has its roots in the Fortran 90 standard (ISO/IEC 1539:1991) and the introduction of "free-form" source (the older FORTRAN 77 punched-card style was renamed "fixed-form" at this time). Unfortunately, compiler vendors began to use the extension .f90 for free-form files (common fixed-form extensions were .f or .for). Of course, the Fortran standard says nothing about file extensions. This is just one of those things the Lords of Fortran pretend are not important, like the ability to access the command line (not standardized until Fortran 2008) or the file system (still waiting for this). So, that left the compiler vendors (and users) free to do what they wanted, and all of them used .f90 for the new file type.
This was fine until the Fortran 95 standard came along (ISO/IEC 1539-1:1997). Fortran 95 was a minor update to the language (mostly cleaning up some of the half-baked bits of Fortran 90, but also adding some new stuff). This is when things started to go wrong. Since Fortran 90 used the .f90 extension, Fortran 95 should use .f95, right? And then Fortran 2003 came along, and then Fortran 2008, etc. The end result is a confusion of file extensions for free-form source such as .f90, .f95, .f03, .f08, and soon .f15. However, it doesn't make any sense to use the file extension to indicate the Fortran standard. Fortran has a high degree of backward compatibility (Python users take note). So, for example, almost all of Fortran 95 code is valid Fortran 90 code. Even a good deal of FORTRAN 77 code is still perfectly valid Fortran 2008 code. In fact, fixed-form source is still valid even in the Fortran 2008 standard (although only a lunatic would write object-oriented code in fixed-form format).
To sum up: .f90 means free-form source, not Fortran 90, and it's the only extension you should to be using. Stop the madness, or we're really going to be in trouble 75 years from now.
See also
- Steve Lionel, Source Form Just Wants to be Free [Intel.com] January 11, 2013.
