How to compile procedural program effectively

Question

Consider these four function in different form:

f1[ls1_] := Append[Rest@ls1, 0.] + Prepend[Most@ls1, 0.]

f2 = Compile[{{ls1, _Real, 1}}, 
   Append[Rest@ls1, 0.] + Prepend[Most@ls1, 0.], 
   CompilationTarget -> "C"];

f3[ls_] := Module[{n = Length@ls, tp = ls},
  tp[[1]] = ls[[2]];
  tp[[n]] = ls[[n - 1]];
  Do[
   tp[[i]] = ls[[i - 1]] + ls[[i + 1]], {i, 2, n - 1}
   ];
  tp]

f4 = Compile[{{ls, _Real, 1}},
   Module[{n = Length@ls, tp = ls},
    tp[[1]] = ls[[2]];
    tp[[n]] = ls[[n - 1]];
    Do[
     tp[[i]] = ls[[i - 1]] + ls[[i + 1]], {i, 2, n - 1}
     ];
    tp], CompilationTarget -> "C"
   ];

Comparing their performance, we can see that the functional version f1 and f2 have about the same performance, but the compiled procedural version f4 is about 2X slower.

ls = Sin /@ (Range[1, 2000]/2000.);

AbsoluteTiming[TimeConstrained[Do[f1[ls];, {20000}], 5]]
AbsoluteTiming[TimeConstrained[Do[f2[ls];, {20000}], 5]]
AbsoluteTiming[TimeConstrained[Do[f3[ls];, {20000}], 5]]
AbsoluteTiming[TimeConstrained[Do[f4[ls];, {20000}], 5]]

(*{0.976474,Null}*)
(*{0.784508,Null}*)
(*{5.000856,$Aborted}*)
(*{1.595374,Null}*)

However, if we test the procedural version in Fortran (code at the end), it runs about 4X faster than f2.

xslittlegrass$ ifort -r8 -O0 test.f90
    xslittlegrass$ ./a.out
      0.232574000000000
xslittlegrass$

Questions:

1. Why f4 is about 2X slower than f1 and f2?
2. Why f4 is about 7X slower than the Fortran version?
3. Can we make a compiled version that has comparable performance to the Fortran version?
4. I thought that a procedural function compiled to C is comparable in performance compared to procedural functions written in C. But this example seems to disagree with it. So in general how large is the gap of the performance between these two?

---

test.f90

program main
implicit none

integer,parameter :: N0 = 2000
integer i,j
real(kind=8) :: list(N0), result(N0),start,finish

do i = 1, N0
    list(i) = sin(i/real(N0))
end do


call CPU_TIME(start)

do j = 1, 20000
    result(1) = list(2)
    result(N0) = list(N0-1)
    do i=2,N0-1
        result(i) = list(i+1) + list(i-1)
    end do
end do

call CPU_TIME(finish)

write(*,*) finish-start

end program

Answer 1

In the fortran version, you have the large loop within the executable and you call the executable only once. Also quite important: You allocate memory for result only once. In order to make the benchmark fair, you should compare to something like

f5 = Compile[{{ls, _Real, 1}},
  Module[{n = Length@ls, tp},
    tp = Table[0., {Length[ls]}];
    Do[
     tp[[1]] = ls[[2]];
     tp[[n]] = ls[[n - 1]];
     Do[tp[[i]] = ls[[i - 1]] + ls[[i + 1]], {i, 2, n - 1}];
     , {20000}];
    tp],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
];

On my machine, calling f5[ls] takes about 0.067 s while Do[f4[ls];, {20000}] needs about 0.58 s. The fortran variant (see below) compiled with gfortran -o bla test2.f90 -O3 needs about 0.012 s. Moreover, the call f6[ls] with the function defined below needs only 0.034 seconds which is not too far away from the fortran timing.

f6 = Compile[{{ls, _Real, 1}},
   Module[{n = Length@ls, tp},
    tp = Table[0., {n}];
    Do[
     tp[[1]] = Compile`GetElement[ls, 2];
     tp[[n]] = Compile`GetElement[ls, n - 1];
     Do[
      tp[[i]] = Compile`GetElement[ls, i - 1] + Compile`GetElement[ls, i + 1],
      {i, 2, n - 1}], 
    {20000}];
    tp],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

test2.f90

program main
implicit none
integer,parameter :: N0 = 2000
integer i,j
real (kind=8) :: list(N0), result(N0),start,finish
do i = 1, N0
    list(i) = sin(i/real(N0))
end do
call CPU_TIME(start)
do j = 1, 20000
    result(1) = list(2)
    result(N0) = list(N0-1)
    do i=2,N0-1
        result(i) = list(i+1) + list(i-1)
    end do
end do
call CPU_TIME(finish)
write(*,*) result
write(*,*) finish-start
end program

PS.: I have still not used parallelization, here. Even with the dullest way to do that within Mathematica, I get (on a Quad Core CPU):

AbsoluteTiming[ParallelDo[f6[ls], {i, 1, $KernelCount}]][[1]]/($KernelCount)
(* 0.0125135 *)

and by even increasing the number of jobs:

AbsoluteTiming[ParallelDo[f6[ls], {i, 1, $KernelCount 10}]][[ 1]]/($KernelCount 10)
(* 0.00855335 *)

Admittedly, it is not exactly fair to compare this to the unparallelized fortran code. I added this only to show another possibility to speed up the Mathematica code.

Answer 2

This is an answer to Q1 and partly Q4, really. I can't test your Fortran version at the moment, but it would be an interesting comparison.

You can improve the performance of f4 compared to f1 and f2 by setting RuntimeOptions -> "Speed". Clearly the change in runtime settings (mainly "CatchMachineIntegerOverflow" it seems...) from the defaults has a different effect on the two functions.

For instance:

f2 = Compile[{{ls1, _Real, 1}}, 
   Append[Rest@ls1, 0.] + Prepend[Most@ls1, 0.], 
   CompilationTarget -> "C", RuntimeOptions -> "Speed"];

f4 = Compile[{{ls, _Real, 1}}, 
   Module[{n = Length@ls, tp = ls}, tp[[1]] = ls[[2]];
    tp[[n]] = ls[[n - 1]];
    Do[tp[[i]] = ls[[i - 1]] + ls[[i + 1]], {i, 2, n - 1}];
    tp], CompilationTarget -> "C", RuntimeOptions -> "Speed"];

AbsoluteTiming[TimeConstrained[Do[f2[ls];, {20000}], 5]]
(* 0.152 seconds *)

AbsoluteTiming[TimeConstrained[Do[f4[ls];, {20000}], 5]]
(* 0.127 seconds *)