FTN95 support for REAL KIND=4 constants and computations.

I'm trying to compile an open source FORTRAN 90 program which uses REAL*16 precision.

FTN95 appears to recognize this (KIND=4) since:

for the command line real*16 arg1,c,darg,x1 there is no error/warning.

for maxint=lnum+3*ndec+iqint(c)+5 there is no warning.

and for a later IQINT(1.23E7) i get the message warning 189 - In a previous call to IQINT, the first argument was of type REAL(KIND=4), it is now REAL(KIND=1).

BUT, the 'q' format for KIND 4 constants is not recognized:

for api=pi/180.q0 i get the message error 773 - Variable Q follows another operand (possible unexpected space?)

and if i try to change the constant to 180.0_4 i get the message error 636 - KIND parameter out of range, permitted KINDs are 1, 2, or 3.

So, my questions are:

1. Is there any way to make the FTN95 support REAL KIND=4?

2. If not, can you recommend a FORTRAN 90 compiler (i am ready to pay for it) that does support REAL KIND=4. (For usage on a Windows 7 system).

Addition: I installed the gFortran and it recognizes REAL KIND=4 constants. However, when building the project i get a link error message

'undefined reference to `qcos''. (as well as to many other intrinsic functions such as qlog10 etc.).

so now my question is:

3. How to include the intrinsic 'q'uadruple precision functions in the build of the gFortran program (in Plato).

Thank you.

FTN95 does not support REAL16. However it does provide support for REAL10, which is an extended precision. While FTN95 might allow the syntax of REAL*16 (or possibly KIND=4) without an error message, as far as I know there is no support for this higher precision.

You might like to try the following program that shows what precision is supported. ! kind test for FTN95 module kind_precision integer, parameter :: int_1byte = selected_int_kind (2) ! 1 integer, parameter :: int_2byte = selected_int_kind (4) ! 2 integer, parameter :: int_4byte = selected_int_kind (9) ! 3 integer, parameter :: int_8byte = selected_int_kind (18) ! 4 integer, parameter :: real_4byte = selected_real_kind (6,37) ! 1 integer, parameter :: real_8byte = selected_real_kind (15,307) ! 2 integer, parameter :: real_long = selected_real_kind (18,4931) ! 3 end module kind_precision

use kind_precision 
! 
   integer ( kind=int_4byte ) :: i,k, last_k 
   character kind_name*10 
   real ( kind=real_8byte ) :: x 
   real ( kind=real_long )   :: y 
! 
!  confirm integer precision 
      write (*,*) ' ' 
      write (*,*) 'Test of integer KIND' 
      last_k = -99 
      do i = 50,0,-1 
         k = selected_int_kind (i) 
         if (k == last_k) cycle 
         kind_name = ' undefined' 
         if ( k == int_1byte ) kind_name = '1 byte' 
         if ( k == int_2byte ) kind_name = '2 byte' 
         if ( k == int_4byte ) kind_name = '4 byte' 
         if ( k == int_8byte ) kind_name = '8 byte' 
         write (*,*) kind_name, ' precision =', i, ' kind =',k 
         last_k = k 
      end do 
! 
!  confirm real precision 
      write (*,*) ' ' 
      write (*,*) 'Test of Real precision KIND' 
      last_k = -99 
      do i = 50,0,-1 
         k = selected_real_kind (i,1) 
         if (k == last_k) cycle 
         kind_name = ' undefined' 
         if ( k == real_4byte  ) kind_name = '4 byte' 
         if ( k == real_8byte  ) kind_name = '8 byte' 
         if ( k == real_long   ) kind_name = '10 byte' 
         write (*,*) kind_name, ' precision =', i, ' kind =',k 
         last_k = k 
      end do 
! 
!  confirm real exponent 
      write (*,*) ' ' 
      write (*,*) 'Test of Real exponent KIND' 
      last_k = -99 
      do i = 5000,0,-1 
         k = selected_real_kind (1,i) 
         if (k == last_k) cycle 
         kind_name = ' undefined' 
         if ( k == real_4byte  ) kind_name = '4 byte' 
         if ( k == real_8byte  ) kind_name = '8 byte' 
         if ( k == real_long   ) kind_name = '10 byte' 
         write (*,*) kind_name, ' exponent =', i, ' kind =',k 
         last_k = k 
      end do 
