This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[Patch, Fortran] Another floating point speedup
- From: Daniel Kraft <d at domob dot eu>
- To: Fortran List <fortran at gcc dot gnu dot org>, gcc-patches <gcc-patches at gcc dot gnu dot org>
- Date: Sat, 03 Jan 2009 15:46:22 +0100
- Subject: [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