fortran 90, precision, and intrinsic functions 2

[billgray1234]

i have a few questions about a program i'm writing in fortran 90. i've written a few explanation notes below, followed by my questions.

PRECISION IN FORTRAN 90
to enable portability between different computers, i'm declaring real variables using the KIND parameter, according to

integer, parameter :: PR = selected_real_kind(p=N)
real(KIND=PR) :: x

where PR stands for 'Precision of Real numbers', x is a real variable, and N could be 4, 8, 15, etc.

INTRINSIC FUNCTIONS
i want to use intrinsic functions (such as SIN, EXP, ABS).

INTRINSIC FUNCTIONS AND PRECISION
to my understanding, each of the above intrinsic functions has several 'versions', depending on the level of precision. these versions include 1) single precision, 2) double precision, and, on some computers, 3) even higher precision. for example:-
SIN(x) -- single precision is SIN(x), double precision is DSIN(x)
EXP(x) -- single precision is EXP(x), double precision is DEXP(x)
and so on.
i think 'single precision' functions are called GENERIC functions, whereas 'double precision' functions are called SPECIFIC functions. also, to my understanding, the 'level of precision' in these intrinsic functions refers to the level of precision of the argument (i.e. x) and the result computed (please correct me if i'm wrong about this).

MY SITUATION
let's say i'm writing my program on computer 1. i know 'in advance' what level of precision computer 1 is capable of. for example, let's say it's capable of higher precision than just single precision (e.g. N=15 above), such that it's capable of using the double precision versions of the intrinsic functions (i.e. DSIN, DEXP, DABS, etc).
but, now i want to compile/run my program on another computer (computer 2, say). i don't know 'in advance' what level of precision computer 2 is capable of -- it might only be capable of a much lower level precision than computer 1 is capable of (e.g. it might only support N=4 above).

MY QUESTIONS
1) when i write my program on computer 1, which 'versions' of the intrinsic functions should i use -- single precision (i.e. SIN, EXP, ABS, etc) or double precision (i.e. DSIN, DEXP, DABS, etc)?
2) in general, if i don't know 'in advance' what level of precision a computer is capable of, then which 'versions' of the intrinsic functions should i use -- single precision or double precision?
3) if x is declared as being high precision (e.g. N=15), but i only use single precision intrinsic functions (i.e. SIN, EXP, ABS, etc), then will each function 'result' be single precision or double precision? in other words, is the level of precision controlled by the variable declaration statement -- or by the intrinsic function?

any help would be appreciated.

 

[xwb]

When you say portability between different computers do you mean Windows/Linux or 16, 32, 64 and possibly 128 bit?

 

[FJacq]

You misunderstand the generic functions SIN, EXP, LOG ...

They are NOT single precision functions but, on the contrary, they are valid for ANY kind of real values.

So I advise to use NEVER the functions DSIN, DEXT ... which are only valid for double precision values.

The rule is simple : SIN(alpha) return a real value having the same kind as its argument.

Example :

program test
  implicit none
  integer,parameter :: sp=selected_real_kind(6)
  integer,parameter :: dp=selected_real_kind(13)
  integer,parameter :: ep=selected_real_kind(18)
  write(*,*) 4*atan(1._sp)
  write(*,*) 4*atan(1._dp)
  write(*,*) 4*atan(1._ep)
end program

With the result :

[lcoul@localhost test]$ ifort t46.f90
[lcoul@localhost test]$ ./a.out
   3.141593    
   3.14159265358979     
   3.14159265358979323846264338327950

François Jacq

 

[mikrom]

FJacq,
Thanks for the explanation.

 

[billgray1234]

sorry it has taken me this long to reply. i've been away from my computer for a week.

firstly, thanks goes to FJacq for that explanation. basically, my confusion on the topic (i.e. my misunderstanding of the intrinsic functions) has been fuelled, in part, by the fact that most of what i said above was recited (word for word) from every book i've read on the topic. so, i guess those books are not accurate! but, FJacq's explanation is both accurate and a valuable lesson learned.

secondly, as for xwb's question about 'portability'. to my understanding (again -- please correct me if i'm wrong), the portability refers solely to PRECISION -- or rather, the PRECISION CAPABILITY of the computer's PROCESSOR. some computers are capable of higher (or lower) precision than others. for example, computer 1 might only be capable of SINGLE precision, whereas computer 2 might be capable of DOUBLE (i.e. higher) precision. therefore, if a code is written on computer 2, say, and then you want to compile/run it on computer 1, then the code needs to be modified (in order to allow for the difference in precision between the two computers). basically, in fortran 90, the portability function SELECTED_REAL_KIND (that i mentioned above) aims to MINIMIZE the total number of changes that need to be made. some examples of how this is so are given below. sorry it's a bit messy!


