Well tested should include:
One can follow one of three paths:
w is the length of the string - then the space would be ignored and the result would be -1.2As we seek (AFAIU) for a fast replacement of read(string, *) x, I’d vote for no. 2, with possible optional argument returning info on whether the whole text was processed or not (as C’s strtod() returns via its second argument).
I think that there is little difference caused by using Gfortran 10.2 instead of 12. Rather, the differences are attributable to the different RTLs used by Cygwin and MinGW64. I compiled the source to an object file using Gfortran Cygwin 11.2. I then linked that object file using Gfortran 12.0 (MinGW64) and the resulting EXE ran in 6.6 seconds (compared to 30 s with Cygwin).
The second block of FortranFan’s code of #4 exposes a performance bug in the MinGW64 library function strtod() on Windows. Building the program with Cygwin and Gfortran 11.2 (Cygwin) gives a run time of 0.26 seconds for the second block. With Gfortran 12.0 (MinGW64) and with NAG Fortran, which uses MINGW, the second block takes “forever”, i.e., I got tired of waiting and aborted the run.
By adding a PRINT statement in the DO loop of the second block, I found that the performance penalty is enormous. The one million iterations of the DO loop take 0.5 s with Gfortran 11.2 (Cygwin), and I estimate 9000 s for NAG FortranMinGW (or Gfortran + MinGW), based on my observed 9 seconds for every 1000 iterations. That is a ratio of MinGW/Cygwin of 18,000 for run times!
Correction/Retraction: Sorry, there is a bug in the example code that I should have caught and corrected before blaming the library routine strtod(). The C library routine expects null-terminated strings. FortranFan’s code in #4 should have done something similar to
strings(i) = trim(strings(i)) // char(0)
or arranged for null-termination of strings(i) in some other way. Or, have I missed some aspects of the C-interoperability features that provides a null termination to strings?
It can take a long time looking for nulls that were not placed in the appropriate places.
Based on your last two messages, it looks like a Windows Fortran programmer trying to optimize should be willing to compile her program with both MinGW and Cygwin versions of gfortran (and also WSL 2 and Intel Fortran).
Equivalencing variables of different types (integers and characters, in your code) has been illegal since Fortran 95. EQUIVALENCE itself is marked “obsolescent” in Fortran 2018
You’re right, but I don’t know any alternative which can compete with the speed of equivalence. transfer is much slower.
Some more data points: the NAIF SPICELIB has routines nparsd (string to double) and nparsi (string to integer) with decades of heritage. To see the code, you have to download and unzip the library here. I have found that nparsd doesn’t always quite give the same answer to the last bit as the built in Fortran read.
Yes. It’s true. I checked.
There are in this post some useful codes converting string to double fast and correctly. I’d like to pack everything in one module. Here I provide the module, testing program for
Herein we can check the elemental attribution, which allows calling routine with array.
Result:
1. with Ifort (2021.6.0 20220226):
N = 1000000 , i0 = 353352
Write: time consumed = 0.3052 seconds
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.3692 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1046 seconds
Read: time consumed = 0.1037 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 3: F st_to_r8
Read: time consumed = 0.0868 seconds
Read: time consumed = 0.0861 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
BLOCK 4: F str2real
Read: time consumed = 0.0395 seconds
Read: time consumed = 0.0390 seconds (array)
A = Max |rval(:)-rref(:)| = 4.424E+06, B = epsilon(1.0d0) = 2.220E-16, B-A = -4.424E+06
*** WARNING: A > B, (A-B)/epsilon = 1.992E+22
2. with Gfortran (8.3.0)
N = 1000000 , i0 = 727798
Write: time consumed = 0.8102 seconds
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.6475 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1046 seconds
Read: time consumed = 0.1038 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 3: F st_to_r8
Read: time consumed = 0.0790 seconds
Read: time consumed = 0.0797 seconds (array)
A = Max |rval(:)-rref(:)| = 1.887E-15, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.665E-15
*** WARNING: A > B, (A-B)/epsilon = 7.500E+00
BLOCK 4: F str2real
Read: time consumed = 0.0416 seconds
Read: time consumed = 0.0412 seconds (array)
A = Max |rval(:)-rref(:)| = 9.649E+05, B = epsilon(1.0d0) = 2.220E-16, B-A = -9.649E+05
*** WARNING: A > B, (A-B)/epsilon = 4.345E+21
Source code:
!------------------------------------------------------------------------------
!
module mod_st_to_dp
!
use, intrinsic :: iso_c_binding, only: c_double, c_char, c_ptr, &
c_null_ptr, c_long
!
use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int64
!
implicit none
!
private
!
public :: f_st_to_r8
public :: c_st_to_r8
public :: str2real
!
!
interface
pure &
function c_strtod ( str, endptr ) result(d) bind(C, name="strtod")
!
!! <stdlib.h> :: double strtod(const char *str, char **endptr)
!
import :: c_char, c_ptr, c_double
!! Argument list
character(kind=c_char,len=1), intent(in) :: str(*)
!! type(c_ptr), intent(inout) :: endptr !<- since we need pure attribution.
type(c_ptr), intent(in) :: endptr
!! function result
real(c_double) :: d
end function
end interface
!
!
! integer, parameter :: dp = selected_real_kind( p=12 )
! FP constants
real(dp), parameter :: ten = 10.0_dp
! Integer constants
integer, parameter :: ascii_negative = 45
integer, parameter :: ascii_period = 46
integer, parameter :: ascii_0 = 48
integer, parameter :: ascii_9 = 57
integer(i8), parameter :: ten_int = 10_i8
!
type(c_ptr),parameter :: endptr = c_null_ptr
!
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
elemental &
subroutine c_st_to_r8 ( str, r )
character(len=*), intent(in) :: str
real(dp), intent(out) :: r
!
r = c_strtod ( str, endptr )
return
end subroutine
!------------------------------------------------------------------------------
!
elemental &
subroutine f_st_to_r8 ( str, r )
character(len=*), intent(in) :: str
real(dp), intent(out) :: r
! local variables
integer(i8) :: n, n_exp
integer :: lens, pos_exp, expnt
!
r = 0.0_dp
n_exp = 0
lens = len_trim( str )
pos_exp = index( str, "E", back=.true. ) ! <-- find the capital 'E'
if ( pos_exp == 0 ) then
pos_exp = index( str, "e", back=.true. )
endif
if ( pos_exp > 0 ) then
call str2dec( str(pos_exp+1:lens), n_exp )
else
pos_exp = lens + 1
endif
call str2dec( str(1:pos_exp-1), n, expnt )
r = real( n, kind=dp ) * ( ten**(n_exp + expnt) )
return
end subroutine
!
elemental &
subroutine str2dec( s, n, expnt )
character(len=*), intent(in) :: s
integer(i8), intent(out) :: n
integer, intent(out), optional :: expnt
integer :: ipos, ic, pos_period, iexponent
!
n = 0_i8
pos_period = 0
iexponent = -1
if ( present(expnt) ) iexponent = 0
do ipos = len(s), 1, -1
ic = ichar( s(ipos:ipos) )
if ( present(expnt) ) then
if ( ic == ascii_period ) then
pos_period = ipos
cycle
endif
endif
if ( ic == ascii_negative ) then
n = -n
exit
endif
if ( (ic < ascii_0) .or. (ic > ascii_9) ) exit
iexponent = iexponent + 1
n = n + (ten_int**iexponent)*(ic - 48)
enddo
if ( present(expnt) ) expnt = -len(s) + pos_period - 1
return
end subroutine
!------------------------------------------------------------------------------
!
elemental &
function str2real ( s ) result ( r )
!
character(*), target, intent(in) :: s
real(dp) :: r
!
integer, parameter :: N = 32
real(dp), parameter :: expo(N) = &
[1e8, 1e7, 1e6, 1e5, 1e4, 1e3, 1e2, 1e1, 1e0, &
1e0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, &
1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 1e-16, &
1e-17, 1e-18, 0., 0., 0., 0.]
character(N) :: equ_s
integer*1 :: equ_i(N)
integer*1 :: mask(N)
equivalence(equ_i, equ_s)
equ_s = s
equ_i = equ_i - 48
where (0 <= equ_i .and. equ_i <= 9)
mask = 1
elsewhere
mask = 0
end where
equ_i = equ_i * mask
r = sum(equ_i * expo) * 10**(equ_i(29)*10+equ_i(30))
!
return
end function str2real
!------------------------------------------------------------------------------
end module
!------------------------------------------------------------------------------
!======================================================================
!
program main
!
use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int64
use mod_st_to_dp
!
implicit none
!
integer,parameter :: n = 1000000 !! number of values
!
character(len=30),dimension(:),allocatable :: strs
real(dp),dimension(:),allocatable :: rval, rref
!
integer :: i, i0
integer :: ierr
real(dp) :: r
integer(i8) :: start, finish, count_rate
real(dp),parameter :: eps = epsilon(1.0_dp)
!
!
! create a list of values to parse
!
allocate( strs(n), rval(n), rref(n) )
!
! Initializing RNG: with a dummy i
r = util_unirand (i)
i0 = int(n*r) + 1
write(*,*) 'N =', n, ', i0 =', i0
do i = 1, n
rref(i) = util_unirand()
enddo
call system_clock(start, count_rate)
do i = 1, n
write(strs(i), '(E30.16)') rref(i)
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Write: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
blk_io: block
write(*,11) "BLOCK 1: formatted read toward string to double"
!
call system_clock(start, count_rate)
do i = 1, n
read(strs(i),fmt=*,iostat=ierr) rval(i)
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during formatted read."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_io
!! endif
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_io
blk_c: block
write(*,11) "BLOCK 2: C st_to_r8"
!
call system_clock(start, count_rate)
do i = 1, n
call c_st_to_r8 ( strs(i), rval(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during st_to_r8."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_c
!! endif
!
call system_clock(start, count_rate)
call c_st_to_r8 ( strs, rval )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_c
blk_f: block
write(*,11) "BLOCK 3: F st_to_r8"
!
call system_clock (start, count_rate)
do i = 1, n
call f_st_to_r8 ( strs(i), rval(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during st_to_r8."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_f
!! endif
!
call system_clock(start, count_rate)
call f_st_to_r8 ( strs, rval )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_f
blk_4: block
write(*,11) "BLOCK 4: F str2real "
!
call system_clock (start, count_rate)
do i = 1, n
rval(i) = str2real ( strs(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
call system_clock(start, count_rate)
rval = str2real ( strs )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_4
!
10 format ( 'i =', i9, ', string(i) =', a, ', zval(i) =', e25.17 )
11 format ( /, 80('-'), /, a )
!
deallocate ( strs, rval, rref )
!
STOP
!
contains
!======================================================================
!
! This is the ran0 in Numerical Rrecipes in Fortran.
! To initialize,
! call system_clock( idummy )
! Do not alter idummy between successive deviates.
!
function util_unirand ( idum ) result( ran0 )
!
implicit none
integer,parameter :: rk = kind(1.0D0), ik = 4
!
! ARGUMENTS:
!
real(rk) :: ran0
integer(ik),intent(in),optional :: idum
!
! local constants and variable
!
integer(ik) :: k
integer(ik),save :: idummy
integer(ik),parameter :: &
ia=16807, im=2147483647, iq=127773, ir=2836, mask=123459876
real(rk),parameter :: am = 1.0_rk / im
!
!
if ( present(idum) ) call system_clock( idummy )
!
idummy = ieor( idummy, mask )
k = idummy / iq
idummy = ia*( idummy - k*iq ) - ir*k
if ( idummy .lt. 0 ) idummy = idummy + im
ran0 = am*idummy
idummy = ieor( idummy, mask )
return
end function
!
!=====
subroutine show_err ( er, ep )
real(dp),intent(in) :: er, ep
real(dp) :: tmp
!
tmp = ep - er
write(*,10) er, ep, tmp
if ( tmp .lt. 0.0_dp ) write(*,11) -tmp/ep
return
10 format( 'A = Max |rval(:)-rref(:)| =', 1pe10.3,',', 3x, &
'B = epsilon(1.0d0) =', 1pe10.3,',', 3x, &
'B-A =', 1pe11.3 )
11 format( '*** WARNING: A > B, (A-B)/epsilon =', 1pe11.3 )
end subroutine
!
!======================================================================
end program
!======================================================================
Conclusion: (obsolete)
==============================================================
UPDATE 1
Based on Carltoffel’s comment, I can check his current version only with Ifort. My gfortran can not compile the code since it needs to be upgraded to the 9.x versions.
Results (with Ifort)
N = 1000000 , i0 = 84194
Write: time consumed = 0.3079 seconds
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.3865 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1036 seconds
Read: time consumed = 0.1049 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 3: F st_to_r8
Read: time consumed = 0.0943 seconds
Read: time consumed = 0.0970 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
BLOCK 4: F str2real
Read: time consumed = 0.0609 seconds
Read: time consumed = 0.0605 seconds (array)
A = Max |rval(:)-rref(:)| = 3.331E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.110E-16
*** WARNING: A > B, A/B = 1.500E+00
Conclusion: (only from my side)
Source code of str2real_m.f90, in which I added the elemental attribution to call the routine with array, as in the test program. To use the new str2real, we must add the line use str2real_m to the test program.
module str2real_m
implicit none
private
public :: str2real
integer, parameter :: wp = kind(1.0d0)
integer, parameter :: ikind = selected_int_kind(2)
integer(kind=ikind), parameter :: digit_0 = ichar('0')
integer(kind=ikind), parameter :: period = ichar('.') - digit_0
integer(kind=ikind), parameter :: minus_sign = ichar('-') - digit_0
real(wp), parameter :: whole_number_base(*) = &
[1d15, 1d14, 1d13, 1d12, 1d11, 1d10, 1d9, 1d8, &
1d7, 1d6, 1d5, 1d4, 1d3, 1d2, 1d1, 1d0]
real(wp), parameter :: fractional_base(*) = &
[1d-1, 1d-2, 1d-3, 1d-4, 1d-5, 1d-6, 1d-7, 1d-8, &
1d-9, 1d-10, 1d-11, 1d-12, 1d-13, 1d-14, 1d-15, 1d-16, &
1d-17, 1d-18, 1d-19, 1d-20, 1d-21, 1d-22, 1d-23, 1d-24]
real(wp), parameter :: period_skip = 0d0
real(wp), parameter :: base(*) = &
[whole_number_base, period_skip, fractional_base]
contains
elemental &
function str2real(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_coefficient
real(wp) :: r_exponent
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer(kind=ikind) :: mask(N)
integer :: period_loc
integer :: exponent_loc
integer :: mask_from
integer :: mask_till
intrinsic :: findloc, scan
integer :: findloc, scan
equivalence(factors, factors_char)
factors_char = s
factors = factors - digit_0
period_loc = findloc(factors, period, 1)
exponent_loc = scan(s, 'eE', back=.true.)
if(exponent_loc == 0) exponent_loc = len(s) + 1
if(period_loc == 0) period_loc = exponent_loc
where (0 <= factors .and. factors <= 9)
mask = 1
elsewhere
mask = 0
end where
mask_from = 18 - period_loc
mask_till = mask_from + exponent_loc - 2
r_coefficient = sum( &
factors(:exponent_loc - 1) * &
base(mask_from:mask_till) * &
mask(:exponent_loc - 1))
r_exponent = sum( &
factors(exponent_loc+1:len(s)) * &
mask(exponent_loc+1:len(s)) * &
base(17-(len(s)-exponent_loc):16))
if(factors(exponent_loc+1) == minus_sign) r_exponent = -r_exponent
if(factors(1) == minus_sign) r_coefficient = -r_coefficient
r = r_coefficient * 10 ** r_exponent
end function str2real
end module mod_str2real
Test program:
program main
!
use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int64
use mod_st_to_dp
use str2real_m
!
implicit none
!
integer,parameter :: n = 1000000 !! number of values
!
character(len=30),dimension(:),allocatable :: strs
real(dp),dimension(:),allocatable :: rval, rref
!
integer :: i, i0
integer :: ierr
real(dp) :: r
integer(i8) :: start, finish, count_rate
real(dp),parameter :: eps = epsilon(1.0_dp)
!
!
! create a list of values to parse
!
allocate( strs(n), rval(n), rref(n) )
!
! Initializing RNG: with a dummy i
r = util_unirand (i)
i0 = int(n*r) + 1
write(*,*) 'N =', n, ', i0 =', i0
do i = 1, n
rref(i) = util_unirand()
enddo
call system_clock(start, count_rate)
do i = 1, n
write(strs(i), '(E30.16)') rref(i)
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Write: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
blk_io: block
write(*,11) "BLOCK 1: formatted read toward string to double"
!
call system_clock(start, count_rate)
do i = 1, n
read(strs(i),fmt=*,iostat=ierr) rval(i)
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during formatted read."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_io
!! endif
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_io
blk_c: block
write(*,11) "BLOCK 2: C st_to_r8"
!
call system_clock(start, count_rate)
do i = 1, n
call c_st_to_r8 ( strs(i), rval(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during st_to_r8."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_c
!! endif
!
call system_clock(start, count_rate)
call c_st_to_r8 ( strs, rval )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_c
blk_f: block
write(*,11) "BLOCK 3: F st_to_r8"
!
call system_clock (start, count_rate)
do i = 1, n
call f_st_to_r8 ( strs(i), rval(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during st_to_r8."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_f
!! endif
!
call system_clock(start, count_rate)
call f_st_to_r8 ( strs, rval )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_f
blk_4: block
write(*,11) "BLOCK 4: F str2real "
!
call system_clock (start, count_rate)
do i = 1, n
rval(i) = str2real ( strs(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds'
!
call system_clock(start, count_rate)
rval = str2real ( strs )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_4
!
10 format ( 'i =', i9, ', string(i) =', a, ', zval(i) =', e25.17 )
11 format ( /, 80('-'), /, a )
!
deallocate ( strs, rval, rref )
!
STOP
!
contains
!======================================================================
!
! This is the ran0 in Numerical Rrecipes in Fortran.
! To initialize,
! call system_clock( idummy )
! Do not alter idummy between successive deviates.
!
function util_unirand ( idum ) result( ran0 )
!
implicit none
integer,parameter :: rk = kind(1.0D0), ik = 4
!
! ARGUMENTS:
!
real(rk) :: ran0
integer(ik),intent(in),optional :: idum
!
! local constants and variable
!
integer(ik) :: k
integer(ik),save :: idummy
integer(ik),parameter :: &
ia=16807, im=2147483647, iq=127773, ir=2836, mask=123459876
real(rk),parameter :: am = 1.0_rk / im
!
!
if ( present(idum) ) call system_clock( idummy )
!
idummy = ieor( idummy, mask )
k = idummy / iq
idummy = ia*( idummy - k*iq ) - ir*k
if ( idummy .lt. 0 ) idummy = idummy + im
ran0 = am*idummy
idummy = ieor( idummy, mask )
return
end function
!
!=====
subroutine show_err ( er, ep )
real(dp),intent(in) :: er, ep
real(dp) :: tmp
!
tmp = ep - er
write(*,10) er, ep, tmp
if ( tmp .lt. 0.0_dp ) write(*,11) er/ep
return
10 format( 'A = Max |rval(:)-rref(:)| =', 1pe10.3,',', 3x, &
'B = epsilon(1.0d0) =', 1pe10.3,',', 3x, &
'B-A =', 1pe11.3 )
11 format( '*** WARNING: A > B, A/B =', 1pe11.3 )
end subroutine
!
!======================================================================
end program
!======================================================================
Welcome to the Fortran community!
I think you are referring to the very first dummy version of str2real. Have you checked out the updated version from GitHub?
No. I will check your version in GitHub and let you know. Thank you.
To Carltoffel:
The problem to my gfortran compiler is the FINDLOC function in your codes. The compiler does not understand that function. So, I replace
period_loc = findloc(factors, period, 1)
by
do period_loc = 1,N
if ( factors(period_loc) == period ) exit
enddo
Your code works. Here are the results from both Ifort and Gfortran.
1. Gfortran
N = 1000000 , i0 = 378457
Write: time consumed = 0.7423 seconds
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.6479 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1044 seconds
Read: time consumed = 0.1035 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 3: F st_to_r8
Read: time consumed = 0.0964 seconds
Read: time consumed = 0.0979 seconds (array)
A = Max |rval(:)-rref(:)| = 1.887E-15, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.665E-15
*** WARNING: A > B, A/B = 8.500E+00
BLOCK 4: F str2real
Read: time consumed = 0.0580 seconds
Read: time consumed = 0.0603 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00
2. Ifort
N = 1000000 , i0 = 284882
Write: time consumed = 0.3148 seconds
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.4024 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1043 seconds
Read: time consumed = 0.1036 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 3: F st_to_r8
Read: time consumed = 0.0871 seconds
Read: time consumed = 0.0865 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
BLOCK 4: F str2real
Read: time consumed = 0.0589 seconds
Read: time consumed = 0.0584 seconds (array)
A = Max |rval(:)-rref(:)| = 3.331E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.110E-16
*** WARNING: A > B, A/B = 1.500E+00
Regarding your code, result from Gfortran is less accurate than the one from Ifort. However,
your code runs faster (with both compilers) and performs more accurate than f_strtod (with Gfortran). It may depend on compilers and on how we compile codes. In my opinion, all are acceptable. Thank you for sharing your amazing code.
Thank you for your effort to test all these routines!
I think the accuracy is lost because of the vector expression that sums the digits. Maybe I can find a way to increase the accuracy without losing much performance.
Which compiler flags did you use? -Ofast? -march=native?
Compiling with Gfortran, I use many flags:
-O3 -fstack-arrays -static -finit-local-zero -Wall -Wextra -Waliasing -Wconversion -pedantic
The flag “-O3” may be the main role.
With Ifort, I use:
-O3 -ipo -fno-alias -zero -warn unused -warn uncalled -static
Again, the flag “-O3” is responsible.
I haven’t tried your suggestions yet. Thank you.
Hi! I got really inspired by this thread and it helped me boost reading huge mesh files!! since this helped me so much I wanted to share my grain of salt in case it can be useful:
module str2real_m
implicit none
private
public :: str2real
integer, parameter :: wp = kind(1.0d0)
integer, parameter :: ikind = selected_int_kind(2)
integer(kind=ikind), parameter :: digit_0 = ichar('0')
integer(kind=ikind), parameter :: period = ichar('.') - digit_0
integer(kind=ikind), parameter :: minus_sign = ichar('-') - digit_0
real(wp), parameter :: whole_number_base(*) = &
[1d15, 1d14, 1d13, 1d12, 1d11, 1d10, 1d9, 1d8, &
1d7, 1d6, 1d5, 1d4, 1d3, 1d2, 1d1, 1d0]
real(wp), parameter :: fractional_base(*) = &
[1d-1, 1d-2, 1d-3, 1d-4, 1d-5, 1d-6, 1d-7, 1d-8, &
1d-9, 1d-10, 1d-11, 1d-12, 1d-13, 1d-14, 1d-15, 1d-16, &
1d-17, 1d-18, 1d-19, 1d-20, 1d-21, 1d-22, 1d-23, 1d-24]
real(wp), parameter :: period_skip = 0d0
real(wp), parameter :: base(*) = &
[whole_number_base, period_skip, fractional_base]
real(wp), parameter :: expbase(*) = [whole_number_base, fractional_base]
contains
elemental function str2real(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_coefficient
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer(kind=ikind) :: mask(N)
integer :: period_loc
integer :: exponent_loc
integer :: exponent_end
integer :: mask_from
integer :: mask_till
intrinsic :: findloc, scan
integer :: findloc, scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent
equivalence(factors, factors_char)
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
period_loc = findloc(factors, period, 1)
exponent_loc = scan(s, 'eE', back=.true.)
if(exponent_loc == 0) exponent_loc = len(s) + 1
if(period_loc == 0) period_loc = exponent_loc
where (0 <= factors .and. factors <= 9)
mask = 1
elsewhere
mask = 0
end where
if(factors(exponent_loc+1) == minus_sign) sig_exp = -1
if(factors(1)== minus_sign) sig_coef = -1
mask_from = 18 - period_loc
mask_till = mask_from + exponent_loc - 2
r_coefficient = sum( &
factors(:exponent_loc - 1) * &
base(mask_from:mask_till) * &
mask(:exponent_loc - 1))
i_exponent = str2int( s(exponent_loc+1:len(s)) )
r = (sig_coef*r_coefficient) * expbase(16-sig_exp*i_exponent)
end function
elemental function str2int(s) result(int)
character(*), intent(in) :: s
integer :: int
integer :: i, val
int = 0
do i = 1, len(s)
val = iachar(s(i:i))-digit_0
if( val >= 0 .and. val <= 9 ) int = int*10 + val
end do
end function
end module str2real_m
There are a few changes with different purposes:
in order to avoid the problem with white spaces or line endings (CR LF) being mixed in the string I replaced reading the exponent with the str2int function (this was also very useful for reading large arrays of integers).
In order to accelerate a bit more the computation:
2.1) Instead of computing the exponential as a sum I took directly the table of bases and constructed the table expbase(*) = [whole_number_base, fractional_base] which gives the value to multiply the coefficient.
2.2) I stored the signs for the coefficient and exponential in two integer(1) values.
With these changes:
I used the test in the GitHub and added the following check:
call check(" 2.5647869e-003 "//char(13)//char(10))
There is one white space before and after the string plus CRLF at the end! this passes without problem
I also ran the tests by @tqviet ( working in Windows11 / Intel 19 / -O3 )
Before the changes I had for the Block 4: time consumed = 0.086 (both lines)
And after: time consumed = 0.073 for the loop, and 0.07 for the array
Results
--------------------------------------------------------------------------------
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.6290 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
--------------------------------------------------------------------------------
BLOCK 2: C st_to_r8
Read: time consumed = 0.2350 seconds
Read: time consumed = 0.2320 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
--------------------------------------------------------------------------------
BLOCK 3: F st_to_r8
Read: time consumed = 0.1130 seconds
Read: time consumed = 0.1140 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
--------------------------------------------------------------------------------
BLOCK 4: F str2real
Read: time consumed = 0.0730 seconds
Read: time consumed = 0.0700 seconds (array)
A = Max |rval(:)-rref(:)| = 3.331E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.110E-16
*** WARNING: A > B, A/B = 1.500E+00
==============================================================
UPDATE
I tried a “simplification”:
elemental function str2real(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_fraction
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer :: period_loc
integer :: exponent_loc
intrinsic :: findloc, scan
integer :: findloc, scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent, i, val
equivalence(factors, factors_char)
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
period_loc = findloc(factors, period, 1)
exponent_loc = scan(s, 'eE', back=.true.)
if(exponent_loc == 0) exponent_loc = len(s) + 1
if(period_loc == 0) period_loc = exponent_loc
if(factors(exponent_loc+1) == minus_sign) sig_exp = -1
if(factors(1)== minus_sign) sig_coef = -1
val = str2int( s(1:period_loc-1) )
r = val
do i = period_loc+1, exponent_loc-1
val = factors(i)
if( val >= 0 .and. val <= 9 ) r = r + val*fractional_base(i-period_loc)
end do
r = sig_coef*r
i_exponent = str2int( s(exponent_loc+1:len(s)) )
if(i_exponent.ne.0) r = r * expbase(16-sig_exp*i_exponent)
end function
By computing the whole number part with the str2int function and then iterating on the fractional part plus avoiding multiplying by the exponential coefficient if not needed it brough my times down from 0.07s to 0.052s
Results
--------------------------------------------------------------------------------
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.6270 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
--------------------------------------------------------------------------------
BLOCK 2: C st_to_r8
Read: time consumed = 0.2370 seconds
Read: time consumed = 0.2320 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
--------------------------------------------------------------------------------
BLOCK 3: F st_to_r8
Read: time consumed = 0.1390 seconds
Read: time consumed = 0.1360 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
--------------------------------------------------------------------------------
BLOCK 4: F str2real
Read: time consumed = 0.0520 seconds
Read: time consumed = 0.0520 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00
==============================================================
UPDATE 2
Just the change w.r.t to previous update of the last lines of the str2real function:
val = str2int( s(1:period_loc-1) )
r = sum( factors(period_loc+1:exponent_loc-1)*fractional_base(1:exponent_loc-period_loc-1) )
r = sig_coef*(r+val)
i_exponent = str2int( s(exponent_loc+1:len(s)) )
if(i_exponent.ne.0) r = r * expbase(16-sig_exp*i_exponent)
Results
--------------------------------------------------------------------------------
BLOCK 1: formatted read toward string to double
Read: time consumed = 0.6390 seconds
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
--------------------------------------------------------------------------------
BLOCK 2: C st_to_r8
Read: time consumed = 0.2370 seconds
Read: time consumed = 0.2350 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
--------------------------------------------------------------------------------
BLOCK 3: F st_to_r8
Read: time consumed = 0.1160 seconds
Read: time consumed = 0.1140 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
--------------------------------------------------------------------------------
BLOCK 4: F str2real
Read: time consumed = 0.0480 seconds
Read: time consumed = 0.0480 seconds (array)
A = Max |rval(:)-rref(:)| = 2.220E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = 0.000E+00
My time went down now to 0.048 and the error is now B-A=0
By the way @Carltoffel I think in line 41 of your test.f90 file the test should be:
if(abs(rel_err) > epsilon(1.0d0)) then
Otherwise I was getting a print with call check(“0.1234E0”), when both values were identical to the last decimal but “abs_err = str2real_out - formatted_read_out” was returning a value smaller than the precision
Thank you @ hkvzjal for your sharing. After checking your codes and comparing them with the previous versions, I conclude that your version with UPDATE 2 runs fastest, while your other versions also show very good performance. It is very interesting!
Herein I’d like to post all modules and the testing program, as follows:
module mod_str2real_02
!
implicit none
private
!
public :: str2real_02
public :: str2real_03
public :: str2real_04
!
integer, parameter :: wp = kind(1.0d0)
integer, parameter :: ikind = selected_int_kind(2)
integer(kind=ikind), parameter :: digit_0 = ichar('0')
integer(kind=ikind), parameter :: period = ichar('.')-digit_0
integer(kind=ikind), parameter :: minus_sign = ichar('-')-digit_0
!
real(wp), parameter :: whole_number_base(*) = &
[1d15, 1d14, 1d13, 1d12, 1d11, 1d10, 1d9, 1d8, &
1d7, 1d6, 1d5, 1d4, 1d3, 1d2, 1d1, 1d0]
!
real(wp), parameter :: fractional_base(*) = &
[1d-1, 1d-2, 1d-3, 1d-4, 1d-5, 1d-6, 1d-7, 1d-8, &
1d-9, 1d-10, 1d-11, 1d-12, 1d-13, 1d-14, 1d-15, 1d-16, &
1d-17, 1d-18, 1d-19, 1d-20, 1d-21, 1d-22, 1d-23, 1d-24]
!
real(wp), parameter :: period_skip = 0d0
real(wp), parameter :: base(*) = &
[whole_number_base, period_skip, fractional_base]
real(wp), parameter :: expbase(*) = &
[whole_number_base, fractional_base]
!
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
! hkvzjal - version 01
!
elemental &
function str2real_02(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_coefficient
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer(kind=ikind) :: mask(N)
integer :: period_loc
integer :: exponent_loc
integer :: exponent_end
integer :: mask_from
integer :: mask_till
!! intrinsic :: findloc
!! integer :: findloc
intrinsic :: scan
integer :: scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent
!
!
equivalence(factors, factors_char)
!
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
do period_loc = 1,N
if ( factors(period_loc) == period ) exit
enddo
!
exponent_loc = scan(s, 'eE', back=.true.)
!
if(exponent_loc == 0) exponent_loc = len(s) + 1
if(period_loc == 0) period_loc = exponent_loc
!
where (0 <= factors .and. factors <= 9)
mask = 1
elsewhere
mask = 0
end where
!
if(factors(exponent_loc+1) == minus_sign) sig_exp = -1
if(factors(1)== minus_sign) sig_coef = -1
!
mask_from = 18 - period_loc
mask_till = mask_from + exponent_loc - 2
!
r_coefficient = sum( &
factors(:exponent_loc - 1) * &
base(mask_from:mask_till) * &
mask(:exponent_loc - 1))
!
i_exponent = str2int( s(exponent_loc+1:len(s)) )
r = (sig_coef*r_coefficient) * expbase(16-sig_exp*i_exponent)
!
end function
!------------------------------------------------------------------------------
elemental &
function str2int(s) result(iou)
character(*), intent(in) :: s
integer :: iou
integer :: i, val
iou = 0
do i = 1, len(s)
val = iachar(s(i:i))-digit_0
if( val >= 0 .and. val <= 9 ) iou = iou*10 + val
end do
end function
!------------------------------------------------------------------------------
! hkvzjal - version 02
!
elemental &
function str2real_03(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_fraction
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer :: period_loc
integer :: exponent_loc
!! intrinsic :: findloc
!! integer :: findloc
intrinsic :: scan
integer :: scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent, i, val
!
equivalence(factors, factors_char)
!
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
do period_loc = 1,N
if ( factors(period_loc) == period ) exit
enddo
!
exponent_loc = scan(s, 'eE', back=.true.)
if (exponent_loc == 0) exponent_loc = len(s) + 1
if (period_loc == 0) period_loc = exponent_loc
if (factors(exponent_loc+1) == minus_sign) sig_exp = -1
if (factors(1)== minus_sign) sig_coef = -1
val = str2int( s(1:period_loc-1) )
r = val
do i = period_loc+1, exponent_loc-1
val = factors(i)
if( val >= 0 .and. val <= 9 ) &
r = r + val*fractional_base(i-period_loc)
end do
r = sig_coef*r
i_exponent = str2int( s(exponent_loc+1:len(s)) )
if(i_exponent.ne.0) r = r * expbase(16-sig_exp*i_exponent)
!
end function
!------------------------------------------------------------------------------
! hkvzjal - version 03
!
elemental &
function str2real_04 (s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_fraction
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer :: period_loc
integer :: exponent_loc
!! intrinsic :: findloc
!! integer :: findloc
intrinsic :: scan
integer :: scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent, i, val
!
equivalence(factors, factors_char)
!
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
do period_loc = 1,N
if ( factors(period_loc) == period ) exit
enddo
!
exponent_loc = scan(s, 'eE', back=.true.)
if (exponent_loc == 0) exponent_loc = len(s) + 1
if (period_loc == 0) period_loc = exponent_loc
if (factors(exponent_loc+1) == minus_sign) sig_exp = -1
if (factors(1)== minus_sign) sig_coef = -1
val = str2int( s(1:period_loc-1) )
r = sum( factors(period_loc+1:exponent_loc-1) * &
fractional_base(1:exponent_loc-period_loc-1) )
r = sig_coef*(r+val)
i_exponent = str2int( s(exponent_loc+1:len(s)) )
if(i_exponent.ne.0) r = r * expbase(16-sig_exp*i_exponent)
!
end function
!------------------------------------------------------------------------------
end module mod_str2real_02
!------------------------------------------------------------------------------
module mod_str2real
!
implicit none
private
!
public :: str2real
!
integer, parameter :: wp = kind(1.0d0)
integer, parameter :: ikind = selected_int_kind(2)
!
integer(kind=ikind), parameter :: digit_0 = ichar('0')
integer(kind=ikind), parameter :: period = ichar('.') - digit_0
integer(kind=ikind), parameter :: minus_sign = ichar('-') - digit_0
!
real(wp), parameter :: whole_number_base(*) = &
[1d15, 1d14, 1d13, 1d12, 1d11, 1d10, 1d9, 1d8, &
1d7, 1d6, 1d5, 1d4, 1d3, 1d2, 1d1, 1d0]
real(wp), parameter :: period_skip = 0d0
real(wp), parameter :: fractional_base(*) = &
[1d-1, 1d-2, 1d-3, 1d-4, 1d-5, 1d-6, 1d-7, 1d-8, &
1d-9, 1d-10, 1d-11, 1d-12, 1d-13, 1d-14, 1d-15, 1d-16, &
1d-17, 1d-18, 1d-19, 1d-20, 1d-21, 1d-22, 1d-23, 1d-24]
real(wp), parameter :: base(*) = &
[whole_number_base, period_skip, fractional_base]
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
elemental &
function str2real(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
!
real(wp) :: r_coefficient
real(wp) :: r_exponent
!
! integer, parameter :: r10 = selected_real_kind ( P=16 ) !<-- not more accurate
! real(r10) :: r_coefficient
! real(r10) :: r_exponent
!
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer(kind=ikind) :: mask(N)
integer :: period_loc
integer :: exponent_loc
integer :: mask_from
integer :: mask_till
! intrinsic :: findloc
! integer :: findloc
intrinsic :: scan
integer :: scan
! integer :: i
equivalence(factors, factors_char)
factors_char = s
factors = factors - digit_0
!! period_loc = findloc(factors, period, 1)
do period_loc = 1,N
if ( factors(period_loc) == period ) exit
enddo
exponent_loc = scan(s, 'eE', back=.true.)
if(exponent_loc == 0) exponent_loc = len(s) + 1
if(period_loc == 0) period_loc = exponent_loc
where (0 <= factors .and. factors <= 9)
mask = 1
elsewhere
mask = 0
end where
!
mask_from = 18 - period_loc
mask_till = mask_from + exponent_loc - 2
!
r_coefficient = sum( &
factors(:exponent_loc - 1) * &
base(mask_from:mask_till) * &
mask(:exponent_loc - 1) )
r_exponent = sum( &
factors(exponent_loc+1:len(s)) * &
mask(exponent_loc+1:len(s)) * &
base(17-(len(s)-exponent_loc):16) )
if ( factors(exponent_loc+1) == minus_sign ) r_exponent = -r_exponent
if ( factors(1) == minus_sign ) r_coefficient = -r_coefficient
r = r_coefficient * (10**r_exponent)
!
end function str2real
!------------------------------------------------------------------------------
end module mod_str2real
and
module mod_st_to_dp
!
! Fortran 2003:
use, intrinsic :: iso_c_binding, only: c_double, c_char, c_ptr, &
c_null_ptr, c_long
!
! Fortran 2008:
use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int64
!
implicit none
!
private
!
public :: f_st_to_r8
public :: c_st_to_r8
!
!
interface
pure &
function c_strtod ( str, endptr ) result(d) bind(C, name="strtod")
!
!! <stdlib.h> :: double strtod(const char *str, char **endptr)
!
import :: c_char, c_ptr, c_double
!! Argument list
character(kind=c_char,len=1), intent(in) :: str(*)
!! type(c_ptr), intent(inout) :: endptr !<- since we need pure attribution.
type(c_ptr), intent(in) :: endptr
!! function result
real(c_double) :: d
end function
end interface
!
!
! integer, parameter :: dp = selected_real_kind( p=12 )
! FP constants
real(dp), parameter :: ten = 10.0_dp
! Integer constants
integer, parameter :: ascii_negative = 45
integer, parameter :: ascii_period = 46
integer, parameter :: ascii_0 = 48
integer, parameter :: ascii_9 = 57
integer(i8), parameter :: ten_int = 10_i8
!
type(c_ptr),parameter :: endptr = c_null_ptr
!
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
elemental &
subroutine c_st_to_r8 ( str, r )
character(len=*), intent(in) :: str
real(dp), intent(out) :: r
!
r = c_strtod ( str, endptr )
return
end subroutine
!------------------------------------------------------------------------------
elemental &
subroutine f_st_to_r8 ( str, r )
character(len=*), intent(in) :: str
real(dp), intent(out) :: r
! local variables
integer(i8) :: n, n_exp
integer :: lens, pos_exp, expnt
!
r = 0.0_dp
n_exp = 0
lens = len_trim( str )
pos_exp = index( str, "E", back=.true. ) ! <-- find the capital 'E'
if ( pos_exp == 0 ) then
pos_exp = index( str, "e", back=.true. )
endif
if ( pos_exp > 0 ) then
call str2dec( str(pos_exp+1:lens), n_exp )
else
pos_exp = lens + 1
endif
call str2dec( str(1:pos_exp-1), n, expnt )
r = real( n, kind=dp ) * ( ten**(n_exp + expnt) )
return
end subroutine
!------------------------------------------------------------------------------
elemental &
subroutine str2dec( s, n, expnt )
character(len=*), intent(in) :: s
integer(i8), intent(out) :: n
integer, intent(out), optional :: expnt
integer :: ipos, ic, pos_period, iexponent
!
n = 0_i8
pos_period = 0
iexponent = -1
if ( present(expnt) ) iexponent = 0
do ipos = len(s), 1, -1
ic = ichar( s(ipos:ipos) )
if ( present(expnt) ) then
if ( ic == ascii_period ) then
pos_period = ipos
cycle
endif
endif
if ( ic == ascii_negative ) then
n = -n
exit
endif
if ( (ic < ascii_0) .or. (ic > ascii_9) ) exit
iexponent = iexponent + 1
n = n + (ten_int**iexponent)*(ic - 48)
enddo
if ( present(expnt) ) expnt = -len(s) + pos_period - 1
return
end subroutine
!------------------------------------------------------------------------------
end module
!------------------------------------------------------------------------------
program main
!
use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int64
use mod_st_to_dp
use mod_str2real
use mod_str2real_02
!
implicit none
!
integer,parameter :: n = 1000000 !! number of values
!
character(len=30),dimension(:),allocatable :: strs
real(dp),dimension(:),allocatable :: rval, rref
!
integer :: i, i0
integer :: ierr
real(dp) :: r
integer(i8) :: start, finish, count_rate
real(dp),parameter :: eps = epsilon(1.0_dp)
real(dp) :: dt(20)
!
!
! create a list of values to parse
!
allocate( strs(n), rval(n), rref(n) )
!
! Initializing RNG: with a dummy i
r = util_unirand (i)
i0 = int(n*r) + 1
write(*,*) 'N =', n, ', i0 =', i0
do i = 1, n
rref(i) = util_unirand()
enddo
call system_clock(start, count_rate)
do i = 1, n
write(strs(i), '(E30.16)') rref(i)
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Write: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
!
blk_io: block
!
write(*,11) "BLOCK 1: Direct reading formatted string to double"
!
call system_clock(start, count_rate)
do i = 1, n
read(strs(i),fmt=*,iostat=ierr) rval(i)
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(1) = (finish-start)/real(count_rate,dp)
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during formatted read."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_io
!! endif
!
call show_err ( maxval(abs(rval-rref)), eps )
!
end block blk_io
!
!
blk_c: block
!
write(*,11) "BLOCK 2: C st_to_r8"
!
call system_clock(start, count_rate)
do i = 1, n
call c_st_to_r8 ( strs(i), rval(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(2) = (finish-start)/real(count_rate,dp)
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during st_to_r8."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_c
!! endif
!
call system_clock(start, count_rate)
call c_st_to_r8 ( strs, rval )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
write(*,'(/,"Time consumed: Direct / This (serial) = ", F7.4)') dt(1) / dt(2)
!
end block blk_c
!
!
blk_f: block
!
write(*,11) "BLOCK 3: F st_to_r8"
!
call system_clock (start, count_rate)
do i = 1, n
call f_st_to_r8 ( strs(i), rval(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(3) = (finish-start)/real(count_rate,dp)
!
!! if ( abs(rval(i0)-rref(i0)) > eps ) then
!! print *, "Warning: mismatch during st_to_r8."
!! write(*,10) i0, trim(strs(i0)), rval(i0)
! exit blk_f
!! endif
!
call system_clock(start, count_rate)
call f_st_to_r8 ( strs, rval )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
write(*,'(/,"Time consumed: Direct / This (serial) = ", F7.4)') dt(1) / dt(3)
!
end block blk_f
!
!
blk_4: block
!
write(*,11) "BLOCK 4: F str2real (Carltoffel)"
!
call system_clock (start, count_rate)
do i = 1, n
rval(i) = str2real ( strs(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(4) = (finish-start)/real(count_rate,dp)
!
call system_clock(start, count_rate)
rval = str2real ( strs )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
write(*,'(/,"Time consumed: Direct / This (serial) = ", F7.4)') dt(1) / dt(4)
!
end block blk_4
!
!
blk_5: block
!
write(*,11) "BLOCK 5: F str2real_02 (hkvzjal - version 01)"
!
call system_clock (start, count_rate)
do i = 1, n
rval(i) = str2real_02 ( strs(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(5) = (finish-start)/real(count_rate,dp)
!
call system_clock(start, count_rate)
rval = str2real_02 ( strs )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
write(*,'(/,"Time consumed: Direct / This (serial) = ", F7.4)') dt(1) / dt(5)
!
end block blk_5
!
!
blk_6: block
!
write(*,11) "BLOCK 6: F str2real_03 (hkvzjal - version 02) "
!
call system_clock (start, count_rate)
do i = 1, n
rval(i) = str2real_03 ( strs(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(6) = (finish-start)/real(count_rate,dp)
!
call system_clock(start, count_rate)
rval = str2real_03 ( strs )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
write(*,'(/,"Time consumed: Direct / This (serial) = ", F7.4)') dt(1) / dt(6)
!
end block blk_6
!
!
blk_7: block
!
write(*,11) "BLOCK 7: F str2real_04 (hkvzjal - version 03) "
!
call system_clock (start, count_rate)
do i = 1, n
rval(i) = str2real_04 ( strs(i) )
enddo
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (serial)'
!
dt(7) = (finish-start)/real(count_rate,dp)
!
call system_clock(start, count_rate)
rval = str2real_04 ( strs )
call system_clock(finish)
write(*,'(A30,1X,F7.4,1X,A)') 'Read: time consumed =', &
(finish-start)/real(count_rate,dp), ' seconds (array)'
!
call show_err ( maxval(abs(rval-rref)), eps )
!
write(*,'(/,"Time consumed: Direct / This (serial) = ", F7.4)') dt(1) / dt(7)
!
end block blk_7
!
10 format ( 'i =', i9, ', string(i) =', a, ', zval(i) =', e25.17 )
11 format ( /, 80('-'), /, a )
!
deallocate ( strs, rval, rref )
!
STOP
!
contains
!======================================================================
!
! This is the ran0 in Numerical Rrecipes in Fortran.
! To initialize,
! call system_clock( idummy )
! Do not alter idummy between successive deviates.
!
function util_unirand ( idum ) result( ran0 )
!
implicit none
integer,parameter :: rk = kind(1.0D0), ik = 4
!
! ARGUMENTS:
!
real(rk) :: ran0
integer(ik),intent(in),optional :: idum
!
! local constants and variable
!
integer(ik) :: k
integer(ik),save :: idummy
integer(ik),parameter :: &
ia=16807, im=2147483647, iq=127773, ir=2836, mask=123459876
real(rk),parameter :: am = 1.0_rk / im
!
!
if ( present(idum) ) call system_clock( idummy )
!
idummy = ieor( idummy, mask )
k = idummy / iq
idummy = ia*( idummy - k*iq ) - ir*k
if ( idummy .lt. 0 ) idummy = idummy + im
ran0 = am*idummy
idummy = ieor( idummy, mask )
return
end function
!
!=====
subroutine show_err ( er, ep )
real(dp),intent(in) :: er, ep
real(dp) :: tmp
!
tmp = ep - er
write(*,10) er, ep, tmp
if ( tmp .lt. 0.0_dp ) write(*,11) er/ep
return
10 format( 'A = Max |rval(:)-rref(:)| =', 1pe10.3,',', 3x, &
'B = epsilon(1.0d0) =', 1pe10.3,',', 3x, &
'B-A =', 1pe11.3 )
11 format( '*** WARNING: A > B, A/B =', 1pe11.3 )
end subroutine
!======================================================================
end program
N = 1000000 , i0 = 847605
Write: time consumed = 0.7626 seconds (serial)
BLOCK 1: Direct reading formatted string to double
Read: time consumed = 0.6423 seconds (serial)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1039 seconds (serial)
Read: time consumed = 0.1060 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16Time consumed: Direct / This (serial) = 6.1831
BLOCK 3: F st_to_r8
Read: time consumed = 0.0922 seconds (serial)
Read: time consumed = 0.0896 seconds (array)
A = Max |rval(:)-rref(:)| = 1.887E-15, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.665E-15
*** WARNING: A > B, A/B = 8.500E+00Time consumed: Direct / This (serial) = 6.9675
BLOCK 4: F str2real (Carltoffel)
Read: time consumed = 0.0525 seconds (serial)
Read: time consumed = 0.0501 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 12.2401
BLOCK 5: F str2real_02 (hkvzjal - version 01)
Read: time consumed = 0.0407 seconds (serial)
Read: time consumed = 0.0408 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 15.7848
BLOCK 6: F str2real_03 (hkvzjal - version 02)
Read: time consumed = 0.0405 seconds (serial)
Read: time consumed = 0.0410 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 15.8451
BLOCK 7: F str2real_04 (hkvzjal - version 03)
Read: time consumed = 0.0370 seconds (serial)
Read: time consumed = 0.0381 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 17.3503
Here, “Time consumed: Direct / This (serial) = 17.3503” means 0.6423/0.0370 = 17.35xxx. It shows the speed of the method with UPDATE 2 is fastest!
Update: Since I replaced the findloc intrinsic function (which is not available in my old compiler), there may be a bug in the following line:
if ( period_loc == 0 ) period_loc = exponent_loc
Hence, I use a logical variable (i.e. not_found) to check if the findloc does not find “period”. Herein I also include all versions in only one module, i.e. remove the module mod_str2real_02, therefore, we should remove the line
use mod_str2real_02
from the main program above. The code runs a bit faster. (But I’m not sure if this change make such good effect.)
Here is the new version:
!------------------------------------------------------------------------------
module mod_str2real
!
implicit none
private
!
public :: str2real
public :: str2real_02
public :: str2real_03
public :: str2real_04
!
integer, parameter :: wp = kind(1.0d0)
integer, parameter :: ikind = selected_int_kind(2)
!
integer(kind=ikind), parameter :: digit_0 = ichar('0')
integer(kind=ikind), parameter :: period = ichar('.') - digit_0
integer(kind=ikind), parameter :: minus_sign = ichar('-') - digit_0
!
real(wp), parameter :: whole_number_base(*) = &
[1d15, 1d14, 1d13, 1d12, 1d11, 1d10, 1d9, 1d8, &
1d7, 1d6, 1d5, 1d4, 1d3, 1d2, 1d1, 1d0]
real(wp), parameter :: period_skip = 0d0
real(wp), parameter :: fractional_base(*) = &
[1d-1, 1d-2, 1d-3, 1d-4, 1d-5, 1d-6, 1d-7, 1d-8, &
1d-9, 1d-10, 1d-11, 1d-12, 1d-13, 1d-14, 1d-15, 1d-16, &
1d-17, 1d-18, 1d-19, 1d-20, 1d-21, 1d-22, 1d-23, 1d-24]
real(wp), parameter :: base(*) = &
[ whole_number_base, period_skip, fractional_base ]
!
real(wp), parameter :: expbase(*) = &
[ whole_number_base, fractional_base ]
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
elemental &
function str2real (s) result(r)
character(*), intent(in) :: s
real(wp) :: r
!
real(wp) :: r_coefficient
real(wp) :: r_exponent
!
! integer, parameter :: r10 = selected_real_kind ( P=16 ) !<-- not more accurate
! real(r10) :: r_coefficient
! real(r10) :: r_exponent
!
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer(kind=ikind) :: mask(N)
integer :: period_loc
integer :: exponent_loc
integer :: mask_from
integer :: mask_till
! intrinsic :: findloc
! integer :: findloc
intrinsic :: scan
integer :: scan
logical :: not_found
!
equivalence(factors, factors_char)
!
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
!! if ( period_loc == 0 ) period_loc = exponent_loc
!
not_found = .true.
do period_loc = 1,N
if ( factors(period_loc) == period ) then
not_found = .false.
exit
endif
enddo
exponent_loc = scan(s, 'eE', back=.true.)
if ( exponent_loc == 0 ) exponent_loc = len(s) + 1
if ( not_found ) period_loc = exponent_loc
!
where (0 <= factors .and. factors <= 9)
mask = 1
elsewhere
mask = 0
end where
!
mask_from = 18 - period_loc
mask_till = mask_from + exponent_loc - 2
!
r_coefficient = sum( &
factors(:exponent_loc - 1) * &
base(mask_from:mask_till) * &
mask(:exponent_loc - 1) )
r_exponent = sum( &
factors(exponent_loc+1:len(s)) * &
mask(exponent_loc+1:len(s)) * &
base(17-(len(s)-exponent_loc):16) )
if ( factors(exponent_loc+1) == minus_sign ) r_exponent = -r_exponent
if ( factors(1) == minus_sign ) r_coefficient = -r_coefficient
r = r_coefficient * (10**r_exponent)
return
end function str2real
!------------------------------------------------------------------------------
! hkvzjal - version 01
!
elemental &
function str2real_02(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_coefficient
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer(kind=ikind) :: mask(N)
integer :: period_loc
integer :: exponent_loc
integer :: exponent_end
integer :: mask_from
integer :: mask_till
!! intrinsic :: findloc
!! integer :: findloc
intrinsic :: scan
integer :: scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent
logical :: not_found
!
equivalence(factors, factors_char)
!
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
!! if ( period_loc == 0 ) period_loc = exponent_loc
!
not_found = .true.
do period_loc = 1,N
if ( factors(period_loc) == period ) then
not_found = .false.
exit
endif
enddo
exponent_loc = scan( s, 'eE', back=.true. )
if ( exponent_loc == 0 ) exponent_loc = len(s) + 1
if ( not_found ) period_loc = exponent_loc
!
where (0 <= factors .and. factors <= 9)
mask = 1
elsewhere
mask = 0
end where
!
if ( factors(exponent_loc+1) == minus_sign ) sig_exp = -1
if ( factors(1) == minus_sign ) sig_coef = -1
!
mask_from = 18 - period_loc
mask_till = mask_from + exponent_loc - 2
!
r_coefficient = sum( &
factors(:exponent_loc - 1) * &
base(mask_from:mask_till) * &
mask(:exponent_loc - 1))
!
i_exponent = str2int( s(exponent_loc+1:len(s)) )
r = (sig_coef*r_coefficient) * expbase(16-sig_exp*i_exponent)
return
end function
!------------------------------------------------------------------------------
elemental &
function str2int(s) result(iou)
character(*), intent(in) :: s
integer :: iou
integer :: i, val
iou = 0
do i = 1, len(s)
val = iachar(s(i:i))-digit_0
if( val >= 0 .and. val <= 9 ) iou = iou*10 + val
end do
return
end function
!------------------------------------------------------------------------------
! hkvzjal - version 02
!
elemental &
function str2real_03(s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_fraction
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer :: period_loc
integer :: exponent_loc
!! intrinsic :: findloc
!! integer :: findloc
intrinsic :: scan
integer :: scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent, i, val
logical :: not_found
!
equivalence(factors, factors_char)
!
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
!! if (period_loc == 0) period_loc = exponent_loc
!
not_found = .true.
do period_loc = 1,N
if ( factors(period_loc) == period ) then
not_found = .false.
exit
endif
enddo
exponent_loc = scan( s, 'eE', back=.true. )
if ( exponent_loc == 0 ) exponent_loc = len(s) + 1
if ( not_found ) period_loc = exponent_loc
!
if ( factors(exponent_loc+1) == minus_sign ) sig_exp = -1
if ( factors(1) == minus_sign ) sig_coef = -1
val = str2int( s(1:period_loc-1) )
r = val
do i = period_loc+1, exponent_loc-1
val = factors(i)
if( val >= 0 .and. val <= 9 ) &
r = r + val*fractional_base(i-period_loc)
end do
r = sig_coef*r
i_exponent = str2int( s(exponent_loc+1:len(s)) )
if(i_exponent.ne.0) r = r * expbase(16-sig_exp*i_exponent)
return
end function
!------------------------------------------------------------------------------
! hkvzjal - version 03 (the fastest one)
!
elemental &
function str2real_04 (s) result(r)
character(*), intent(in) :: s
real(wp) :: r
real(wp) :: r_fraction
integer, parameter :: N = 32
character(N) :: factors_char
integer(kind=ikind) :: factors(N)
integer :: period_loc
integer :: exponent_loc
!! intrinsic :: findloc
!! integer :: findloc
intrinsic :: scan
integer :: scan
integer(1) :: sig_coef, sig_exp
integer :: i_exponent, i, val
logical :: not_found
!
equivalence(factors, factors_char)
!
sig_coef = 1; sig_exp = 1
factors_char = s
factors = factors - digit_0
!
!! period_loc = findloc(factors, period, 1)
!! if (period_loc == 0) period_loc = exponent_loc
!
not_found = .true.
do period_loc = 1,N
if ( factors(period_loc) == period ) then
not_found = .false.
exit
endif
enddo
exponent_loc = scan( s, 'eE', back=.true. )
if ( exponent_loc == 0 ) exponent_loc = len(s) + 1
if ( not_found ) period_loc = exponent_loc
!
if ( factors(exponent_loc+1) == minus_sign ) sig_exp = -1
if ( factors(1) == minus_sign ) sig_coef = -1
val = str2int( s(1:period_loc-1) )
r = sum( factors(period_loc+1:exponent_loc-1) * &
fractional_base(1:exponent_loc-period_loc-1) )
r = sig_coef*(r+val)
i_exponent = str2int( s(exponent_loc+1:len(s)) )
if(i_exponent.ne.0) r = r * expbase(16-sig_exp*i_exponent)
return
end function
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
end module mod_str2real
!------------------------------------------------------------------------------
and new result:
N = 1000000 , i0 = 55962
Write: time consumed = 0.7679 seconds (serial)
BLOCK 1: Direct reading formatted string to double
Read: time consumed = 0.6444 seconds (serial)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16
BLOCK 2: C st_to_r8
Read: time consumed = 0.1050 seconds (serial)
Read: time consumed = 0.1047 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-17, B = epsilon(1.0d0) = 2.220E-16, B-A = 1.665E-16Time consumed: Direct / This (serial) = 6.1357
BLOCK 3: F st_to_r8
Read: time consumed = 0.0910 seconds (serial)
Read: time consumed = 0.0915 seconds (array)
A = Max |rval(:)-rref(:)| = 1.887E-15, B = epsilon(1.0d0) = 2.220E-16, B-A = -1.665E-15
*** WARNING: A > B, A/B = 8.500E+00Time consumed: Direct / This (serial) = 7.0826
BLOCK 4: F str2real (Carltoffel)
Read: time consumed = 0.0576 seconds (serial)
Read: time consumed = 0.0589 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 11.1965
BLOCK 5: F str2real_02 (hkvzjal - version 01)
Read: time consumed = 0.0400 seconds (serial)
Read: time consumed = 0.0407 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 16.1034
BLOCK 6: F str2real_03 (hkvzjal - version 02)
Read: time consumed = 0.0375 seconds (serial)
Read: time consumed = 0.0381 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 17.1991
BLOCK 7: F str2real_04 (hkvzjal - version 03)
Read: time consumed = 0.0346 seconds (serial)
Read: time consumed = 0.0361 seconds (array)
A = Max |rval(:)-rref(:)| = 5.551E-16, B = epsilon(1.0d0) = 2.220E-16, B-A = -3.331E-16
*** WARNING: A > B, A/B = 2.500E+00Time consumed: Direct / This (serial) = 18.6250
@hkvzjal and @Carltoffel ,
A few weeks ago I made an attempt at porting the C++ fast_float library to Fortran (which you can find here). My goal is to have the interface support reporting an invalid number. However, without the unsigned integers and defined bit representations, the exact algorithm from the C++ library would not be portable Fortran. I initially still tried to follow their algorithm in spirit and found it didn’t end up always being faster than list-directed read. I then tried something (I think) similar to what you’ve done in str2real. I still found it wasn’t always faster. Would either of you be willing to try using your algorithm with my interface and test suite?
Hey @tqviet , thanks a lot for putting everything together!!
Just for information, I ran the tests 6 times with the findloc and with your boolean version
Compiler: Intel 19.1.3.311 (x64 architecture)
CPU: Intel(R) Xeon(R) Gold 6226R CPU @ 2.90GHz 2.89 GHz
(results for the last version)
With findloc:
With boolean:
All runs with the intrinsic gave me B-A =0,
Just 1 of the runs of the boolean implementation gave me B-A=-1.1106E-16 (on rref(i)=0.678125232773891 → rval(i) =0.678125232773892 )
@everythingfunctional I’ll take a look and try to come back soon