This is the mail archive of the gcc-patches@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[Patch, Fortran] Another floating point speedup


Hi all,

again a patch to rework libgfortran's read_f for formatted reading of floating point numbers. This one should be equivalent to the old code with respect to precision, as we're still using convert_real and utimatly strtod for the actual conversion.

This just reworks the way that read_f reformats the string to handle the Fortran/format specific issues. It is not as "good" as the other patch posted here, but it still brings down execution times from 33.4s to 25s for the test attached on my system. Profiling shows that most of the time (5s) is spent in strtod, about 3s in read_f for the reformatting and some 1s and below in other helper routines; so I'd guess that there's not much further speed-up possible (at least no dramatic one, except we could somehow get rid of this preprocessing completely) with still using strtod; but that's just a wild guess.

io/list_read.c contains two other similar routines (parse_float and read_float) that sound as if we could do some cleaning up and maybe speed-up for those three (like combining common code or at least parse_float and read_float for unformatted IO, as I see it), I can work on this as a new patch or extension to this one. What do you think?

Regression-testing at the moment on GNU/Linux-x86-32. Maybe I've got introduced some bugs with BLANK_ZERO/BLANK_NULL again, but those should not interfer much with the timing figure. Ok for 4.5 (at least from the main implementation and rather than my other patch)? I can then work on the others, if it is (or we take this one on hold and I work for a seperate patch on the rest).

Daniel

--
Done:  Arc-Bar-Cav-Rog-Sam-Tou-Val-Wiz
To go: Hea-Kni-Mon-Pri-Ran
2009-01-03  Daniel Kraft  <d@domob.eu>

	PR fortran/38654
	* io/read.c (read_f): Reworked to speed up floating point parsing.
