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:
- information about compiler:
GNU Fortran (GCC for Simply Fortran) 11.2.0
Copyright (c) 2021 Free Software Foundation, Inc. - 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