SLATEC

slatecThe SLATEC Common Mathematical Library (CML) is written in FORTRAN 77 and contains over 1400 general purpose mathematical and statistical routines. SLATEC is an acronym for the “Sandia, Los Alamos, Air Force Weapons Laboratory Technical Exchange Committee”, an organization formed in 1974 by the computer centers of these organizations. In 1977, it was decided to build a FORTRAN library to provide portable, non-proprietary, mathematical software for member sites’ supercomputers. Version 1.0 of the CML was released in April 1982.

An example SLATEC routine is shown below, which computes the inverse hyperbolic cosine:

*DECK ACOSH
      FUNCTION ACOSH (X)
C***BEGIN PROLOGUE  ACOSH
C***PURPOSE  Compute the arc hyperbolic cosine.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4C
C***TYPE      SINGLE PRECISION (ACOSH-S, DACOSH-D, CACOSH-C)
C***KEYWORDS  ACOSH, ARC HYPERBOLIC COSINE, ELEMENTARY FUNCTIONS, FNLIB,
C             INVERSE HYPERBOLIC COSINE
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ACOSH(X) computes the arc hyperbolic cosine of X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  ACOSH
      SAVE ALN2,XMAX
      DATA ALN2 / 0.6931471805 5994530942E0/
      DATA XMAX /0./
C***FIRST EXECUTABLE STATEMENT  ACOSH
      IF (XMAX.EQ.0.) XMAX = 1.0/SQRT(R1MACH(3))
C
      IF (X .LT. 1.0) CALL XERMSG ('SLATEC', 'ACOSH', 'X LESS THAN 1',
     +   1, 2)
C
      IF (X.LT.XMAX) ACOSH = LOG (X + SQRT(X*X-1.0))
      IF (X.GE.XMAX) ACOSH = ALN2 + LOG(X)
C
      RETURN
      END

acoshOf course, this routine is full of nonsense from the point of view of a modern Fortran programmer (fixed-form source, implicitly typed variables, SAVE and DATA statements, constants computed at run time rather than compile time, hard-coding of the numeric value of ln(2), etc.). And you don’t even want to know what’s happening in the R1MACH() function. A modern version would look something like this:

pure elemental function acosh(x) result(acoshx)

!! arc hyperbolic cosine

use iso_fortran_env, only: wp => real64

implicit none

real(wp),intent(in) :: x
real(wp) :: acoshx

real(wp),parameter :: aln2 = &
  log(2.0_wp)
real(wp),parameter :: r1mach3 = &
  real(radix(1.0_wp),wp)**(-digits(1.0_wp)) !! largest relative spacing
real(wp),parameter :: xmax = &
  1.0_wp/sqrt(r1mach3)

if ( x<1.0_wp ) then
  error stop 'slatec : acosh : x less than 1'
else
  if ( x<xmax ) then
    acoshx = log(x + sqrt(x*x - 1.0_wp))
  else
    acoshx = aln2 + log(x)
  end if
end if

end function acosh

Of course, the point is moot, since ACOSH() is now an intrinsic function in Fortran 2008, so we don’t need either of these anymore. Unfortunately, like many great Fortran codes, SLATEC has been frozen in amber since the early 1990’s. It does contain a great many gems, however, and any deficiencies in the original Fortran 77 code can be easily remedied using modern Fortran standards. See also, for example, my own DDEABM project, which is a complete refactoring and upgrade of the Adams-Bashforth-Moulton ODE solver from SLATEC.

The GNU Scientific Library was started in 1996 to be a “modern version of SLATEC“. Of course, rather than refactoring into modern Fortran (Fortran 90 had been published several years earlier, but I believe it was some time before a free compiler was available) they decided to start over in C. There’s now even a Fortran interface to GSL (so now we can use a Fortran wrapper to call a C function that replaced a Fortran function that worked perfectly well to begin with!) The original SLATEC was a public domain work of the US government, and so can be used without restrictions in free or proprietary applications, but GSL is of course, GPL, so you’re out of luck if you are unable to accept the restrictions (I mean freedom) of that license.

Note that John Burkardt has produced a version of SLATEC translated to free-form source. The license for this modified version is unclear. Many of his other codes are LGPL (still more restrictive than the original public domain codes). It’s not a full refactoring, either. The codes still contain non-standard extensions and obsolescent and now-unnecessary features such as COMMON blocks. What we really need is a full modern refactoring of SLATEC (and other unmaintained codes like the NIST CMLIB and NSWC, etc.)

References

  1. Fong, Kirby W., Jefferson, Thomas H., Suyehiro, Tokihiko, Walton, Lee, “Guide to the SLATEC Common Mathematical Library“, July 1993.
  2. SLATEC source code at Netlib.
Posted in Programming Tagged with: , , ,
3 comments on “SLATEC
  1. I agree that FORTRAN codes should be used. I also agree that these codes should be modified for Fortran standards. However I do not think that they are to be refactored to OOP. I’ve seen many examples when refactored code lost speed and robustness.

  2. Jacob says:

    I think it depends on the situation. OOP makes a lot of sense for some things, and doesn’t make sense for others.

  3. Wadud Miah says:

    I think OOP is only really useful for very large codes, e.g. multi-physics and multi-scale codes. For simple codes, OOP has a large overhead in terms of maintenance as well as performance. If OOP is used, it has to be used carefully and not blindly just because the language provides OOP. I think the computer science community use OOP a lot and sometimes properly and sometimes not very well. In computational science, I don’t think it really provides any real advantage and Fortran 95 with modules or even Fortran 2008 submodules is usually sufficient.

Leave a Reply

Your email address will not be published. Required fields are marked *

*