Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Jun 19, 2022

D1MACH Re-Revisited

slatec

NightCafe AI generated image for "old computers"

The classic Fortran routines r1mach (for single precision real numbers), d1mach (for double precision real numbers), and i1mach (for integers) were originally written in the mid-1970s for the purpose of returning basic machine or operating system dependent constants in order to provide portability. Typically, these routines had a bunch of commented out DATA statements, and the user was required to uncomment out the ones for their specific machine. The original versions included constants for the following systems:

Over the years, various versions of these routines have been created:

In the Fortran 90 versions, the original routines were modified to use various intrinsic functions available in Fortran 90. However, some hard-coded values remained in i1mach. With the advent of more recent standards, it is now possible to provide modern and completely portable versions which are described below. Starting with the Fortran 90 versions, the following changes were made:

  • The routines are now in a module. We declare two module parameters which will be used in the routines:
integer, parameter :: sp = kind(1.0)    ! single precision kind
integer, parameter :: dp = kind(1.0d0)  ! double precision kind

Of course, putting them in a module means that they are not quite drop-in replacements for the old routines. If you want to use this module with legacy code, you'd have to update your code to add a use statement. But that's the way it goes. Nobody should be writing Fortran anymore that isn't in a module, so I'm not going to encourage that.

  • The comment strings have been updated to Ford-style. I find the old SLATEC comment block style extremely cluttered. Various cosmetic changes were also made (e.g., changing the old .GE. to >=, etc.)

  • All three routines were made pure. The error stop construct was used for the out of bounds input condition (Fortran 2008 allows error stop in pure procedures).

  • New intrinsic constants available in the Fortran 2008 standard were employed (see below).

The new routines can be found on GitHub. They are shown below without the comment blocks.

The new routines

The new d1mach routine is:

D1MACH

  pure real(dp) function d1mach(i)

    integer, intent(in) :: i

    real(dp), parameter :: x = 1.0_dp
    real(dp), parameter :: b = real(radix(x), dp)

    select case (i)
    case (1); d1mach = b**(minexponent(x) - 1) ! the smallest positive magnitude.
    case (2); d1mach = huge(x)                 ! the largest magnitude.
    case (3); d1mach = b**(-digits(x))         ! the smallest relative spacing.
    case (4); d1mach = b**(1 - digits(x))      ! the largest relative spacing.
    case (5); d1mach = log10(b)
    case default
        error stop 'Error in d1mach - i out of bounds'
    end select

  end function d1mach

This is mostly just a reformatting of the Fortran 90 routine. We make x and b parameters, and use error stop for the error condition rather than this monstrosity:

          WRITE (*, FMT = 9000)
 9000     FORMAT ('1ERROR    1 IN D1MACH - I OUT OF BOUNDS')
          STOP

R1MACH

The updated R1MACH is similar, except the real values are real(sp) (default real, or single precision):

  pure real(sp) function r1mach(i)

    integer, intent(in) :: i

    real(sp), parameter :: x = 1.0_sp
    real(sp), parameter :: b = real(radix(x), sp)

    select case (i)
    case (1); r1mach = b**(minexponent(x) - 1) ! the smallest positive magnitude.
    case (2); r1mach = huge(x)                 ! the largest magnitude.
    case (3); r1mach = b**(-digits(x))         ! the smallest relative spacing.
    case (4); r1mach = b**(1 - digits(x))      ! the largest relative spacing.
    case (5); r1mach = log10(b)
    case default
      error stop 'Error in r1mach - i out of bounds'
    end select

  end function r1mach

I1MACH

The new i1mach routine is:

  pure integer function i1mach(i)

    integer, intent(in) :: i

    real(sp), parameter :: x = 1.0_sp
    real(dp), parameter :: xx = 1.0_dp

    select case (i)
    case (1);  i1mach = input_unit
    case (2);  i1mach = output_unit
    case (3);  i1mach = 0  ! Punch unit is no longer used
    case (4);  i1mach = error_unit
    case (5);  i1mach = numeric_storage_size
    case (6);  i1mach = numeric_storage_size / character_storage_size
    case (7);  i1mach = radix(1)
    case (8);  i1mach = numeric_storage_size - 1
    case (9);  i1mach = huge(1)
    case (10); i1mach = radix(x)
    case (11); i1mach = digits(x)
    case (12); i1mach = minexponent(x)
    case (13); i1mach = maxexponent(x)
    case (14); i1mach = digits(xx)
    case (15); i1mach = minexponent(xx)
    case (16); i1mach = maxexponent(xx)
    case default
      error stop 'Error in i1mach - i out of bounds'
    end select

  end function i1mach

In this one, we make heavy use of constants now in the intrinsic iso_fortran_env module, including input_unit, output_unit, error_unit, numeric_storage_size, and character_storage_size. In the Fortran 90 code, the values for I1MACH(1:4) and I1MACH(6) were totally not portable. This new one is now (almost) entirely portable.

The only remaining parameter that cannot be expressed portably is I1MACH(3), the standard punch unit. If you look at the original i1mach.f from SLATEC, there were a range of values of this parameter used by various machines (e.g. 102 for the Cray, or 7 for the IBM 360/370). Hopefully, this is not something that anybody needs nowadays, so we leave it set to 0 in the updated routine. If the Fortran committee ever adds punch_unit to iso_fortran_env, then we can update it then, and this routine will finally be finished, and the promise of a totally portable machine constants routine will finally be fulfilled.

Another interesting thing to note are I1MACH(7) (which is RADIX(1)) and I1MACH(10) (which is RADIX(1.0)). There was no I1MACH parameter for RADIX(1.0D0). There were machine architectures where I1MACH(7) /= I1MACH(10). For example, the Data General Eclipse S/200 had I1MACH(7)=2 and I1MACH(10)=16. But, it seems unlikely that there would ever be an architecture where the radix would be different for different real kinds. So maybe we are safe?

Results

On my laptop (Apple M1) with gfortran, the three routines produce the following values:

i1mach( 1) =                5
i1mach( 2) =                6
i1mach( 3) =                0
i1mach( 4) =                0
i1mach( 5) =               32
i1mach( 6) =                4
i1mach( 7) =                2
i1mach( 8) =               31
i1mach( 9) =       2147483647
i1mach(10) =                2
i1mach(11) =               24
i1mach(12) =             -125
i1mach(13) =              128
i1mach(14) =               53
i1mach(15) =            -1021
i1mach(16) =             1024

r1mach( 1) =             1.17549435E-038    (            800000 )
r1mach( 2) =             3.40282347E+038    (          7F7FFFFF )
r1mach( 3) =             5.96046448E-008    (          33800000 )
r1mach( 4) =             1.19209290E-007    (          34000000 )
r1mach( 5) =             3.01030010E-001    (          3E9A209B )

d1mach( 1) =   2.22507385850720138E-0308    (    10000000000000 )
d1mach( 2) =   1.79769313486231571E+0308    (  7FEFFFFFFFFFFFFF )
d1mach( 3) =   1.11022302462515654E-0016    (  3CA0000000000000 )
d1mach( 4) =   2.22044604925031308E-0016    (  3CB0000000000000 )
d1mach( 5) =   3.01029995663981198E-0001    (  3FD34413509F79FF )

References

posted at 01:00  ·   ·  Fortran  SLATEC  PORT  BLAS