Hello,

Following my question about what <0.0 means for real numbers, here is a
little more ifnormation and another plea for help from those that can help.

Here is a segment of code in a simulation which calculates the collision
time of two expanding spheres moving along linear paths. It is simply a
quadratic equation that needs to be solved, and the answer is C /
(-B+Sqrt(B**2-A*C)) if the two spheres collide (i.e. if this result is
nonnegative and exists). I need to be very careful that this gives a
very accurate answer, and tricks like replacing ">0" with
">EPSILON(1.0)" have consistently lead to the simulation breaking down
(i.e. a misprediction of some collision), to which the simulation is
very sensitive, since it is an "event driven" molecular dynamics which
relies on proper ordering and prediction of the sequence of sphere
collisions.

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
         B = DOT_PRODUCT (r_ij, v_ij) - expansion_rate * current_diameter
         C = SUM (r_ij**2) - current_diameter ** 2
         IF (((B <= 0.0_r_wp) .OR. (A < 0.0_r_wp)) .AND. (B**2-A*C >=
0.0_r_wp)) THEN
            IF (C <= 0.0_r_wp) THEN
               IF (C <-Max(100.0_r_wp*EPSILON(1.0_r_wp), 1E-5_r_wp)) THEN
                  WRITE (ERROR_UNIT,*) "ERROR:", "Overlapping spheres ",
particle, " and ", partner, " detected at time ", &
                 & current_time, " at ", lsd_packing%lsd_time, ".
Overlap=", (Sqrt(SUM(r_ij**2))-current_diameter) / &
                 & current_diameter
                  lsd_in_error = .TRUE.
                  lsd_problem = (/ particle, real_partner /)
               END IF
               collision_time = current_time
            ELSE IF ((B < 0.0_r_wp) .OR. (A < 0.0_r_wp)) THEN
               collision_time = current_time + C / (-B+Sqrt(B**2-A*C))
            ELSE
               WRITE (OUTPUT_UNIT,*) "Suspicios collision of particles
", particle, " and ", partner, " not expected to collide: &
              &A, B, C=", A, B, C
               collision_time = HUGE (1.0_r_wp)
            END IF
         ELSE
            collision_time = HUGE (1.0_r_wp)
         END IF

However, as NAG is telling me, suspicious IEEE operations are happening:

Warning: Floating divide by zero occurred
Warning: Floating underflow occurred

Debugging shows that the problem is at:

Program received signal SIGFPE, Arithmetic exception.
661                    collision_time = current_time + C /
(-B+Sqrt(B**2-A*C))

Unfortunately I could not get the debugger (dbx90) to give me the values
of A, B and C, and so I am not sure what is happenning here. Other
compilers never complain about anything even if I use their switches to
enable IEEE signal trapping.

How should I write something like this in Fortran 90? Or is using the
IEEE TR module and its facilities the only way to go?

Thanks,
Aleksandar

--
__________________________________
Aleksandar Donev
Complex Materials Theory Group (http://cherrypit.princeton.edu/)
Princeton Materials Institute & Program in Applied and Computational Mathematics
@ Princeton University
Address:
   419 Bowen Hall, 70 Prospect Avenue
   Princeton University
   Princeton, NJ 08540-5211
E-mail: [log in to unmask]
WWW: http://atom.princeton.edu/donev
Phone: (609) 258-2775
Fax: (609) 258-6878
__________________________________