Compile time detection of extended-double or quad precision

Is there a Fortran-conforming way to detect whether a processor supports precision higher than the default double-precision at compile time without resorting to the processor? More generally, can this be done without resorting to the preprocessor?

A partial (non-working) example of what I have in mind is the following:

function sum_extra_prec(a) result(asum)
   implicit none
   integer, parameter :: dp = kind(1.0d0)
   real(dp), intent(in) :: a(:)
   real(dp) :: asum

#if __NVCOMPILER
   integer, parameter :: hp = -1
#else
#if __GFORTRAN__
   integer, parameter :: xdp = 10
#else
   integer, parameter :: xdp = -1
#endif
   integer, parameter :: qp = 16
   integer, parameter :: hp = merge(xdp,qp,xdp > 0)
#endif

   if (hp > 0) then
      asum = real(sum(real(a,hp)),dp)  ! Error with nvfortran
   else
      block:
         use ddreal
         ! ... use double-double ...
      end block
   end if

end function

This works with gfortran and the Intel Fortran compilers, however it fails with nvfortran since hp = -1 is an illegal kind selector in the expression real(a,hp). This would seem to imply I need to replace the Fortran if-else block with preprocessor fences too:

#if __NVCOMPILER
! double double
! ...
#else
asum = real(sum(real(a,hp)),dp)
#endif

Since I’m already using the preprocessor in this example to select the precision, that’s fine. But assume the precision selection were done in a precision module. What can be done in this case to avoid misusing a negative kind parameter?

In C++, I believe it would be possible to use an if constexpr:

if constexpr (quad_supported) {
  // ...
else {
  // double-double
  // ...
}
1 Like

Maybe something like (untested):

module m
use iso_fortran_env, only: real_kinds, real128

logical, parameter :: has_qp = any(real_kinds==real128)
end module

1 Like

I was testing that maybe naive solution:

program main
    use, intrinsic :: iso_fortran_env

    implicit none
    integer, parameter :: wp = max(real64, real128)
    real(wp) :: x

    print *, kind(x)
end program main

It prints 16 with GFortran in Ubuntu. I don’t know if it is OK with every compiler.

If you assume real128 is equal to quad precision, you may as well write:

logical parameter :: has_qp = real128 > 0
integer, parameter :: qp = real128

But how would you pick gfortran’s real(10) assuming you want the next supported thing after double?

Does the standard mandate that the kind specifier of real128 is an integer larger than real64?

Your proposed solution does provide a workaround to the negative kind parameter issue:

program main
    use, intrinsic :: iso_fortran_env

    implicit none

    integer, parameter :: hp = real128
    integer, parameter :: dp = kind(1.0d0)

    integer, parameter :: wp = max(dp,hp)
    real(wp) :: x

    real(dp) :: a(100), asum

    call random_number(a)

    if (wp == dp) then
        ! work precision is double, which is insufficient
        ! use custom double-double instead
        error stop "Not implemented."
    else
        ! higher precision kind is available
        asum = real(sum(real(a,wp)),dp)  ! okay to reference wp
    endif
    print *,asum

end program main

No… :face_with_raised_eyebrow:

16.10.2.25 REAL32, REAL64, and REAL128
The values of these default integer scalar named constants shall be those of the kind type parameters that specify a REAL type whose storage size expressed in bits is 32, 64, and 128 respectively. If, for any of these constants, the processor supports more than one kind of that size, it is processor dependent which kind value is provided. If the processor supports no kind of a particular size, that constant shall be equal to −2 if the processor supports kinds of a larger size and −1 otherwise.

1 Like

We use this check for stdlib

In test-drive I turned it into a compile time check with

1 Like

Thank you all for the proposed solutions, piece by piece I’ve constructed something that works (with an implicit assumption that the default double precision has less than 18 decimal digits of precision).

! has_quad.f90
program main

    implicit none

    ! Default calculation precision
    integer, parameter :: dp = kind(1.0d0)

    ! 1) Check for extended double precision
    integer, parameter :: xdp = merge(-1,selected_real_kind(18), &
        selected_real_kind(18) == selected_real_kind(33))
    logical, parameter :: has_xdp = xdp > 0

    ! 2) Check for true quadruple precision
    integer, parameter :: qp = selected_real_kind(33)
    logical, parameter :: has_qp = qp > 0

    ! 2) Assign speculative higher precision
    integer, parameter :: hp = merge(xdp,merge(qp,-1,has_qp),has_xdp)
    logical, parameter :: has_hp = hp > 0
    
    ! 3) Assign work precision; if higher precision not available, 
    !    fall back to default
    integer, parameter :: wp = merge(hp,dp,has_hp)

    real(dp) :: a(100), asum
    integer :: i

    a = [(i,i=1,100)]

    if (wp == dp) then
        ! work precision equals default, so 
        ! use custom double-double instead
        error stop "Not implemented."
    else
        ! use the higher precision
        asum = real(sum(real(a,wp)),dp)
    endif
    print *, wp, asum