! 
!   test real constants 
! 
      x = 1.0e300_real_8byte    ; write (*,*) 'x=',x 
      y = 1.0e3000_real_long    ; write (*,*) 'y=',y 
      x = 1.0d300               ; write (*,*) 'x=',x 
!      x = 1.0e300               ; write (*,*) 'x=',x 
!      y = 1.0e3000              ; write (*,*) 'y=',y 
! 
      end 

or real16 x real10 y real8 z real4 a ! write (,) 'kind x16 =', kind(x) write (,) 'kind y10 =', kind(y) write (,) 'kind z8 =', kind(z) write (,) 'kind a4 =', kind(a) end

I have looked further and it appears that FTN95 does offer some compatibility support for REAL*16, but not the precision. I modified the second program above and ran it with FTN95 and gFortran 4.9.1

!  Program run on gFortran includes sizeof intrinsic
      real*16 x(2), pi_16, one_16
      real*10 y(2), pi_10, one_10
      real*8 z(2),  pi_8,  one_8
      real*4 a(2),  pi_4,  one_4
!
      write (*,*) 'kind x*16 =', kind(x), precision(x), loc(x(2))-loc(x(1)), sizeof(x(1))
      write (*,*) 'kind y*10 =', kind(y), precision(y), loc(y(2))-loc(y(1)), sizeof(y(1))
      write (*,*) 'kind z*8  =', kind(z), precision(z), loc(z(2))-loc(z(1)), sizeof(z(1))
      write (*,*) 'kind a*4  =', kind(a), precision(a), loc(a(2))-loc(a(1)), sizeof(a(1))
!
      one_16 = 1   
      pi_16  = 4 * atan (one_16)
      write (*,*) pi_16
      write (*,'(f30.25)') pi_16
!
      one_10 = 1   
      pi_10  = 4 * atan (one_10)
      write (*,*) pi_10
      write (*,'(f30.25)') pi_10
!
      one_8  = 1   
      pi_8   = 4 * atan (one_8 )
      write (*,*) pi_8 
      write (*,'(f30.25)') pi_8 
!
      one_4  = 1   
      pi_4   = 4 * atan (one_4 )
      write (*,*) pi_4 
      write (*,'(f30.25)') pi_4 
   end

Ftn95 output
[FTN95/Win32 Ver. 7.10.0 Copyright (c) Silverfrost Ltd 1993-2014]
    NO ERRORS, 4 WARNINGS  [<main program> FTN95/Win32 v7.10.0]
Creating executable: c:\\temp\\lgotemp@.exe
Program entered
 kind x*16 =           3          18          16
 kind y*10 =           3          18          10
 kind z*8  =           2          15           8
 kind a*4  =           1           6           4
      3.1415926535897932    
   3.1415926535897932386000000
      3.1415926535897932    
   3.1415926535897932386000000
           3.14159265359    
   3.1415926535897931160000000
      3.14159    
   3.1415927410125732423000000

gFortran output
gcc version 4.9.1 (x86_64-posix-seh-rev0, Built by MinGW-W64 project) 
 kind x*16 =          16          33                   16                   16
 kind y*10 =          10          18                   16                   16
 kind z*8  =           8          15                    8                    8
 kind a*4  =           4           6                    4                    4
   3.14159265358979323846264338327950280      
   3.1415926535897932384626434
   3.14159265358979323851      
   3.1415926535897932385128090
   3.1415926535897931     
   3.1415926535897931159979635
   3.14159274    
   3.1415927410125732421875000

The FTN95 output shows that it allocates the correct number of bytes for real16 but only provides real10 precision The gFortran indicates it provides 16 bytes for real*10 ? but does provide the correct precision.

It appears that with FTN95, you can test the memory allocation of REAL16 but not get the appropriate precision. As is required, both REAL10 and REAL*16 have the same exponent range.

It would be a good addition if SIZEOF could be supported or BIT_SIZE could be extended to Real and complex data types.

John

John's summary appears to be correct. It looks like FTN95 should not accept real16 and maybe it is being confused with complex16 which is OK.

I will log this for investigation.

As far as I know, FTN95 is one of the few compilers that accept real*10.

