On Fri, 31 Jan 2003, Dick Hendrickson wrote:
> Aleksandar Donev wrote:
>>
>>
>> The code given below, which is a translation of the mathematical breakup
>> of the cases that can occur (the case C<=0 means the spheres are already
>> touching or overlapping and is a by-product of finite accuracy), runs
>> fine if I enable full IEEE handling:
>>
>> A = SUM (v_ij**2) - expansion_rate ** 2
> I'm no expert here, but that's never stopped me before. Here you
> are letting the compiler do the sums in whatever order and then
> wanting the final accuracy to be smaller than epsilon. Is that
> reasonable? Wouldn't you be better off doing these sums in
> double precision (or double whatever r_wp is) to try to
> control the roundoff? Or maybe sort the values before doing the
> sums? Or maybe trying Kahan's trick of doing the additions
> to minimize round-off problems? It just seems odd to me to
> expect essentially perfect results from this longish series of
> floating point calculations.
>
> Dick Hendrickson
Hi Dick,
I sent a pointer to Aleksander privately regarding his quadratic equation
solution, but the sum point is very good also. I would recommend both Kahan
summation *AND* sorting before summation, but Kahan summation is really just a
trick for doubling the working precision. Another way to get this effect
would be to simply use 16B precision in this situation. That may or may not
be a speed issue.
I think the main thing that we could agree on at this point is that what he is
running into is not a fortran problem per se but rather an issue of the
numerical stability of one or more parts of his algorithm. How he chooses to
deal with this would of course bring us eventually back to compiler
implementations, though ;-).
Ted
--
Ted Stern Engineering Applications
Cray Inc. office: 206-701-2182
411 First Avenue South, Suite 600 cell: 206-383-1049
Seattle, WA 98104-2860 FAX: 206-701-2500
Debuggers' motto: Frango ut patefaciam -- I break in order to reveal
|
