Hey guys, I’ve been working with Fortran about 6 years now but only started playing with it for some personal projects recently. Today I came across a very strange bug, which seems to only affect gfortran. I have a simple math module which implements sum, average, and variance calculations as below:
math module with assumed rank functions
module math
implicit none
private
public :: sum_reduce
public :: avg
public :: var
contains
function var(x) result(y)
implicit none
real, intent(in) :: x(..)
real :: y
select rank (x)
rank (0)
y = 0.0
rank (1)
write(*,*) 'size(inner): ',size(x - avg(x))
y = sum_reduce((x - avg(x))**2)/real(size(x))
rank default
error stop '- VAR only defined for scalar or vector <x>'
end select
end function var
pure function avg(x) result(y)
implicit none
real, intent(in) :: x(..)
real :: y
select rank (x)
rank (0)
y = x
rank (1)
y = sum_reduce(x)/real(size(x))
rank default
error stop '- AVG only defined for scalar or vector <x>'
end select
end function avg
pure recursive function sum_reduce(x) result(y)
implicit none
real, intent(in) :: x(..)
real :: y
integer :: n
select rank (x)
rank (0)
y = x
rank (1)
n = size(x)
if (n.gt.2) then
y = sum_reduce(x(1:(n/2))) + sum_reduce(x(((n/2)+1):n))
else
y = sum(x)
end if
rank default
error stop '- SUM_REDUCE only defined for scalar or vector <x>'
end select
end function sum_reduce
end module math
My test program is as follows:
test program
program main
use math
implicit none
integer, parameter :: n = 2000000
real :: x(n), x_sum, x_avg, x_var
call random_number(x)
write(*,'(a,i0,2(a,e13.6))') 'x has ',size(x),' elements, minval: ',minval(x),', maxval: ',maxval(x)
x_sum = sum_reduce(x)
write(*,*) 'x_sum: ',x_sum
x_avg = avg(x)
write(*,*) 'x_avg: ',x_avg
x_var = var(x)
write(*,*) 'x_var: ',x_var
end program main
Which yields output:
assumed rank variance calculation fails with array out of bounds
fpm: Entering directory '/home/tyranids/projects/math-error'
Project is up to date
x has 2000000 elements, minval: 0.119209E-06, maxval: 0.100000E+01
x_sum: 1000636.81
x_avg: 0.500318408
size(inner): 2000000
At line 51 of file app/math.f90
Fortran runtime error: Index '2000000' of dimension 1 of array '__tmp_REAL_4_rank_1' outside of expected range (1999999:0)
Error termination. Backtrace:
#0 0x14d0046b4ad0 in ???
#1 0x14d0046b5649 in ???
#2 0x14d0046b5c46 in ???
#3 0x5575e14e7c7f in __math_MOD_sum_reduce
at app/math.f90:51
#4 0x5575e14e8433 in __math_MOD_var
at app/math.f90:20
#5 0x5575e14e7739 in MAIN__
at app/main.f90:18
#6 0x5575e14e77ff in main
at app/main.f90:2
<ERROR> Execution failed for object " math-error "
<ERROR>*cmd_run*:stopping due to failed executions
STOP 1
fpm appears to run with array bounds checking enabled by default, which is what caught this. I can confirm that with ifort, even with -check all, this does not happen. Why is gfortran indexing the temporary array it creates starting from 0? From the debug print added in var, I can see that the size of the arguments provided to sum_reduce is the expected 2000000. I can also see that sum_reduce works correctly for both the direct call and the one from avg, both print the expected value earlier than the call to var, and neither trips the array out of bounds check.
After additional investigation, I have copied the math module to a new math2 version, which uses all assumed shape arrays of dimension 1, since that is my currently failing use case. With this new code, I get no errors and it works as expected. Is there a known problem with chaining assumed rank procedures together in gfortran, but not Intel’s compilers?
math module with assumed shape arrays, seems fine with all compilers
module math2
implicit none
private
public :: sum_reduce
public :: avg
public :: var
public :: std
contains
pure function std(x) result(y)
implicit none
real, intent(in) :: x(:)
real :: y
y = sqrt(var(x))
end function std
pure function var(x) result(y)
implicit none
real, intent(in) :: x(:)
real :: y
integer :: n
n = size(x)
if (n.gt.1) then
y = sum_reduce((x - avg(x))**2)/(real(n) - 1.0)
else
y = 0.0
end if
end function var
pure function avg(x) result(y)
implicit none
real, intent(in) :: x(:)
real :: y
integer :: n
n = size(x)
if (n.gt.1) then
y = sum_reduce(x)/real(n)
else
y = x(1)
end if
end function avg
pure recursive function sum_reduce(x) result(y)
implicit none
real, intent(in) :: x(:)
real :: y
integer :: n
n = size(x)
if (n.gt.2) then
y = sum_reduce(x(1:(n/2))) + sum_reduce(x(((n/2)+1):n))
else
y = sum(x)
end if
end function sum_reduce
end module math2
==========
I originally thought it was coming from here. Upon further investigation, this was found to not be the case. The following code implements a simple box-muller transform to produce normal random numbers with mean v_mu and standard deviation v_sig:
box-muller transform for normal random numbers
subroutine normrnd(vals, n, v_mu, v_sig)
implicit none
real, parameter :: twopi = 2.0*acos(-1.0)
real, intent(out) :: vals(*)
integer, intent(in) :: n
real, intent(in) :: v_mu, v_sig
real :: u(2), r
integer :: i
do i=1,n
if (mod(i,2).eq.1) then
call random_number(u)
do while (u(1).eq.0.0)
call random_number(u(1))
end do
r = v_sig*sqrt(-2.0*log(u(1)))
vals(i) = r*cos(twopi*u(2)) + v_mu
else
vals(i) = r*sin(twopi*u(2)) + v_mu
end if
end do
end subroutine normrnd
) Yes, that was what I was asking in the above post.