Speeding up floating point number reading

[Daniel Kraft <d at domob dot eu>]

Hi all,

as promised on PR 38654, I gave floating point reading a look.  I used 
the attached program for testing with -O3 -static.  It takes around 31s 
on my system with an unpatched gfortran 4.4.

First, I tried to replace the strtod call in 
libgfortran/io/read.c:convert_real with a custom made parser according 
to the one I implemented for FreeWRL, and even got execution time to 
around 29s (a bit surprising, I though strtod must be well optimized).

However, I then looked further, and discovered that read_f, the caller 
of convert_real, is already doing most of the work!  And then, it prints 
the number back to a string just to parse it again using convert_real... 
With the simple (experimental and rough) patch attached, the final 
floating point value can be calculated with little additional effort in 
read_f, without ever calling convert_real.  This cuts the execution time 
in half to 15.5s!

The only weak point is that I'm not sure if we may experience rounding 
errors in parsing the float due to all the *10 and + operations, and 
that could be something strtod takes care of.  On the other hand, I 
think this is rather unlikely.  What do you think about this one?

I'm going to work out a real patch, that handles all kinds, not just 
double, and in addition cleans up all now-useless stuff.  Comments welcome!

Yours,
Daniel

--
Done:  Arc-Bar-Cav-Rog-Sam-Tou-Val-Wiz
To go: Hea-Kni-Mon-Pri-Ran

Index: read.c
===================================================================
--- read.c	(revision 142940)
+++ read.c	(working copy)
@@ -792,10 +792,14 @@ read_f (st_parameter_dt *dtp, const fnod
   char *p, *buffer;
   char *digits;
   char scratch[SCRATCH_SIZE];
+  double result;
+  int dec_exp; /* Additional exponent from decimal point.  */
 
   val_sign = 1;
   seen_dp = 0;
   wu = f->u.w;
+  result = 0.0;
+  dec_exp = 0;
 
   p = gfc_alloca (wu);
 
@@ -834,7 +838,8 @@ read_f (st_parameter_dt *dtp, const fnod
   digits = p;
   ndigits = 0;
 
-  /* Scan through the string to find the exponent.  */
+  /* Scan through the string to find the exponent.  Calculate value of
+     digits on the way.  */
   while (w > 0)
     {
       switch (*p)
@@ -862,6 +867,11 @@ read_f (st_parameter_dt *dtp, const fnod
 	case '7':
 	case '8':
 	case '9':
+          result *= 10;
+          result += (*p - '0');
+          if (seen_dp)
+            --dec_exp;
+          /* Fall through */
 	case ' ':
 	  ndigits++;
 	  p++;
@@ -1008,6 +1018,8 @@ read_f (st_parameter_dt *dtp, const fnod
   if (!seen_dp)
     exponent -= f->u.real.d;
 
+  return val_sign * result * pow (10.0, exponent + dec_exp);
+
   if (exponent > 0)
     {
       edigits = 2;

! Time floating-point IO.

PROGRAM floatIO
  IMPLICIT NONE

  ! Floating-point numbers to read/write
  INTEGER, PARAMETER :: wp = 8
  REAL(KIND=wp), PARAMETER :: nums(6) = (/ 1.0_wp, .023_wp, -100.005_wp, &
                                           -.5e7_wp, 1e10_wp, &
                                           123456789000.00001_wp /)

  ! Number of repetitions.
  INTEGER, PARAMETER :: rep = 4000000

  CHARACTER(len=255) :: buffer
  REAL(KIND=8) :: readNums(SIZE (nums))
  INTEGER :: i

100 FORMAT (6F24.6)
  WRITE (buffer, 100) nums
  PRINT *, buffer

  DO i = 1, rep
    READ (buffer, 100) readNums
    !IF (i == 1 .AND. ANY (readNums /= nums)) THEN
    !  PRINT *, nums
    !  PRINT *, readNums
    !  CALL abort ()
    !END IF
  END DO
END PROGRAM floatIO