SINGLE PRECISION
declaring a real number --------------------------------> REAL :: X
assigning a variable value -----------------------------> X = 1.0
intrinsic functions ------------------------------------> SIN(X)
converting integer number (I) into a real number (X) ---> X = REAL(I)


DOUBLE PRECISION
declaring a real number --------------------------------> DOUBLE PRECISION :: X
assigning a variable value -----------------------------> X = 1.0D0
intrinsic functions ------------------------------------> DSIN(X)
converting integer number (I) into a real number (X) ---> X = DBLE(I)


ANY PRECISION (FORTRAN 90)
declaring a real number --------------------------------> REAL(KIND=PR) :: X
assigning a variable value -----------------------------> X = 1.0_PR
intrinsic functions ------------------------------------> SIN(X) (as FJacq clarified, above)
converting integer number (I) into a real number (X) ---> X = REAL(I,KIND=PR)


where PR=SELECTED_REAL_KIND(P=N), with the value of N representing the number of decimal digits of precision. for example, N=8 for single precision, N=15 for double precision, N=32 for higher precision, etc. (note that these are just 'hypothetical' values of N -- i don't know for certain if N=8 does represent single precision, or if N=15 does represent double precision, or if N=32 does represent higher precision. you would have to look them up).

anyway, in the FORTRAN 90 code above, in order to 'transport' the program between the different computers, all that needs to be changed is the value of N (because this then automatically changes the value of PR). therefore, only minimal changes need to be made to the code.

 

[FJacq]

Your question about portability has no easy answer. One has to use SELECTED_REAL_KIND to ask for the precision the algorithms need. If a particular couple processor/compiler cannot provides that precision, then the program will not compile at all. Example :

program test
  integer,parameter :: ep=selected_real_kind(50)
  real(ep) :: a
  a=5
  write(*,*) a
end program

Result

[lcoul@localhost test]$ ifort t50.f90
t50.f90(3): error #6684: This is an incorrect value for a kind type parameter in this context.   [EP]
  real(ep) :: a
-------^
compilation aborted for t50.f90 (code 1)

The reason is simple : if a compiler cannot insure a real kind with the expected precision, then the function selected_real_kind returns the value -1, which cannot be used to define variables or constants.

Now, if two couples processor/compiler are able to compile a program, it does not mean that thy will give exactly the same answer :
- the algorithm may be wrong,
- the programmer may be not aware about the exact precision needed by the algorithm,
- all codes have bugs,
- rounding errors are inherent to the use of real values and the result may also depend on the optimization level chosen when compiling the program (an optimizer may change the order of operations),
- the real precision used by a couple (processor,compiler) may be greater than another one because selected_real_kind must just return a kind able to provide AT LEAST the expected precision.

Example of the last point :

program test
  implicit none
  integer,parameter :: ep=selected_real_kind(18)
  write(*,*) ep
  write(*,*) 4*atan(1._ep)
end program
[code]

Two compilers on the same processor (x86) :

[code]
[lcoul@localhost test]$ ifort t51.f90
[lcoul@localhost test]$ ./a.out
          16
   3.14159265358979323846264338327950      
[lcoul@localhost test]$ gfortran t51.f90
[lcoul@localhost test]$ ./a.out
          10
   3.1415926535897932385

The precision of ifort is greater because ifort has used a REAL(16) (128bits) whereas gfortran has used a REAL(10) (80bits) available natively on the x86 processor.

Few other remarks :

- the precision of a single precision real depends on the couple processor/compiler. But on most usual systems today, all single precision real values are codes on 4 bytes with a precision of 6-7 digits. They are exceptions like Cray computers but very few people may use them...

- when writing an instruction like "x=i", one rarely use the function REAL because integers may be represented exactly by a real value if this integer is not too big. In particular, I don't now any couple processor/compiler which cannot store a default kind integer value into a double precision real value (no rounding error).

François Jacq

 

[ArkM]

>I don't now any couple processor/compiler which cannot store a default kind integer value into a double precision real value (no rounding error).

Strictly speaking, most of 64-bit compilers allow to redefine (promote) default INTEGER precision to 64-bit. It's impossible to store all such integer values into a double precision w/o rounding error.

 

[FJacq]

Strictly speaking, most of 64-bit compilers allow to redefine (promote) default INTEGER precision to 64-bit. It's impossible to store all such integer values into a double precision w/o rounding error.

Except that such promotion is not standard conforming. The Fortran norm clearly precises that :
- the default integer occupies the same memory of the single precision real,
- the double precision real occupies twice the memory of a single precision real.

In addition, Fortran2003 warrants the existence of an extended integer (selected_int_kind(10)) which makes this trick useless.

François Jacq