Dr W.W. Schulz wrote:
>
> On Tue, 23 Jun 1998, Peter Shenkin wrote:
>
> > On Jun 23, 6:22pm, Dr W.W. Schulz wrote:
> > > 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 ! <--
> > > ...
> > >
> > > So, obviously DEC thinks DPROD(x,y) is defined as DBLE(x)*DBLE(y)
> > ...
> >
> > DEC may or may not treat DPROD this way; I don't know the Alpha
> > architecture. But it is *not* obvious that it does. There is
> > nothing to suggest that sp and se are cast to double before
> > they are multiplied; on machines that support mulitplication
> > of REALs into a DOUBLE register, the cast is not needed. If
> > DEC supports this, and if their DPROD carries out that instruction
> > then this would be the expected result *** even though DPROD(x,y)
> > would not then be defined as DBLE(x)*DBLE(y) ***.
> <snip>
>
> My results above show clearly that the DEC Alpha treats the multiplication
> of single and double precision differently. Other computers may not do that
> since they treat single precision internally as double precision (which is
> legal) and do all operations in double precision internally for whatever
> reasons.On such a machine the whole discussion about DPROD is mute since
> there is no need for it.
>
> Only machines like DEC Alpha and the HP (and others, thanks to the posted
> results) it matters. And the given results seem to be consistent with the
> equivalence of DPROD(x,y) and DBLE(x)*DBLE(y). How this is achieved in detail
> is largely irrelevant. If DPROD(x,y) and DBLE(x)*DBLE(y) always give the
> same result then for all what users can know and care they can assume
> that DPROD(x,y) is "defined as" DBLE(x)*DBLE(Y). This is good enough for
> a definition.
>
> Besides this there is the question of the importance of having DPROD.
> M Olagnon wrote earlier (19 Jun 98):
> > ... we have huge sets of (single precision) data that we spend lots of
> > time correlating in every possible manner (double precision dot products).
> > If one multiplies 2 single-precision numbers, one gets twice the number
> > of significant bits.
> ^^^^^^^^^^^
> This is wrong, one gets twice the number of bits but half of them are
> insignificant. (A similar result is known to people using Fast Fourier
> Transforms.)
I am sorry, the result of the multiplication of 2 n-bit integer values
has 2n-1 bits, of which all are significant (or else, I wish I had
known that in kindergarten, it would have spared me half my efforts ;-)).
If I represent my data as real values because the products exceed
integer representation capacity, I don't see any reason why these
bits should suddenly become insignificant.
> If one starts with single precision input one cannot squeeze out
> more useful "information" than is given by a single precision number.
> One may want to restrict error propagation internally by switching to
> double precision but the final results still can only be quoted with
> single precision accuracy. The intrinsic EPSILON function is very valuable
> here.
>
> I can only repeat, the value of DPROD is zero. Where people have worry
> seriously about accuracy Fortran 90 has the relevant tools ready, and
> they should know the basic results of error analysis.
> And where the multiplication of singles in double precision registers
> is supported I am sure the compilers knows how to optimise DBLE(x)*DBLE(y).
>
I'd like to be so sure. My experience is that many compilers are not so good
at optimizing these expressions when x and y are arrays. But of course,
they should soon be all excellent ...
Michel
Michel OLAGNON email: [log in to unmask]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|