Index: libgfortran/io/read.c
===================================================================
--- libgfortran/io/read.c	(revision 142940)
+++ libgfortran/io/read.c	(working copy)
@@ -785,73 +785,79 @@ read_f (st_parameter_dt *dtp, const fnod
 {
   size_t wu;
   int w, seen_dp, exponent;
-  int exponent_sign, val_sign;
-  int ndigits;
-  int edigits;
-  int i;
-  char *p, *buffer;
-  char *digits;
-  char scratch[SCRATCH_SIZE];
+  int exponent_sign;
+  char* p;
+  char* buffer;
+  char* out;
+  int seen_int_digit; /* Seen a digit before the decimal point?  */
+  int seen_dec_digit; /* Seen a digit after the decimal point?  */
 
-  val_sign = 1;
   seen_dp = 0;
+  seen_int_digit = 0;
+  seen_dec_digit = 0;
+  exponent_sign = 1;
+  exponent = 0;
   wu = f->u.w;
 
+  /* Read in the next block.  */
   p = gfc_alloca (wu);
-
   if (read_block_form (dtp, p, &wu) == FAILURE)
     return;
-
   w = wu;
 
   p = eat_leading_spaces (&w, p);
   if (w == 0)
     goto zero;
 
-  /* Optional sign */
+  /* In this buffer we're going to re-format the number cleanly to be parsed
+     by convert_real in the end; this assures we're using strtod from the
+     C library for parsing and thus probably get the best accuracy possible.
+     This process may add a '+0.0' in front of the number as well as change the
+     exponent because of an implicit decimal point or the like.  Thus allocating
+     strlen ("+0.0e-100") == 9 characters plus one for NUL more than the
+     original buffer had should be enough.  */
+  buffer = gfc_alloca (w + 10);
+  out = buffer;
 
+  /* Optional sign */
   if (*p == '-' || *p == '+')
     {
       if (*p == '-')
-        val_sign = -1;
-      p++;
-      w--;
+	*(out++) = '-';
+      ++p;
+      --w;
     }
 
-  exponent_sign = 1;
   p = eat_leading_spaces (&w, p);
   if (w == 0)
     goto zero;
 
-  /* A digit, a '.' or a exponent character ('e', 'E', 'd' or 'D')
-     is required at this point */
-
-  if (!isdigit (*p) && *p != '.' && *p != ',' && *p != 'd' && *p != 'D'
-      && *p != 'e' && *p != 'E')
-    goto bad_float;
-
-  /* Remember the position of the first digit.  */
-  digits = p;
-  ndigits = 0;
-
-  /* Scan through the string to find the exponent.  */
+  /* Process the mantissa string.  */
   while (w > 0)
     {
       switch (*p)
 	{
 	case ',':
-	  if (dtp->u.p.current_unit->decimal_status == DECIMAL_COMMA
-               && *p == ',')
-	    *p = '.';
-	  else
+	  if (dtp->u.p.current_unit->decimal_status != DECIMAL_COMMA)
 	    goto bad_float;
-	  /* Fall through */
+	  /* Fall through.  */
 	case '.':
 	  if (seen_dp)
 	    goto bad_float;
+	  if (!seen_int_digit)
+	    *(out++) = '0';
+	  *(out++) = '.';
 	  seen_dp = 1;
-	  /* Fall through */
+	  break;
 
+	case ' ':
+	  if (dtp->u.p.blank_status == BLANK_ZERO)
+	    *p = '0';
+	  else if (dtp->u.p.blank_status == BLANK_NULL)
+	    break;
+	  else
+	    goto done;
+	  /* Fall through.  */
 	case '0':
 	case '1':
 	case '2':
@@ -862,207 +868,165 @@ read_f (st_parameter_dt *dtp, const fnod
 	case '7':
 	case '8':
 	case '9':
-	case ' ':
-	  ndigits++;
-	  p++;
-	  w--;
+	  *(out++) = *p;
+	  if (!seen_dp)
+	    seen_int_digit = 1;
+	  else
+	    seen_dec_digit = 1;
 	  break;
 
 	case '-':
 	  exponent_sign = -1;
-	  /* Fall through */
-
+	  /* Fall through.  */
 	case '+':
-	  p++;
-	  w--;
 	  goto exp2;
 
-	case 'd':
 	case 'e':
-	case 'D':
 	case 'E':
-	  p++;
-	  w--;
+	case 'd':
+	case 'D':
+	  ++p;
+	  --w;
 	  goto exp1;
 
 	default:
 	  goto bad_float;
 	}
-    }
-
-  /* No exponent has been seen, so we use the current scale factor */
-  exponent = -dtp->u.p.scale_factor;
-  goto done;
-
- bad_float:
-  generate_error (&dtp->common, LIBERROR_READ_VALUE,
-		  "Bad value during floating point read");
-  next_record (dtp, 1);
-  return;
-
-  /* The value read is zero */
- zero:
-  switch (length)
-    {
-      case 4:
-	*((GFC_REAL_4 *) dest) = 0;
-	break;
-
-      case 8:
-	*((GFC_REAL_8 *) dest) = 0;
-	break;
-
-#ifdef HAVE_GFC_REAL_10
-      case 10:
-	*((GFC_REAL_10 *) dest) = 0;
-	break;
-#endif
 
-#ifdef HAVE_GFC_REAL_16
-      case 16:
-	*((GFC_REAL_16 *) dest) = 0;
-	break;
-#endif
-
-      default:
-	internal_error (&dtp->common, "Unsupported real kind during IO");
+      ++p;
+      --w;
     }
-  return;
+  
+  /* No exponent has been seen, so we use the current scale factor.  */
+  exponent = - dtp->u.p.scale_factor;
+  goto done;
 
-  /* At this point the start of an exponent has been found */
- exp1:
-  while (w > 0 && *p == ' ')
+  /* At this point the start of an exponent has been found.  */
+exp1:
+  p = eat_leading_spaces (&w, p);
+  if (*p == '-' || *p == '+')
     {
-      w--;
-      p++;
+      if (*p == '-')
+	exponent_sign = -1;
+      ++p;
+      --w;
     }
 
-  switch (*p)
-    {
-    case '-':
-      exponent_sign = -1;
-      /* Fall through */
-
-    case '+':
-      p++;
-      w--;
-      break;
-    }
+  /* At this point a digit string is required.  We calculate the value
+     of the exponent in order to take account of the scale factor and
+     the d parameter before explict conversion takes place.  */
+exp2:
 
   if (w == 0)
     goto bad_float;
 
-  /* At this point a digit string is required.  We calculate the value
-     of the exponent in order to take account of the scale factor and
-     the d parameter before explict conversion takes place. */
- exp2:
-  /* Normal processing of exponent */
-  exponent = 0;
   if (dtp->u.p.blank_status == BLANK_UNSPECIFIED)
     {
       while (w > 0 && isdigit (*p))
-        {
-          exponent = 10 * exponent + *p - '0';
-          p++;
-          w--;
-        }
-        
+	{
+	  exponent *= 10;
+	  exponent += *p - '0';
+	  ++p;
+	  --w;
+	}
+	
       /* Only allow trailing blanks */
-
       while (w > 0)
-        {
-          if (*p != ' ')
-  	  goto bad_float;
-          p++;
-          w--;
-        }
+	{
+	  if (*p != ' ')
+	    goto bad_float;
+	  ++p;
+	  --w;
+	}
     }    
   else  /* BZ or BN status is enabled */
     {
-      while (w > 0 && (isdigit (*p) || *p == ' '))
-        {
-          if (*p == ' ')
-            {
-	      if (dtp->u.p.blank_status == BLANK_ZERO) *p = '0';
-	      if (dtp->u.p.blank_status == BLANK_NULL)
-                {
-                  p++;
-                  w--;
-                  continue;
-                }
-            }
-          else if (!isdigit (*p))
-            goto bad_float;
-
-          exponent = 10 * exponent + *p - '0';
-          p++;
-          w--;
-        }
-    }
+      while (w > 0)
+	{
+	  if (*p == ' ')
+	    {
+	      if (dtp->u.p.blank_status == BLANK_ZERO)
+		*p = '0';
+	      else if (dtp->u.p.blank_status == BLANK_NULL)
+		{
+		  ++p;
+		  --w;
+		  continue;
+		}
+	    }
+	  else if (!isdigit (*p))
+	    goto bad_float;
 
-  exponent = exponent * exponent_sign;
+	  exponent *= 10;
+	  exponent += *p - '0';
+	  ++p;
+	  --w;
+	}
+    }
 
- done:
+done:
   /* Use the precision specified in the format if no decimal point has been
      seen.  */
   if (!seen_dp)
     exponent -= f->u.real.d;
 
-  if (exponent > 0)
-    {
-      edigits = 2;
-      i = exponent;
-    }
-  else
-    {
-      edigits = 3;
-      i = -exponent;
+  /* Output a trailing '0' after decimal point if not yet found.  */
+  if (seen_dp && !seen_dec_digit)
+    *(out++) = '0';
+
+  /* Print out the exponent to finish the reformatted number.  */
+  if (exponent != 0)
+    {
+      *(out++) = 'e';
+      if (exponent_sign < 0)
+	*(out++) = '-';
+      while (exponent != 0)
+	{
+	  *(out++) = (char) ('0' + (exponent % 10));
+	  exponent /= 10;
+	}
     }
+  *(out++) = '\0';
 
-  while (i >= 10)
-    {
-      i /= 10;
-      edigits++;
-    }
+  /* Do the actual conversion.  */
+  convert_real (dtp, dest, buffer, length);
 
-  i = ndigits + edigits + 1;
-  if (val_sign < 0)
-    i++;
-
-  if (i < SCRATCH_SIZE) 
-    buffer = scratch;
-  else
-    buffer = get_mem (i);
-
-  /* Reformat the string into a temporary buffer.  As we're using atof it's
-     easiest to just leave the decimal point in place.  */
-  p = buffer;
-  if (val_sign < 0)
-    *(p++) = '-';
-  for (; ndigits > 0; ndigits--)
+  return;
+
+  /* The value read is zero.  */
+zero:
+  switch (length)
     {
-      if (*digits == ' ')
-        {
-	  if (dtp->u.p.blank_status == BLANK_ZERO) *digits = '0';
-	  if (dtp->u.p.blank_status == BLANK_NULL)
-            {
-              digits++;
-              continue;
-            } 
-        }
-      *p = *digits;
-      p++;
-      digits++;
-    }
-  *(p++) = 'e';
-  sprintf (p, "%d", exponent);
+      case 4:
+	*((GFC_REAL_4 *) dest) = 0.0;
+	break;
 
-  /* Do the actual conversion.  */
-  convert_real (dtp, dest, buffer, length);
+      case 8:
+	*((GFC_REAL_8 *) dest) = 0.0;
+	break;
 
-  if (buffer != scratch)
-     free_mem (buffer);
+#ifdef HAVE_GFC_REAL_10
+      case 10:
+	*((GFC_REAL_10 *) dest) = 0.0;
+	break;
+#endif
 
+#ifdef HAVE_GFC_REAL_16
+      case 16:
+	*((GFC_REAL_16 *) dest) = 0.0;
+	break;
+#endif
+
+      default:
+	internal_error (&dtp->common, "Unsupported real kind during IO");
+    }
+  return;
+
+bad_float:
+  generate_error (&dtp->common, LIBERROR_READ_VALUE,
+		  "Bad value during floating point read");
+  next_record (dtp, 1);
+  return;
 }
 
 
! 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

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]