question about "exact" constants

Accelerated Computing HPC Compilers Legacy PGI Compilers
1

Hi All,

I have a set of subroutines that solves certain types of differential equations. The algorithm makes use of some constants that are rational numbers, and are therefore exact. These constants are set in a subroutine that is called many times:

    IMPLICIT NONE
!                                                                                                   
    INTEGER, INTENT(IN)                            :: ns
    INTEGER, INTENT(IN)                            :: nsd
    REAL(KIND=dp), INTENT(OUT)               :: a(nsd,ns)
    REAL(KIND=dp), INTENT(OUT)               :: b(ns)
.
.
.
   IF(ns == 6)THEN
      b(1)=88.0_dp/1024.0_dp
      b(2)=185.0_dp/1024.0_dp
      b(3)=239.0_dp/1024.0_dp
      a(1,1)=44.0_dp/1024.0_dp
      a(1,2)=-15.0_dp/1024.0_dp
      a(1,3)=10.0_dp/1024.0_dp
      a(1,4)=-6.0_dp/1024.0_dp
.
.
.
  END IF

I have a module which selects the KIND:

!
!  gives REAL*16
!
! integer, parameter :: dp = selected_real_kind(33,4931)
!
!  gives REAL*10
!
!  integer, parameter :: dp = selected_real_kind(18,4931)
!
!  gives double precision, REAL*8
!
  integer, parameter :: dp = selected_real_kind(15,307)
!
!   gives single precision, REAL*4
!
  integer, parameter :: sp = selected_real_kind(6,37)
!

The KIND is usually set to the equivalent of REAL8, but the equivalents to REAL10 and REAL*16 are possible for some compilers.

So the questions is: Is it better to leave the rational expressions, or should I put in the decimal forms, like

   b(3) = 0.2333984375_dp

As noted above, this routine is called many times. I realize there might be on the order of microseconds in the difference of the execution time, but it is the principle of the matter.

Also related, another routine returns the weights and abscissa for Gaussian quadrature. These can be precomputed to large number of digits. For example

IF(norder == 4)THEN
  ww(1)= 6.52145154862546142626936050778001E-01_dp
  xx(1)=-3.39981043584856264802665759103245E-01_dp
  ww(2)= 6.52145154862546142626936050778001E-01_dp
  xx(2)= 3.39981043584856264802665759103245E-01_dp
  ww(3)= 3.47854845137453857373063949221999E-01_dp
  xx(3)=-8.61136311594052575223946488892809E-01_dp
  ww(4)= 3.47854845137453857373063949221999E-01_dp
  xx(4)= 8.61136311594052575223946488892809E-01_dp
  RETURN
END IF

These coefficients are suitable up to KIND=16. I assume the compiler has enough sense to take the leading digits as needed when KIND=4 or KIND=8.

Jerry

2

It’s my understanding that at compile-time, if these are constant:

b(3)=239.0_dp/1024.0_dp

Then it should be the same as

b(3) = 0.2333984375_dp

The compiler is smart enough to compute constants at compile-time, both forms should have the same representation (and as accurate as floating point numbers can be at some arbitrary precision) at run-time.

3

Thanks for the reply. My concern was more about speed, and it appears the computer is not constantly evaluating the divisions.

Jerry