A strange problem met when using Fortran

When I’m trying to do llinearfitting by Fortran, I met a strange problem, my code is as follows:

program fit
use numerical_s
implicit none
integer, parameter :: n1 = 31, n2 = 34, n3 = 11
real(kind=8) :: lower = 0.0, upper1 = 30.0, upper2 = 33.0, upper3 = 10.0
real(kind=8), dimension(n1) :: y1, x11, x12
real(kind=8), dimension(n2) :: y2, x21, x22
real(kind=8), dimension(n3) :: y3, x31, x32
real(kind=8), dimension(6) :: a, b, r
integer :: i

y1 = Linespace(lower, upper1, n1)
y2 = Linespace(lower, upper2, n2)
y3 = Linespace(lower, upper3, n3)

x11 = (/00.0, 01.7, 03.0, 04.4, 05.4, 06.6, 07.7, 08.7, &
        09.8, 10.8, 11.6, 12.4, 13.5, 14.5, 15.5, 16.4, &
        17.3, 18.2, 19.0, 19.8, 20.5, 21.4, 22.3, 23.4, &
        24.0, 25.6, 26.0, 27.0, 28.1, 28.8, 29.6/)
x12 = (/00.3, 01.3, 02.0, 02.9, 03.6, 04.4, 05.3, 06.0, &
        06.8, 07.7, 08.5, 09.5, 10.3, 11.3, 12.1, 13.0, &
        14.0, 15.0, 15.8, 17.0, 17.9, 18.8, 20.1, 21.0, &
        22.0, 23.2, 24.3, 25.9, 27.0, 28.6, 30.0/)

x21 = (/00.0, 02.2, 03.2, 04.4, 05.2, 06.4, 07.4, 09.0, 09.7, &
        11.3, 13.0, 14.6, 16.1, 16.7, 17.1, 17.6, 18.4, 18.9, &
        19.7, 20.3, 21.0, 21.7, 22.7, 23.4, 24.2, 24.9, 25.8, &
        26.4, 27.0, 27.7, 28.0, 28.6, 29.2, 29.8/)
x22 = (/00.3, 00.8, 01.5, 02.3, 03.1, 03.8, 04.4, 05.2, 06.2, &
        07.2, 08.1, 09.2, 10.0, 10.7, 11.5, 12.5, 13.4, 14.2, &
        14.8, 15.7, 16.7, 17.6, 18.3, 19.4, 20.5, 21.6, 22.4, &
        23.3, 24.5, 25.7, 26.9, 28.1, 29.4, 30.0/)

x31 = (/0.0, 8.0, 15.0, 21.0, 27.6, 32.0, 37.0, 42.0, 47.0, 52.0, 57.0/)
x32 = (/3.0, 8.0, 13.0, 18.0, 23.0, 28.0, 33.0, 40.0, 46.0, 52.0, 60.0/)

call Lin_fitting(x11, y1, n1, a(1), b(1), r(1))
call Lin_fitting(x12, y1, n1, a(2), b(2), r(2))
call Lin_fitting(x21, y2, n2, a(3), b(3), r(3))
call Lin_fitting(x22, y2, n2, a(4), b(4), r(4))
print *, ' '
call Lin_fitting(x31, y3, n3, a(5), b(5), r(5))
call Lin_fitting(x32, y3, n3, a(6), b(6), r(6))

open(unit = 101, file = 'data.txt', status = 'replace')
201 format(3(F015.8))
202 format(31(F015.8))
203 format(34(F015.8))
204 format(11(F015.8))

write(101, 202) (y1(i), i = 1, n1)
write(101, 202) (x11(i), i = 1, n1)
write(101, 202) (x12(i), i = 1, n1)
write(101, 203) (y2(i), i = 1, n2)
write(101, 203) (x21(i), i = 1, n2)
write(101, 203) (x22(i), i = 1, n2)
write(101, 204) (y3(i), i = 1, n3)
write(101, 204) (x31(i), i = 1, n3)
write(101, 204) (x32(i), i = 1, n3)
write(101, 201) a(1), b(1), r(1)
write(101, 201) a(2), b(2), r(2)
write(101, 201) a(3), b(3), r(3)
write(101, 201) a(4), b(4), r(4)
write(101, 201) a(5), b(5), r(5)
write(101, 201) a(6), b(6), r(6)
!write(*, 203) (yt(i), i = 1, 2*n)

end program fit

function and subroutine used:
subroutine Lin_fitting(x, y, dim, a, b, r)
integer, intent(in) :: dim
real(kind=8), dimension(dim), intent(in) :: x, y
real(kind=8), intent(out) :: a, b, r
real(kind=8) :: average_x, average_y, sum1, sum2, sum3
integer :: i
call Average(x, dim, average_x)
call Average(y, dim, average_y)
do i = 1, dim
sum1 = sum1 + (x(i) - average_x) ** 2
sum2 = sum2 + (x(i) - average_x) * (y(i) - average_y)
sum3 = sum3 + (y(i) - average_y) ** 2
end do
b = sum2 / sum1
a = average_y - b * average_x
r = sum2 ** 2 / (sum1 * sum3)
end subroutine Lin_fitting

subroutine Average(in, dim, out)
integer, intent(in) :: dim
real(kind=8), dimension(dim), intent(in) :: in
real(kind=8), intent(out) :: out
real(kind=8) :: sum_
call Sum(in, dim, sum_)
out = sum_ / dim
end subroutine Average

