1 files changed, 60 insertions, 0 deletions
diff --git a/libs/numeric/odeint/performance/fortran_lorenz.f90 b/libs/numeric/odeint/performance/fortran_lorenz.f90
new file mode 100644
index 000000000..26869973c
--- /dev/null
+++ b/libs/numeric/odeint/performance/fortran_lorenz.f90
@@ -0,0 +1,60 @@
+program main
+    implicit none
+
+    integer, parameter :: dp = 8
+    real(dp), dimension(1:3) :: x
+    integer, parameter :: nstep = 20000000
+    real(dp) :: t = 0.0_dp
+    real(dp) :: h = 1.0e-10_dp
+    integer, parameter :: nb_loops = 21
+    integer, parameter :: n = 3
+    integer :: k
+    integer :: time_begin
+    integer :: time_end
+    integer :: count_rate
+    real(dp) :: time
+    real(dp) :: min_time = 100.0
+
+    do k = 1, nb_loops
+        x = [ 8.5_dp, 3.1_dp, 1.2_dp ]
+        call system_clock(time_begin, count_rate)
+        call rk4sys(n, t, x, h, nstep)
+        call system_clock(time_end, count_rate)
+        time = real(time_end - time_begin, dp) / real(count_rate, dp)
+        min_time = min(time, min_time)
+        write (*,*) time, x(1)
+    end do
+    write (*,*) "Minimal Runtime:", min_time
+contains
+    subroutine xpsys(x,f)
+        real(dp), dimension(1:3), intent(in) :: x
+        real(dp), dimension(1:3), intent(out) :: f
+        f(1) = 10.0_dp * ( x(2) - x(1) )
+        f(2) = 28.0_dp * x(1) - x(2) - x(1) * x(3)
+        f(3) = x(1) * x(2) - (8.0_dp / 3.0_dp) * x(3)
+    end subroutine xpsys
+
+    subroutine rk4sys(n, t, x, h, nstep)
+        integer, intent(in) :: n
+        real(dp), intent(in) :: t
+        real(dp), dimension(1:n), intent(inout) :: x
+        real(dp), intent(in) :: h
+        integer, intent(in) :: nstep
+        ! Local variables
+        real(dp) :: h2
+        real(dp), dimension(1:n) :: y, f1, f2, f3, f4
+        integer :: i, k
+
+        h2 = 0.5_dp * h
+        do k = 1, nstep
+            call xpsys(x, f1)
+            y = x + h2 * f1
+            call xpsys(y, f2)
+            y = x + h2 * f2
+            call xpsys(y, f3)
+            y = x + h * f3
+            call xpsys(y, f4)
+            x = x + h * (f1 + 2.0_dp * (f2 + f3) + f4) / 6.0_dp
+        end do
+    end subroutine rk4sys
+end program main