Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

compatibility of results Windows XP/7/8.1/10

Chernov__Sergey
Beginner
681 Views

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?

0 Kudos
11 Replies
andrew_4619
Honored Contributor II
681 Views

I think you should give us some numeric examples of 'significantly different'.

0 Kudos
mecej4
Honored Contributor III
681 Views

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.

0 Kudos
Chernov__Sergey
Beginner
681 Views
In the source code below, there was a significant simplification of my original program in order to reveal the difference between its results and the version of the Windows operating system.
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:
/nologo /Od /warn:declarations /real_size:64 /align:dcommons /module:"x64\Release\\" /object:"x64\Release\\" /Fd"x64\Release\vc140.pdb" /check:all /libs:static /threads /c
 
 
 

 

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 !
   integer(4) :: NOR, NI, NS, NBS, NF
   real(8) :: TSTART, TEND, YNOW(NEQ), HSTART, TOL
end module a
program Everhart
  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(TS.EQ.TF) RETURN
      IF((NOR.LT.2) .OR. (NOR.GT.15)) STOP 'Inadmissible Order!'
      K=NOR/2
      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 I=1,K
        DO J=1,I-1
          SMS(I,J)=1.D0/(SPC(I)-SPC(J))
        END DO
      END DO
      C=0.D0
      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.D0
      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
      DO I=1,K
        FLT=DFLOAT(I+1)
        C(:,I)=C(:,I)/FLT
        D(I,:)=D(I,:)*FLT
        E(I,:)=E(I,:)*FLT
      END DO
      DEG=1.D0/DFLOAT(K+1)
      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
met_2: DO
      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
      IF(NS.LE.1) BL=B
      Q=1.D0
      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
      DO L=1,K
        P=0.D0
        DO M=L,K
          P=P+D(M,L)*B(M,:)
        END DO
        A(L,:)=P
      END DO
      CALL FUN(T,X,F0)
      NF=NF+1
      DO IT=1,IPS
        DO I=1,K
          S=SPC(I)
          P=0.D0
          DO J=K,1,-1
            P=S*(B(J,:)+P)
          END DO
          Y=X+H*S*(F0+P)
          CALL FUN(T+H*S,Y,F); NF=NF+1
          P=(F-F0)/S
          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
        END DO
        IF(ALL(DABS(Y-XK).LE.DABS(EPS*Y))) EXIT
        XK=Y
      END DO
      IF(IT.GT.IPS0) THEN
        IF(CNV.AND.NS.NE.0) NBS=NBS+1
      END IF
      IF(VS) THEN
        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
      P=0.D0
      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
      IF(CNV) THEN
        IPS=IPS0
      ELSE
        IPS=NI
      END IF
      IF(LS) EXIT met_3 ! GOTO 2
      CYCLE met_1 ! GOTO 1
  END DO met_1
END DO met_3
      DEALLOCATE(A,B,C,D,E,BL,SMS,SPC)
      END SUBROUTINE GAUSS_15
!***********************************************************************
0 Kudos
mecej4
Honored Contributor III
681 Views

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        
       
0 Kudos
Chernov__Sergey
Beginner
681 Views

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?

0 Kudos
Steve_Lionel
Honored Contributor III
681 Views

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

0 Kudos
mecej4
Honored Contributor III
681 Views

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.

0 Kudos
JVanB
Valued Contributor II
681 Views

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!

 

0 Kudos
John_Campbell
New Contributor II
682 Views

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.

0 Kudos
Chernov__Sergey
Beginner
676 Views
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
!***********************************************************************

 

0 Kudos
John_Campbell
New Contributor II
676 Views

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.

0 Kudos
Reply