function Linespace(low, up, step) result(retval)
real(kind=8), intent(in) :: low, up
integer, intent(in) :: step
real(kind=8), dimension(step):: retval
real(kind=8) :: lenos
integer :: i
lenos = (up - low) / (step - 1)
do i = 1, step - 1
retval(i) = low + lenos * (i-1)
end do
retval(step) = up
end function Linespace

when I seek the result, I found that a(5), a(6), b(5), b(6), r(5), r(6) are not at a right number:
a(5) = -21.08166659, b(5) = 0.84730754, r(5) = 0.81409002
a(6) = -15.87139223, b(6) = 0.70859665, r(6) = 0.71153680

but the strange thing is when I add the sentence: print *, ’ ’ , all the results get ideal:
a(5) = -0.52758477, b(5) = 0.17957304, r(5) = 0.99320218
a(6) = -0.26949830, b(6) = 0.17890272, r(6) = 0.99404857

further test shows that, if I put the sentence: print *, ‘\n’ between the setences like this:
call Lin_fitting(x31, y3, n3, a(5), b(5), r(5))
print *, ’ ’
call Lin_fitting(x32, y3, n3, a(6), b(6), r(6))
the result will become:
a(5) = -21.08166659, b(5) = 0.84730754, r(5) = 0.81409002
a(6) = -0.26949830, b(6) = 0.17890272, r(6) = 0.99404857

This problem really puzzled me alot, I had no idea about how could it appears. Could anyone tells me what caused this problem?

The additional information is as follows:

  1. information about compiler:
    GNU Fortran (GCC for Simply Fortran) 11.2.0
    Copyright (c) 2021 Free Software Foundation, Inc.
  2. command used to complie:
    gfortran.exe -c numerical_f.f90 -o numerical_f.o

gfortran.exe -c numerical_s.f90 -o numerical_s.o

gfortran.exe fit.f90 -o fit.exe numerical_s.o numerical_f.o

.\fit.exe

If my understanding is correct, the Average subroutine should be the following:

  subroutine Average(in, dim, out)
      implicit none
      integer, intent(in) :: dim
      real(kind=8), dimension(dim), intent(in) :: in
      real(kind=8), intent(out) :: out
      real(kind=8) :: sum_tmp
      sum_tmp = Sum(in, dim=1)
      out = sum_tmp / dim
   end subroutine Average

The keyword argument dim of the Sum function should correctly specify the number of the dimension of the array, so it would not be appropriate to pass the variable dim used in this code to the SUM function.

In short: in Lin_fitting you are using the sum1, sum2, sum3 variables without having initialized them before. If you initialize them to zero before the loop, you always get the right result.

Understanding why inserting print *,'' between the calls “fixes” the problem is another story, but the point is: the behavior of any program when using an uninitialized variable is undefined (unpredictable results).

Welcome to the forum!

I have merely glanced at your program’s code, but the symptoms you describe may be caused by uninitialised variables (among many other things, of course). However, in routine lin_spacing you use sum1, sum2 and sum3 without first setting them to zero.

In general: if a variable is not explicitly assigned a value, it can have any value as Fortran does not initialise them to zero. So, try the program again with these three local variables set to zero and see if the results are now correct.

There may be more errors lurking in the code, but these are certainly candidates.

Thank you very much, the result get ideal after initialize sum1, sum2, and sum3.

Additional comments

Could you please edit your post and enclose the whole code between triple backticks (```) to have it properly formatted? As it is, it a bit difficult to read.

Golden rule #1: always put implicit none at the top of your programs and modules, particularly when debugging. It helps catching undeclared variables, which are a significant source of bugs (even if it was not the problem in this code)

Golden rule #2: when debugging, always enable the runtime checking options of your compiler. This is -fcheck=all with gfortran. That said, gfortran doesn’t catch the error here, I had to use the ifort complier with the -check option, and it caught the initialized variable issue.

If possible, when asking debugging help, try posting a code that can be directly copy/pasted and compiled, so that people here can try running it on their own.

Fortran has a built-in sum function, you don’t need to redefined it. In you code you may write:
average_x = sum(x) / dim

Rather avoid using a variable names dim: it’s perfectly legal, but it introduces some confusion, because dim is also an optional keyword of several intrinsic functions (including sum(), by the way).

real(kind=8) is not portable in theory, even if it works on most compilers. If you want a 64 bits real, better use real(kind=real64), real64being a constant provided by the iso_fortan_env module. Or define the kind by yourself:

integer, parameter :: wp = selected_real_kind(p=12)
real(kind=wp) :: ...
1 Like

Yes, you are right, thank you very much. I just forget to initiallize the variable before using them, after doing like this, the problem is solved. What’s more, thankyou for your advice, it’s my first thime to use this forum, next time, I will pay more attention about these details. At last, could you tell me why inserting print *,'' between the calls “fixes” the problem? Thankyou in advance.

Local variables such as sum1 only “live” during the execution of the lin_fitting routine. However, when exiting the routine the value is still in the memory slot that was used. When entering again the routine the same memory slot is used, and the value at the end of the previous call is still there. So, if you don’t initialize the variable you get wrong results.

When you insert print *,'' between the 2 calls, the print statement likely uses the memory slot that was attributed to sum1 before, to store its own local variables, and for some reason sets it to zero. This is by pure chance, because it could as well set it to any other value…

1 Like

gfortran has had the option -finit-real=snan for some years now to signal when there is an attempt to use an unitialized real. As far as I can tell it isn’t at all costly, so I leave it on. No similar option for integers, although I would gladly give up one value in the range of integers for such a feature.