The real16 does appear to act like a real10, but with 16 bytes of memory. This is not what I would expect of complex.

I think there is a significant performance penalty with Real*10. I would guess it is emulated as I don't think the 8087(?) maths coprocessor capability is available on recent processors.

One of the advantages of the apparent support of real16 is that it allows programs with real16 to run with the expected exponent range, although the precision is not what is expected.

John

John and Paul, thank you very much for your detailed and quick response. you are fantastic.

John, from your demo gFortran program it looks like the atan function is equivalent to the qtan function. Is this correct? does it support the REAL*16 precision?

Also, can you help me with my third question - how can i tell the Plato linker to access the 'libqmath...' libraries when linking a gFortran compiled program (in case the 'axxx' functions do not give the required precision?

Thanks again.

zeevBP,

I have not checked the values returned by gF real*16, but I would expect they achieve the correct precision.

I have not yet attempted to compile or link gFortran from in Plato. I'd like to know myself. I use batch files to define the appropriate path for gFortran.exe then other batch files for compile and load. With 4.9.1 I think that the library paths are defined internally relative to the bin\.. path.

John

ps: Here are comparisons of FTN95 (Ver 7.1) , gFortran (Ver 4.9.1) and LaheyFortran (Ver 5.55j). All appear to give reasonable precision for the KIND claimed. FTN95 provides fewer significant figures when writing with write (,) although the F30.25 indicates more precision might be available.

3.1415926535897932384626433832795028   *16(LF95)
3.14159265358979323846264338327950280  *16(gFortran)
3.14159265358979323851                 *10(gFortran)
3.1415926535897932386000000
3.1415926535897932                     *10(FTN95)
3.141592653589793                      *8 (LF95)
3.1415926535897931                     *8 (gFortran)
3.1415926535897931160000000
3.14159265359                          *8 (FTN95)
3.14159274                             *4 (LF95)
3.14159274                             *4 (gFortran)
3.1415927410125732423000000
3.14159                                *4 (FTN95)

John, thank you for your support.

We (my son) found out in the documentation of gFortran that although it accepts the quadruple precision functions (such as qcos) as a valid part of the program syntax, there is currently no implementation of these functions - and that is why i get the error messages in the link phase.

Since i think i really need the REAL*16 precision, i am now looking at some commercial compilers which claim support for these function (specifically a Fortran compiler running under Intel Parallel Studio XE 2015).

Thanks again - if it turns out that REAL*10 support is good enough - i am sure i will come back to FTN95!

FTN95 REAL*10 gives 18 significant figures and exponents up to + or - 4932.

zeevBP,

I am not sure what your son was reading, but why do you need QCOS. Simply use the generic COS and provide a real*16 argument. The attached program uses ATAN and ACOS and provides the precision required.

I am also unsure of your need for 128 bit reals. I have always struggled to identify then confirm this precision. You are basically saying you can provide data inputs to the calculation to better than 18 significant figures, which is a big ask. I would be looking at the numerical algorithm to see if there is a better way that doesn't loose so much precision and can work with real*8, or *10 for larger exponents.

John

      real*16 x(2), pi_16, one_16, x_16
      real*10 y(2), pi_10, one_10
      real*8 z(2),  pi_8,  one_8
      real*4 a(2),  pi_4,  one_4
!
      write (*,*) 'kind x*16 =', kind(x), precision(x), loc(x(2))-loc(x(1)), sizeof(x(1))
      write (*,*) 'kind y*10 =', kind(y), precision(y), loc(y(2))-loc(y(1)), sizeof(y(1))
      write (*,*) 'kind z*8  =', kind(z), precision(z), loc(z(2))-loc(z(1)), sizeof(z(1))
      write (*,*) 'kind a*4  =', kind(a), precision(a), loc(a(2))-loc(a(1)), sizeof(a(1))
!
      one_16 = 1   
      pi_16  = 4 * atan (one_16)
      write (*,*) pi_16
      write (*,'(f30.25)') pi_16
!
      one_10 = 1   
      pi_10  = 4 * atan (one_10)
      write (*,*) pi_10
      write (*,'(f30.25)') pi_10
!
      one_8  = 1   
      pi_8   = 4 * atan (one_8 )
      write (*,*) pi_8 
      write (*,'(f30.25)') pi_8 
!
      one_4  = 1   
      pi_4   = 4 * atan (one_4 )
      write (*,*) pi_4 
      write (*,'(f30.25)') pi_4 
!
      x_16 = 2 * acos (one_16)
      write (*,*) x_16
      write (*,'(f30.25)') x_16
!
      x_16 = 2 * acos (x_16)
      write (*,*) x_16
      write (*,'(f30.25)') x_16
!
   end

I have looked into this and FTN95 deliberately accepts REAL16 (at least in many contexts) but goes on to treat it as REAL10.

I am not sure if we need to do anything about this. We could report a comment when compiling.

REAL16 and COMPLEX32 were perfectly acceptable in Digital VAX Fortran-77 but there were two versions of this that were switched on by compiler switches. These were: 'HUGE' which gave an huge range, i.e. the exponent and 'GRAND' wich restricted the exponent but increased the precision.

Sadly no longer easily available but I do have a Microvax 3200 hiding in a cupboard at home. I could look up the Fortran manual to provide the exact details.

As to whether it is needed, I remember one program that failed on REAL8 and this was a transfer matrix analysis of the whirling and lateral vibration of shafts. A series of matrices were multiplied together as you progress along a shaft from one end to the other multiplying 'point' matrices which contain lumped mass and inertia terms multiplied by a frequency in radians per sec squared and a 'shaft' matrtix which added in the effect of the bending stiffness of the shaft and its own distributed mass and inertia. After doing this for the entire length of the shaft, the determination of whether the guessed frequency was in fact a natural frequency was by calculating the determinant which was a series of large figures multiplied together minus a second series of large numbers multiplied together. If this is zero the it is the natural frequency. This failed on REAL8 and suceeded on REAL16 due to precision. I have sucessfully operated this program on REAL10 in FTN95. However, I think that there other ways to extract the eigenvalues of a system rather than transfer matrices which are only an step up from a Holzer frequency table.

Ian

In the early days of finite elements, Brits were able to solve problems with ill-conditioned matrices and other numerically poor algorithms that our US cousins could not, because most British computers used two 24-bit words for REALs (e.g. Elliot 4120, ICL 1900) making in effect REAL6, whereas standard precision at the time in many US computers was 32 bit (REAL4). Whether or not this promoted the drive for double and quad precision is a moot point, but precise numbers of bytes instead of words certainly came later. The CDC 6000 and 7000 computers used 60-bit words.

The original design of the 8088/8086 chip with an 8087 numeric coprocessor implemented REAL10 in the latter as a means to control roundoff in REAL8 calculations as I understand it, rather than to implement a system using REAL*10 throughout.

As Ian points out, algorithms that demand extraordinary precision may not be the best ones to use. Indeed, extraordinary precision is not required for every variable - I doubt if the diameter and length of the shaft needed specifying in REAL*8 millimetres, let alone anything higher.

OK, Now I'm home I've looked up the VAX REAL16 and found that COMPLEX32 did not exist. The REAL16 variables permitted the magnitude of values in the range of 0.84Q-4932 to 0.59Q+4932 where the Q replaces the E or D for REAL4 and REAL8 respectively. It states that 33 decimal digit precision is provided. The variation in precision of REAL8 is actually referred to as 'D_floating' and 'G_floating'. It states that the G_floating implementation offers a greater range but is less precise, having a smaller number of significant digits. D_floating is the default implementation of REAL8 and COMPLEX16, the '/G_floating' compilation switch was used to choose the optional lower precision.

The summary of the VAX Floating point number implementations is

Type    Sign  Exponent Fraction    Precision  Range                 DEC's Description
        bits  bits     bits        dec dig     
REAL*4  1     8        23          7         0.29E-38 to 1.7E+38      F_floating
REAL*8  1     8        55          16        0.29E-38 to 1.7E+38      D_Floating
REAL*8  1     11       52          15        0.56D-308 to 0.9D+308    G_floating
REAL*16 1     15       112         33        0.84Q-4932 to 0.59Q+4932 H_floating

Hence the G_floating is 'Grand' an H_floating is 'Huge' precision. FTN95 'REAL*10' has approximately the same range as H_floating, but only has 18 decimal digits precision.

I hope this helps. Ian