I would like to add my two cents to this discussion.
First, a test. Maybe all the participants could try this piece of
code to see wether there is actually a consistent implementation
of DPROD across compilers:
-----------------------------------------------------
program prod_test
double precision :: dp, de
real :: sp, se
dp = 4.0D0 *datan( 1.0D0 ) ! Double Precision Pi
de = dexp( 1.0D0 ) ! Double Precision e
sp = 4.0 *atan( 1.0 ) ! and here in Single Precision
se = exp( 1.00 )
write(*,'(A/)') " Double Precision Results:"
write(*,*) "dp*de =", dp*de
write(*,*) "dble(sp*se) =", dble(sp*se )
write(*,*) "dprod(sp*se) =", dprod(sp,se)
write(*,*) "dble(sp)*dble(se) =", dble(sp)*dble(se)
write(*,'(/A/)') " Single Precision Results:"
write(*,*) "sp*se =", sp*se
write(*,*) "real(dp*de) =", real(dp*de )
write(*,*) "real(dble(sp)*dble(se)) =", real(dble(sp)*dble(se) )
write(*,*) "real(dprod(sp,se)) =", real(dprod(sp,se))
write(*,'(/A,E10.2)') " Single Precision Epsilon ", epsilon(sp)
write(*,*) "Single Precision Nearest ", nearest(sp*se,2.0)
end program prod_test
-----------------------------------------------------
My results on a DEC Alpha with f90, version DIGITAL Fortran 90 X4.1-575:
Double Precision Results:
dp*de = 8.53973422267357
dble(sp*se) = 8.53973388671875
dprod(sp*se) = 8.53973420097986 ! <--
dble(sp)*dble(se) = 8.53973420097986 ! <--
Single Precision Results:
sp*se = 8.539734
real(dp*de) = 8.539734
real(dble(sp)*dble(se)) = 8.539734
real(dprod(sp,se)) = 8.539734
Single Precision Epsilon 0.12E-06
Single Precision Nearest 8.539735
-----------------------------------------------------
So, obviously DEC thinks DPROD(x,y) is defined as DBLE(x)*DBLE(y),
but not DBLE(x*y). BUT, the difference between both results is small
and, of course, within the accuracy that one can expect from using
single precision numbers x and y. Anyone who tries to make use of such
a difference in actual code should visit a course on numerical analysis.
The results of NEAREST for single precision reveal that the computer
can only represent the following numbers:
..., 8.539733, sp*se = 8.539734, 8.539735, ...
and the Single Precision Results are consistent with that, as one would
expect.
DPROD is an oddball among the intrinsic procedures in Fortran. No
definition in the standard seems to exist, except that x and y are
REAL and the result is DOUBLE PRECISION, a rather empty statement
for real practical use.
My solution: Make DPROD obsolete asap for the sake of a cleaner, leaner
Fortran.
Cheers,
Werner
-----------------------------------------------------------------------
| Werner W Schulz |
| Dept of Chemistry email: [log in to unmask] |
| University of Cambridge Phone: (+44) (0)1223 336 502 |
| Lensfield Road Secretary: 1223 336 338 |
| Cambridge CB2 1EW Fax: 1223 336 536 |
| United Kingdom WWW: |
-----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|