- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
There is a program of numerical modeling (exe-file), written on fortran and compiled by the Intel Visual Fortran Compiler 17.0.2.187 [Intel 64] in the MS Visual Studio 2013 environment for Windows 7 pro X64 (Intel Core 2 Duo).
When you run this program under Windows 10 pro X64 (Intel i3 6100) you get the result significantly different from the result of work under OS Windows 7 pro X64. In this case, the results the program under Windows 8.1 X64 and Windows XP X64 are the same on different computers (Intel Core 2 Duo, Intel Xeon 5130, Intel i5 2400).
What can be wrong?
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I think you should give us some numeric examples of 'significantly different'.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Such things can happen when, for example, (i) the program uses uninitialized variables or exceeds array bounds, (ii) subroutines are called with incorrect arguments, (iii) numerical instability, (iv) differences in compiler versions. The combinations are too many, so you will need to show us a test program and instructions to verify that it does what you say it does.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Here, the difference in the calculation results is observed in the eighth decimal place.
The program was compiled by the Intel Visual Fortran Compiler 17.0.2.187 compiler [Intel 64] in MS Visual Studio 2013 for Windows 7 pro X64 (Intel Core 2 Duo) with the following options:
real(8), parameter :: T0 = 1.405824485291162d0 !
real(8), parameter :: L0 = 6.371110d3 ![km]
real(8), parameter :: V0 = L0/(T0*3.6d3) ![km/s]
!******************************************************************************************
integer(4), parameter :: NEQ = 6
real(8) :: TOut !
integer(4) :: NOR, NI, NS, NBS, NF
real(8) :: TSTART, TEND, YNOW(NEQ), HSTART, TOL
end module a
use a
external FCN, GAUSS_15
!*** Everhart
TSTART = 0.0d0
TOut = 1.7520d4
TEND = TOut/T0
YNOW(1) = 0.668595076910870d0
YNOW(2) = 0.565985496508811d0
YNOW(3) = 1.55503237139677d0
YNOW(4) = -4.36065622514795d0
YNOW(5) = 0.602577749972889d0
YNOW(6) = 1.65556876151629d0
NOR = 15
NI = 0
HSTART = 0.0d0
TOL = 2.3d-16
call GAUSS_15(YNOW, TSTART, TEND, HSTART, TOL, NEQ, NOR, NI, NS, NBS, NF, FCN)
OPEN (10, file = 'result_Everhart.txt', action = 'WRITE', status = 'REPLACE')
WRITE (10,*) 'TEND = ', TEND
WRITE (10,*) 'YNOW(1) = ', YNOW(1)
WRITE (10,*) 'YNOW(2) = ', YNOW(2)
WRITE (10,*) 'YNOW(3) = ', YNOW(3)
WRITE (10,*) 'YNOW(4) = ', YNOW(4)
WRITE (10,*) 'YNOW(5) = ', YNOW(5)
WRITE (10,*) 'YNOW(6) = ', YNOW(6)
WRITE (10, *) 'NS = ', NS
WRITE (10, *) 'NBS = ', NBS
WRITE (10, *) 'NF = ', NF
CLOSE (10)
!*** Everhart
end program Everhart
subroutine FCN(T, Y, YDOT)
use a
real(8) T, Y(6), YDOT(6), R
real(8), parameter :: Pi = 3.141592653589793238462643d0
real(8), parameter :: Cg = Pi*Pi*4.0d0
!**********************************************************************************
R = dsqrt(Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3))
YDOT(1) = Y(4)
YDOT(2) = Y(5)
YDOT(3) = Y(6)
!********************
YDOT(4) = -Cg*Y(1)/R/R/R
!********************
YDOT(5) = -Cg*Y(2)/R/R/R
!********************
YDOT(6) = -Cg*Y(3)/R/R/R !Центральное гравитационное поле Земли
end subroutine FCN
!***********************************************************************
SUBROUTINE GAUSS_15(X,TS,TF,STEP,ERR,N,NOR,NI,NS,NBS,NF,FUN)
! Gauss-Everhart Integrator (Radau & Lobatto Spacings)
! Written by Avdyushev V. (2006)
IMPLICIT REAL*8(A-H,O-Z)
INTEGER(4), PARAMETER :: IPS0 = 100
REAL(8), PARAMETER :: EPS = 1.0D-15
REAL(8), PARAMETER :: SGM = 3.16D0
REAL(8), PARAMETER :: SPAC(56) = &
(/1.00000000000000000000000000000000D0, &
0.66666666666666666666666666666667D0, &
0.50000000000000000000000000000000D0, &
1.00000000000000000000000000000000D0, &
0.35505102572168219018027159252941D0, &
0.84494897427831780981972840747059D0, &
0.27639320225002103035908263312687D0, &
0.72360679774997896964091736687313D0, &
1.00000000000000000000000000000000D0, &
0.21234053823915294397475811012400D0, &
0.59053313555926528913507374793117D0, &
0.91141204048729605260445385623054D0, &
0.17267316464601142810085377187657D0, &
0.50000000000000000000000000000000D0, &
0.82732683535398857189914622812343D0, &
1.00000000000000000000000000000000D0, &
0.13975986434378055215208708112488D0, &
0.41640956763108317994330233133708D0, &
0.72315698636187617231995400231437D0, &
0.94289580388548231780687880744588D0, &
0.11747233803526765357449851302033D0, &
0.35738424175967745184292450297956D0, &
0.64261575824032254815707549702044D0, &
0.88252766196473234642550148697967D0, &
1.00000000000000000000000000000000D0, &
0.09853508579882642612349889788775D0, &
0.30453572664636390548538517627883D0, &
0.56202518975261385599498747999477D0, &
0.80198658212639182746420786320470D0, &
0.96019014294853125765919330990667D0, &
0.08488805186071653506398389301626D0, &
0.26557560326464289309811405904562D0, &
0.50000000000000000000000000000000D0, &
0.73442439673535710690188594095438D0, &
0.91511194813928346493601610698373D0, &
1.00000000000000000000000000000000D0, &
0.07305432868025888514812603418031D0, &
0.23076613796994549908311663988435D0, &
0.44132848122844986791860665819448D0, &
0.66301530971884570090294702791922D0, &
0.85192140033151570815002314750402D0, &
0.97068357284021510802794972308684D0, &
0.06412992574519669233127711938966D0, &
0.20414990928342884892774463430102D0, &
0.39535039104876056561567136982732D0, &
0.60464960895123943438432863017268D0, &
0.79585009071657115107225536569898D0, &
0.93587007425480330766872288061033D0, &
1.00000000000000000000000000000000D0, &
0.05626256053692214646565219103231D0, &
0.18024069173689236498757994280918D0, &
0.35262471711316963737390777017124D0, &
0.54715362633055538300144855765235D0, &
0.73421017721541053152321060830661D0, &
0.88532094683909576809035976293249D0, &
0.97752061356128750189117450042915D0/)
REAL(8), ALLOCATABLE :: A(:,:),B(:,:),C(:,:),D(:,:),E(:,:)
REAL(8), ALLOCATABLE :: SMS(:,:),SPC(:),BL(:,:)
INTEGER(4) :: N, NOR, K, L, I, J, NI, NS, NBS, NST, NF, IPS, M, IT
REAL(8), DIMENSION(N) :: X,Y,F,F0,P,XK
LOGICAL LS,VS,INC,CNV
REAL(8) :: TS, TF, FLT, DEG, RMX, RMN, DIR, H_, STEP, ERR, T, H, R, Q, B_, S
EXTERNAL FUN
IF((NOR.LT.2) .OR. (NOR.GT.15)) STOP 'Inadmissible Order!'
ALLOCATE(SPC(0:K),SMS(K,K))
ALLOCATE(A(K,N),B(K,N),BL(K,N))
ALLOCATE(C(0:K,0:K),D(0:K,0:K),E(0:K,0:K))
SPC(0)=0.D0; L=K*(NOR-K-1); SPC(1:K)=SPAC(L+1:L+K)
DO J=1,I-1
SMS(I,J)=1.D0/(SPC(I)-SPC(J))
END DO
END DO
C(0,0)=1.D0
D=0.D0
D(0,0)=1.D0
DO J=1,K
DO I=J,K
C(I,J)=C(I-1,J-1)-SPC(I-1)*C(I-1,J)
D(I,J)=D(I-1,J-1)+SPC(J )*D(I-1,J)
END DO
END DO
E(0:K,0)=1.D0
DO I=1,K
DO J=1,I
E(I,J)=E(I-1,J-1)+E(I-1,J)
END DO
END DO
FLT=DFLOAT(I+1)
C(:,I)=C(:,I)/FLT
D(I,:)=D(I,:)*FLT
E(I,:)=E(I,:)*FLT
END DO
RMX=SGM**DEG
RMN=1.D0/RMX
DIR=(TF-TS)/DABS(TF-TS)
H_=DABS(STEP)
VS=ERR.NE.0.D0
CNV=NI.LE.0
INC=.TRUE.
T=TS
NS=0
NBS=0
NST=0
NF=0
IPS=IPS0
met_3: DO
H=H_*DIR
R=1.D0
LS=.FALSE.
B=0.D0
XK=0.D0
IF(H_.EQ.0.D0) THEN
IF(.NOT.VS) STOP 'Zero Constant Step!'
CALL FUN(TS,X,F0)
NF=NF+1
H_=EPS
B_=0.D0
DO WHILE(B_.EQ.0.D0)
H_=H_*10.D0
H=H_*DIR
Y=X+H*F0
IF(DABS(TF-TS).LE.H_) EXIT met_2 !GOTO 1
CALL FUN(TS+H,Y,F)
NF=NF+1
DO J=1,N
B_=B_+(F(J)-F0(J))**2
END DO
END DO
H_=DSQRT(2.D0*ERR*H_/DSQRT(B_))
H=H_*DIR
END IF
EXIT met_2
END DO met_2
met_1: DO
IF(DABS(TF-T).LE.H_) THEN
STEP=H_
H=TF-T
H_=DABS(H)
R=R*H_/STEP
LS=.TRUE.
END IF
DO L=1,K
P=0.D0
DO M=L,K
P=P+E(M,L)*B(M,:)
END DO
B(L,:)=B(L,:)-BL(L,:)
Q=Q*R
BL(L,:)=P*Q/DFLOAT(L+1)
B(L,:)=BL(L,:)+B(L,:)
END DO
P=0.D0
DO M=L,K
P=P+D(M,L)*B(M,:)
END DO
A(L,:)=P
END DO
NF=NF+1
P=0.D0
DO J=K,1,-1
P=S*(B(J,:)+P)
END DO
Y=X+H*S*(F0+P)
DO J=1,I-1
P=(P-A(J,:))*SMS(I,J)
END DO
DO J=1,I
B(J,:)=B(J,:)+C(I,J)*(P-A(I,:))
END DO
A(I,:)=P
XK=Y
IF(CNV.AND.NS.NE.0) NBS=NBS+1
END IF
B_=0.D0
DO J=1,N
B_=B_+B(K,J)**2
END DO
IF(B_.NE.0.D0) THEN
R=(ERR/DSQRT(B_)/H_)**DEG
ELSE
R=RMX
END IF
met_4: DO
IF(NS.EQ.0) THEN
INC=.NOT.LS.AND.INC
IF(INC.AND.R.GT.RMX.OR.R.LT.RMN) THEN
NST=NST+1
IF(NST.GT.100.AND.R.GT.RMX) EXIT met_4 !GOTO 4
H_=H_*R
CYCLE met_3 !GOTO 3
END IF
END IF
EXIT met_4
END DO met_4
IF(R.GT.RMX) R=RMX
END IF
DO J=K,1,-1
P=B(J,:)+P
END DO
X=X+H*(F0+P)
T=T+H
NS=NS+1
H=H*R
H_=DABS(H)
XK=0.D0
IPS=IPS0
ELSE
IPS=NI
END IF
CYCLE met_1 ! GOTO 1
END DO met_1
END DO met_3
!***********************************************************************
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The discrepancies are inherent in the problem that you are solving with your chosen integrator.
The program integrates the equation of motion of a particle attracted by gravity to a center of mass with an initial velocity vector that does not pass through the center. The integrator adjusts the step-size based on intermediate results. All the usual vagaries of floating-point calculation are at play. The small differences in the results are not really related to the OS version. On the Windows 10-64 laptop that I am writing this response on, I get results that agree to about 7 digits when I run the program using different compilers (Ifort 17.0.6 and Gfortran 6.2).
TEND = 12462.4376537100 TEND = 12462.43765370996 YNOW(1) = 0.700009814515993 YNOW(1) = 0.7000098372565887 YNOW(2) = -0.561524885579132 YNOW(2) = -0.5615248822633300 YNOW(3) = -1.54277694359061 YNOW(3) = -1.542776934480667 YNOW(4) = 4.32628928299755 YNOW(4) = 4.326289257444321 YNOW(5) = 0.630890584685439 YNOW(5) = 0.6308906051805072 YNOW(6) = 1.73335763557065 YNOW(6) = 1.733357691879451 NS = 48544224 NS = 48381581 NBS = 0 NBS = 0 NF = 728233634 NF = 725792778
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
When different integrators or different compilers are used, I can understand minor differences in the results.
But why with the same integrator with the same compiler exe-file on WIN 7 and WIN 10 gives different results?
It seemed to me that the work process of the exe-file is deterministic. Or am I wrong?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Different processors may use different paths through the math library. The OS is not the only variable here. See http://sc13.supercomputing.org/sites/default/files/WorkshopsArchive/pdfs/wp129s1.pdf
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Chernov, Sergey wrote:
When different integrators or different compilers are used, I can understand minor differences in the results.But why with the same integrator with the same compiler exe-file on WIN 7 and WIN 10 gives different results?
It seemed to me that the work process of the exe-file is deterministic. Or am I wrong?
It is deterministic, but there is more than one path! These days, more of the Fortran library resides in a DLL rather than in the EXE itself. If the different OSes that you used had different Intel Fortran RTL DLLs, that would be part of the reason for the differences. Windows has some control over which instruction sets are allowed to execute, and sometimes this choice can be modified by the user using the BCDEDIT command.
That the Windows version is different is the most visible difference, but may be unimportant. There are so many cogs in the machine, and it is not easy to pinpoint the most significant cause of the differences.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Just for fun I wrote a little program that checked the entries of the SPAC array. Turns out SPAC(31). SPAC(37), and SPAC(43) were low by 1 in the last place, i.e. 1.0e-32, confirmed by comparing with http://keisan.casio.com/exec/system/15041456505628 . So there is no problem with the SPAC array!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Out of interest, I rewrote subroutine FUN to the following (to direct some vectorizing)
subroutine FCN(T, Y, YDOT) use a real*8 T, Y(6), YDOT(6), RRR, rr, x real*8, parameter :: Pi = 3.141592653589793238462643d0 real*8, parameter :: Cg = Pi*Pi*4.0d0 !********************************************************************************** YDOT(1) = Y(4) YDOT(2) = Y(5) YDOT(3) = Y(6) !******************** !z R = dsqrt(Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)) !z YDOT(6) = -Cg*Y(3)/R/R/R ! RR = ( Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3) ) RRR = SQRT ( RR ) * RR x = -Cg/RRR YDOT(4) = Y(1)*x YDOT(5) = Y(2)*x YDOT(6) = Y(3)*x !******************** end subroutine FCN
This produces a change in values, of the order that you describe. Any change to use of vector instructions (with a change of processor) would produce a similar result.
{old FUN} cycle_met_1 = 39278734 cycle_met_2 = 12 cycle_met_3 = 12 cycle_met_4 = 39278734 TEND = 12462.437653709958 YNOW(1) = 0.70000985149761608 YNOW(2) = -0.56152488018610514 YNOW(3) = -1.5427769287744608 YNOW(4) = 4.3262892414438259 YNOW(5) = 0.63089061801514912 YNOW(6) = 1.7333577271451566 NS = 39278723 NBS = 0 NF = 589242241 {new FUN} cycle_met_1 = 36575171 cycle_met_2 = 12 cycle_met_3 = 12 cycle_met_4 = 36575171 call_fun = 548685093 TEND = 12462.437653709958 YNOW(1) = 0.70000984092209018 YNOW(2) = -0.56152488172903547 YNOW(3) = -1.5427769330129555 YNOW(4) = 4.3262892533261299 YNOW(5) = 0.63089060848288026 YNOW(6) = 1.7333577009535404 NS = 36575160 NBS = 0 NF = 548685093 sec = 107.47313891611338
There is considerable variation in the number of iterations being reported in the different results above, which suggests that the convergence tests could be reviewed.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The replacement of the integrator by the one given below leads, under the same initial conditions, to a complete coincidence of the results in Windows 7 and Windows 10: module a real(8), parameter :: T0 = 1.405824485291162d0 !real(8), parameter :: L0 = 6.371110d3 ![km] real(8), parameter :: V0 = L0/(T0*3.6d3) ![km/s] !****************************************************************************************** integer(4), parameter :: NEQ = 6 real(8) :: TOut ! real(8) :: TSTART, TNOW, TEND, YNOW(NEQ), H_step_0, H_step, TOL integer(4) :: N_step end module a program RKM use a external FCN, One_Step TSTART = 0.0d0 TOut = 1.7520d4 TEND = TOut/T0 YNOW(1) = 0.668595076910870d0 YNOW(2) = 0.565985496508811d0 YNOW(3) = 1.55503237139677d0 YNOW(4) = -4.36065622514795d0 YNOW(5) = 0.602577749972889d0 YNOW(6) = 1.65556876151629d0 TOL = 2.3D-16 H_step_0 = 1.0D-2 H_step = H_step_0 N_step = 0 !*** RKM TNOW = TSTART integr: do call One_Step(TNOW, YNOW, H_step, TOL) N_step = N_step + 1 if (TNOW >= TEND ) exit integr end do integr OPEN (10, file = 'result_RKM.txt', action = 'WRITE', status = 'REPLACE') WRITE (10,*) 'TOL = ', TOL WRITE (10,*) 'H_step_0 = ', H_step_0 WRITE (10,*) ' H_step = ', H_step WRITE (10,*) 'TNOW = ', TNOW WRITE (10,*) 'YNOW(1) = ', YNOW(1) WRITE (10,*) 'YNOW(2) = ', YNOW(2) WRITE (10,*) 'YNOW(3) = ', YNOW(3) WRITE (10,*) 'YNOW(4) = ', YNOW(4) WRITE (10,*) 'YNOW(5) = ', YNOW(5) WRITE (10,*) 'YNOW(6) = ', YNOW(6) WRITE (10,*) 'N_step = ', N_step CLOSE (10) end program RKM subroutine FCN(T, Y, YDOT) use a real(8) T, Y(6), YDOT(6), R real(8), parameter :: Pi = 3.141592653589793238462643d0 real(8), parameter :: Cg = Pi*Pi*4.0d0 !********************************************************************************** R = dsqrt(Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)) YDOT(1) = Y(4) YDOT(2) = Y(5) YDOT(3) = Y(6) !******************** YDOT(4) = -Cg*Y(1)/R/R/R !******************** YDOT(5) = -Cg*Y(2)/R/R/R !******************** YDOT(6) = -Cg*Y(3)/R/R/R !Центральное гравитационное поле Земли end subroutine FCN !*********************************************************************** subroutine One_Step(T, Y, H, E) external FCN integer(4) :: j real(8) :: T, T_temp, Y(6), Y_temp(6), H, E, H3, F(6), K(6, 5), Err(6), MaxErr do H3 = H/3.0d0 T_temp = T Y_temp = Y call FCN(T_temp, Y_temp, F) K(1,1) = F(1)*H3 K(2,1) = F(2)*H3 K(3,1) = F(3)*H3 K(4,1) = F(4)*H3 K(5,1) = F(5)*H3 K(6,1) = F(6)*H3 T_temp = T + H3 Y_temp = Y + K(1:6,1) call FCN(T_temp, Y_temp, F) K(1,2) = F(1)*H3 K(2,2) = F(2)*H3 K(3,2) = F(3)*H3 K(4,2) = F(4)*H3 K(5,2) = F(5)*H3 K(6,2) = F(6)*H3 T_temp = T + H3 Y_temp = Y + (K(1:6,1) + K(1:6,2))/2.0d0 call FCN(T_temp, Y_temp, F) K(1,3) = F(1)*H K(2,3) = F(2)*H K(3,3) = F(3)*H K(4,3) = F(4)*H K(5,3) = F(5)*H K(6,3) = F(6)*H T_temp = T + H/2.0d0 Y_temp = Y + 0.375d0*(K(1:6,1) + K(1:6,3)) call FCN(T_temp, Y_temp, F) K(1,4) = K(1,1) + 4.0d0*H3*F(1) K(2,4) = K(2,1) + 4.0d0*H3*F(2) K(3,4) = K(3,1) + 4.0d0*H3*F(3) K(4,4) = K(4,1) + 4.0d0*H3*F(4) K(5,4) = K(5,1) + 4.0d0*H3*F(5) K(6,4) = K(6,1) + 4.0d0*H3*F(6) T_temp = T + H Y_temp = Y + 1.5d0*(K(1:6,4) - K(1:6,3)) call FCN(T_temp, Y_temp, F) K(1,5) = F(1)*H3 K(2,5) = F(2)*H3 K(3,5) = F(3)*H3 K(4,5) = F(4)*H3 K(5,5) = F(5)*H3 K(6,5) = F(6)*H3 do j = 1, 6 Err(j) = dabs((2.0d0*K(j,4) - 3.0d0*K(j,3) - K(j,5))/10.0d0) end do MaxErr = Err(1) do j = 2, 6 if (Err(j) > MaxErr) then MaxErr = Err(j) end if end do if (MaxErr <= E) then T = T + H do j = 1, 6 Y(j) = Y(j) + (K(j,4) + K(j,5))/2.0d0 end do exit else H = H/2.0d0 end if end do end subroutine One_Step !***********************************************************************
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I think you choice of an alternate integrator is a good approach.
Previously you were using a high order polynomial approximation to get a result, but calling the function about 700 million times. This does not appear to be a good result.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page