end program main
$ ifx -warn all has_quad_v2.f90 
$ ./a.out
          16   5050.00000000000     
$ gfortran -Wall has_quad_v2.f90 
$ ./a.out
          10   5050.0000000000000     
$ nvfortran has_quad_v2.f90 
$ ./a.out
ERROR STOP Not implemented.

4 Likes

With a conforming processor, you can use the following standard-conforming code to fetch the “list” of floating-point precisions supported by said processor:

   use, intrinsic :: iso_fortran_env, only : .., real_kinds, ..
   ..
   integer, parameter :: precisions(*) = [( precision(real(1.0, kind=real_kinds(i))),               &
                                            integer :: i=1, size(real_kinds) )]
   ..

You can then employ the “list” in your subsequent compile-time actions.

That is, once you can get a compiler to support the above, I don’t think any of them do so at the moment.

Intel might be the first to fix their compilers in oneAPI toolkit, Intel has support requests toward the same.

2 Likes

Enhanced compile-time computing in Fortran with a few better tools toward compile-time “reflective programming” as well as CONSTEXPR functions that can support the use of certain user-defined functions in constant expressions has long been part of my vision for Fortran.

2 Likes

Interestingly, while gfortran-12 fails to compile the implicit do in precisions specification, Intel’s ifort ver. 2021.7 does compile it but the result is wrong. The following code:

program t
  use, intrinsic :: iso_fortran_env
  implicit none
  integer, parameter :: precisions(*) = &
      [ (precision(real(1.0, kind=real_kinds(i))),  &
      integer :: i=1, size(real_kinds)) ]
  
  print *, precisions
  print *, real_kinds
end program t

gives all precisions equal to 6:

           6           6           6
           4           8          16
1 Like

FYI, with Absoft af95 on Ubuntu 22.04, this is what happened.

af95 --version && af95 has_quad.f90 
Absoft 64-bit Pro Fortran 22.0.3
Absoft ERROR 587: MAIN: has_quad.f90: 10, 29 
  The initialization expression must be a constant to be used with PARAMETER assignment for object "XDP".
Absoft ERROR 113: MAIN: has_quad.f90: 34, 5 
  IMPLICIT NONE is specified in the local scope, therefore an explicit type must be specified for data object "ERROR".
Absoft ERROR 724: MAIN: has_quad.f90: 34, 11 
  Unknown statement.  Expected assignment statement but found "s" instead of "=" or "=>".

Absoft Pro Fortran 22.0: 3 Errors, 0 Warnings, 0 Other messages, 0 ANSI
f90fe reported errors.

I am not sure which standard af95 follows in this regard. Maybe the constraint that “The initialization expression must be a constant to be used with PARAMETER assignment for object XXX” has been removed in recent standards? Does anyone here know about this constraint?

The Absoft compiler has been discontinued recently, yet I suppose it is still a good compiler for the time being.

Interesting test problem: is it Fortran- conforming code to have inputs to a function with different kinds? What you’re trying to do here is equivalent to:

integer :: available_prec(*) = [precision(1.0_real32),precision(1.0_real64),precision(1.0_real128)]

the processors evidently have a problem unrolling a loop that contains different real kinds, but the point is: is that allowed by the standard?

1 Like

