Steve Lionel of Dr Fortran fame posted a blog (
http://intel.ly/2nZSxoE) recently at the Intel software site on the KIND facilities in Fortran and he concluded, "Most applications should use the SELECTED_xxx_KIND intrinsics as I described..."
Dr Fortran blogpost was in reaction to suggestions made on this forum by some expert(s) on another recent thread which the 'doctor' didn't stop to examine fully, consequently found them objectionable, and he wrote, "we get to the point of contention I mentioned at the start of this post... Fortran 2008 extended intrinsic module ISO_FORTRAN_ENV to include named constants INT8, INT16, INT32, INT64, REAL32, REAL64 and REAL128 whose values correspond to the kinds of integer and real kinds that occupy the stated number of bits. More than one response in that newsgroup thread recommended using these instead of hard-coding integer and real kinds. In my view, this is little better than the old *n extension in that it tells you that a type fits in that many bits, but nothing else about it. As an example, there's one compiler where REAL128 is stored in 128 bits but is really the 80-bit "extended precision" real used in the old x86 floating point stack registers. If you use these you might think you're using a portable feature, but really you're not and may get bitten when the kind you get doesn't have the capabilities you need. My advice is to not use those constants unless you truly understand what they do and do not provide."
The blog thus raised a couple of alarms - quite loudly one might add - for beginning and casual or otherwise unsuspecting coders, but it did not offer any complete example or a traceable reference that illustrated how someone can actually "get bitten" or get themselves into "trouble down the road". Moreover, the blog did not present any additional insight into the nature or the probability or the frequency of the problem.
The main takeaways from the blog appeared to be that with REAL128 named constant "you might think you're using a portable feature, but really you're not and may get bitten " and that "Most applications should use the SELECTED_xxx_KIND intrinsics", the value of the recommendation to use SELECTED_REAL_KIND seems to be with portability, though it is not stated as such in the blog.
With this background, I would like to bring attention to Fortranners using a fully worked out example of a portability issue that can occur with SELECTED_REAL_KIND also. The example uses a simple computation explained in a paper by Ghazi et al., Computing in Science & Engineering, Volume 12, Issue 3, p5, May-June 2010 (ref*) in which they try "to determine 10 decimal digits
of the constant d = 173746a + 94228b − 78487c, where a = sin(1E22), b = log(17.1), and c = exp(0.42)."
[ref*
http://ieeexplore.ieee.org/document/5457292/ ]
Note the quantities and the functions involved in the example computation present a challenge in terms of "double precision" arithmetic and a coder will quickly realize the need for extended precision, say something beyond the 15 that is offered by IEEE-754 binary64 floating-point type.
Now you can try to envision a scenario where a coder using Intel Fortran compiler determines a precision of 18 is good for the need at hand and writes code as shown below and calculates a result of d = -0.1341818958E-011 (rounded to 10 significant digits) which is the same as the expected value. The programmer can justifiably conclude the code is satisfactory, it appears standard-conforming, it seems to adhere to the advice of the good doctor at least re: portability with the use of SELECTED_REAL_KIND, and the processor agrees and it does not throw any warning or errors during compilation.
-- begin code --
program p
! see for reference:
http://ieeexplore.ieee.org/document/5457292/
use, intrinsic :: iso_fortran_env, only : output_unit
implicit none
! Constants
integer, parameter :: wp = selected_real_kind( p=18 )
! Say a computational scientist guessetimates a precision of 18 is sufficient
! Calculation Constants
real(wp), parameter :: alpha = 1E22_wp
real(wp), parameter :: beta = 17.1_wp
real(wp), parameter :: gamma = 0.42_wp
real(wp), parameter :: c1 = 173746.0_wp
real(wp), parameter :: c2 = 94228.0_wp
real(wp), parameter :: c3 = 78487.0_wp
real(wp), parameter :: d_expected = -0.1341818958E-11_wp
! format strings
character(len=*), parameter :: fmt_gen = "(*(g0))"
character(len=*), parameter :: fmt_hex = "(g0,z0)"
character(len=*), parameter :: fmt_pct = "(g0,1pg16.4)"
character(len=*), parameter :: NL = new_line("")
! variables
real(wp) :: a
real(wp) :: b
real(wp) :: c
real(wp) :: d
!
write( output_unit, fmt=fmt_gen ) ProgTitle() // NL
a = sin( alpha )
b = log( beta )
c = exp( gamma )
! Computation
d = c1*a + c2*b - c3*c
! Output results
write( output_unit, fmt=fmt_gen ) "Results for working precision defined by " // &
"SELECTED_REAL_KIND( P=18 )" // NL
write( output_unit, fmt=fmt_gen ) "d = ", d
write( output_unit, fmt=fmt_hex ) "d in HEX: ", d
write( output_unit, fmt=fmt_pct ) "% error in d = ", (d/d_expected - 1.0_wp)*100.0_wp
write( output_unit, fmt=fmt_gen ) "d (expected) = ", d_expected
write( output_unit, fmt=fmt_hex ) "d (expected) in HEX: ", d_expected
stop
contains
function ProgTitle() result( Title )
!
integer :: i
integer, parameter :: Ititle(*) = &
[ 80,114,111,103,114,97,109,32,116,111,32,99,111,109,112,117,116,101,32,100,32,61, &
32,49,55,51,55,52,54,97,32,43,32,57,52,50,50,56,98,32,45,32,55,56,52,56,55,99,10, &
119,104,101,114,101,32,97,32,61,32,115,105,110,40,49,69,50,50,41,44,32,98,32,61, &
32,108,111,103,40,49,55,46,49,41,44,32,97,110,100,32,99,32,61,32,101,120,112,40, &
48,46,52,50,41,46,10,65,115,115,117,109,101,32,119,101,32,119,97,110,116,32,116, &
111,32,100,101,116,101,114,109,105,110,101,32,49,48,32,115,105,103,110,105,102, &
105,99,97,110,116,115,32,100,105,103,105,116,115,32,105,110,32,116,104,101,32,97, &
110,115,119,101,114,46,10,84,104,101,32,101,120,112,101,99,116,101,100,32,114, &
101,115,117,108,116,32,105,115,32,100,32,61,32,45,48,46,49,51,52,49,56,49,56,57, &
53,56,69,45,49,49 ]
! Function result
character(len=size(Ititle)) :: Title
forall (i=1:size(Ititle)) Title(i:i) = achar( Ititle(i) )
return
end function ProgTitle
end program p
-- end code --
Now say a need has been presented where the code needs to be ported to another compiler, say gfortran. Processing the code with gfortran will also show no warnings or errors and runtime execution will be without any exceptions. However, one will have to pay attention to notice the result using gfortran of d = -0.1314504061E-11 only provides an agreement of 4 significant digits with the expected value, compare this to 10 significant figures achieved with Intel Fortran.
Those familiar with gfortran will quickly realize the use of SELECTED_REAL_KIND with a precision of 18 matches the KIND result to the so-called KIND(10) which indeed has a precision of 18. But with Intel Fortran, the selected value is of KIND(16) with a precision of 33, for it is the only kind that can meet the minimum precision requested by the coder. It is the added precision in the chosen kind with the Intel Fortran case that allows the computed result to carry the extra accuracy in terms of significant digits in the result, but this can occur unbeknownst to the developer. A coder needs to be paying close attention to detail to notice such subtle differences, it is an aspect that indeed does get overlooked easily by many developers. Imagine now a situation where such a calculation is buried thousands of lines deep into a large, complex application and where the loss of 6 significant digits out of desired 10 in an intermediate result are impactful. It can consume countless hours of debugging exercise to uncover the problem. Of course a quick change involving only a couple of keystrokes in the SELECTED_REAL_KIND statement in the code to request a higher precision can resolve the matter, but that can only happen if the root cause can be so easily identified. I contend the situation shown here is similar to that alluded to in the above-mentioned blog with the use of the named constants from ISO_FORTRAN_ENV but here it is occurring with the touted cure, not the acursed disease in the blog.
So note in the particular example given here, if the coder had simply taken the 'shortcut' and made use of the named constant of REAL128 from ISO_FORTRAN_ENV intrinsic module, an approach that was the point of contention in the above-mentioned blog, it is highly likely the result would have been the same with both processors with d = -0.13418189578E-011 (full disclosure: this was only checked with a couple of different versions of Intel Fortran and gfortran compilers). Thus the portability exercise of going from Intel Fortran to gfortran would have been uneventful (as is desired during porting), at least for the computation at hand. This is because the computations would have obtained the few extra digits of precision required for accuracy in the calculations, regardless of how REAL128 is implemented by the processor and even if it were in a restricted manner with 80-bit 'extended precision'. Those who disagree with this are invited to provide some workable examples with sufficient detail, so readers can further their learning with precision and computations.
Note in computations that require precision beyond the so-called double precision, it's not often the number of digits of precision that can truly suffice for one's needs are known. Coders will then either pick the highest precision possible or resort to a newer paradigm involving "arbitrary precision", both of which only go in the direction of limiting the applicability of SELECTED_REAL_KIND option in end user code relative to the flexibility a coder might think there exists with precision and range selection.
I would then like readers to make a note of the following:
* recommendations come with qualifications and it is up to each "doctor" to explain the qualifications fully, for each reader to take the time and effort to notice and appreciate the qualifications an expert is striving to impart, often in a language different from one of their proficiency, and for each coder to develop a good understanding of them,
* while portability is a term bandied about a lot, one will struggle to find a fully detailed and workable definition of portability that can satisfy all the stakeholders of an application (particularly the pesky powers-that-be). Ultimately you may be responsible for defining it for your own needs, to determine the success factors for it, and how to validate your code against them.
* SELECTED_REAL_KIND is but one aspect that helps with portability, but even this can prove insufficient, as explained with the example above. Like with most things, paying attention of detail is critical to success, as can be seen here. Developing minimal working examples of crucial components of one's calculations and evaluating them thoroughly can make it easier to port code. More unit tests the better when it comes to validating code.