Well, it is not. @FortranFan’s version fills precision array with all available kinds. @FedericoPerini’s - only three (at most) defined in iso_fortran_env module (real32, real64, real128). Obvious example of non-equivalence is gfortran with its 80-bit extended precision real kind

Sorry if I was unclear: my question was not about whether real(10) was in the real_kinds array, but whether the following is legal code. It seems like the compiler is interpreting the implicit loop that you’re trying to compute as:

precision([1.0,1.d0,1.0_blabla])

which I understand is different than

[precision(1.0),precision(1.d0),precision(1.0_whatever)]

Is the compiler getting it wrong? What does the standard say about implicit loops with variables of different kind? This is what I was trying to ask.

I’m not sure what you mean. The implicit do should unroll to something similar to your second example, making it an array constructor of all integer values. The first example is IMHO wrong as you cannot have values of different types/kinds in an array constructor.

Nevermind, I think I found my answer in this old post:

https://groups.google.com/g/comp.lang.fortran/c/LY2ftkwzLbM/m/4fHDXrQDZ2sJ

My understanding per the standard is the code I posted above conforms, hence my assertion upthread.

The thread at comp.lang.fortran linked by @FedericoPerini mentions NAG Fortran compiler whose lead developer is also the Fortran standard Editor and whose compiler implementation strives to maximize conformance. So can anyone with access to latest NAG compiler try the following and post the results here? It is a somewhat simplified variant of my previous snippet, here I have removed some of the other Fortran 2008 and 2018 facilities that NAG may not support yet:

   use, intrinsic :: iso_fortran_env, only : real_kinds
   integer, parameter :: n = size(real_kinds)
   integer :: i
   integer, parameter :: precisions(n) = [( precision(real(1., real_kinds(i))), i=1,n )]
   print *, precisions
end

Going by the comp.lang.fortran thread as well as NAG documentation online, I expect the following output:

3 6 15 31

Note the type and kind of each ac-value in the implied-DO here is the result of the intrinsic function PRECISION and each of them conform as default integer, so there is no issue here.

Thank you, that’s spot on -

I also looked at the integer result at the first glance; however, the implied loop variable i affects the real(1.0,my_kind(i)) function if we assume that our loop logic should start at the innermost location. That would mean that return variables from real have different kinds, hence non-conforming.

Not trying to make a case for either one - just trying to better understand the issue

Most of the above posts are using the intrinsic array real_kinds(:) to answer the question. Here is an excerpt from some of my code from about 2004. It works with algebraic expressions because the iso_fortran_env module did not exist at that time.

   integer, parameter :: COLUMBUS_WP            = selected_real_kind(14)  
   integer, parameter :: COLUMBUS_EP_REQUESTED  = selected_real_kind(17) 
   integer, parameter :: COLUMBUS_EEP_REQUESTED = selected_real_kind(30) 

   ! the following initialization expressions are equivalent to:
   !  if ( COLUMBUS_EP_REQUESTED < 0 ) then
   !     COLUMBUS_EP = COLUMBUS_WP
   !  else
   !     COLUMBUS_EP = COLUMBUS_EP_REQUESTED
   !  endif
   integer, parameter :: COLUMBUS_EP = (1+sign(1,COLUMBUS_EP_REQUESTED))/2 * COLUMBUS_EP_REQUESTED + & 
        &                              (1-sign(1,COLUMBUS_EP_REQUESTED))/2 * COLUMBUS_WP
   integer, parameter :: COLUMBUS_EEP = (1+sign(1,COLUMBUS_EEP_REQUESTED))/2 * COLUMBUS_EEP_REQUESTED + & 
        &                               (1-sign(1,COLUMBUS_EEP_REQUESTED))/2 * COLUMBUS_EP

This works as required, but of course if new real kinds are added by some compiler, then this code would need to be extended manually. The long parameter names were used intentionally so that programmers would override them with more reasonable parameter names in their codes.

There are similar algebraic expressions using merge() that can used instead. All in all, it is nice that the KIND facility in fortran allows these types of expressions to be used to define and propagate KIND values throughout a code, but I have to think that the standards committee could have made all of this much easier on programmers than they did. In the code above, I explain what it does in the comments, but wouldn’t the whole thing have been easier if programmers had been allowed to use the if-then-else in the first place?

2 Likes