Re: [Caml-list] Closing the performance gap to C

["Christoph Höger" <christoph.hoeger@tu-berlin.de>]

[-- Attachment #1.1: Type: text/plain, Size: 14418 bytes --]

Dear all,

thanks a lot for your valuable input. It seems like there still is
something to gain for OCaml, but it is not quite simple to achieve. From
your comments I reckon the main culprit is the use of a tail-recursive
function instead of a loop, which causes boxing (which causes spilling).
I will repeat my experiment with a new setup asap.

@Ivan: Your point about the optimization is well taken, but please
consider that the runge-kutta integrator is a sophisticated, general
numerical method. Usually, the user of such a method has neither the
knowledge nor the time to manually adapt it, but rather implements a
standard interface. Therefore, the unused arguments and
algebraic/trigonometric identities have to be resolved by the compiler.
Otherwise, we could as well compute the exact solution, anyways ;).

Am 19.12.2016 um 17:59 schrieb Ivan Gotovchits:
>>  Flambda should be able to do the kind of optimizations Ivan mentions,
> at least ideally.
> 
> I've used 4.03+flambda and it didn't perform the optimizations. I've
> played a lot with the inline options and 
> `inline` directive and even noticed a strange (on a first glance) thing,
> that the performance is better, when 
> I explicitly _disable_ of the `trk4_step` function (with `[@@inline never]`)
> 
> As a side note, not about compilers, but about mathematics. 
> Since the three cosine functions, that are computed in the
> loop are basically cos(t), cos(t+h) and cos(t + h/2) and `h` is a
> constant, you can use the  `cos(t+x) = cos(t)cos(x) - sin(t)sin(x)`
> trigonometric
> identity to factor out the constant factors (that are only depend on
> `h`), so finally, you will need to compute one cosine and one sine.
> This will boost the performance of the program, leaving even the SSE
> corei7 specific version far behind.
> 
> 
> On Mon, Dec 19, 2016 at 11:44 AM, Yotam Barnoy <yotambarnoy@gmail.com
> <mailto:yotambarnoy@gmail.com>> wrote:
> 
>     Christoph, you don't mention which version of OCaml you're compiling
>     with, and whether FLambda is used. Flambda should be able to do the
>     kind of optimizations Ivan mentions, at least ideally. If it's not
>     able to do them yet, it will be when the rewritten version is
>     deployed.
> 
>     -Yotam
> 
> 
>     On Mon, Dec 19, 2016 at 10:48 AM, Ivan Gotovchits <ivg@ieee.org
>     <mailto:ivg@ieee.org>> wrote:
>     > Hi Christoph,
>     >
>     > The problem is that your function definitions, like `loop` and
>     `rk4_step`,
>     > have too many parameters
>     > and OCaml is not able to eliminate them and is actually not
>     trying. It was
>     > always a feature of OCaml
>     > that poorly written code will be compiled into a poorly written
>     program.
>     > OCaml never tries to fix
>     > programmer's errors. It will try to minimize the abstraction
>     overhead (often
>     > to the zero). But if the abstractions,
>     > on the first hand, were chosen incorrectly, then it won't fix the
>     code.
>     >
>     > In your particular example, the C compiler was quite aggressive
>     and was
>     > capable of eliminating unnecessary computations.
>     > I wouldn't, however, assume that the C compiler will be always
>     able to do
>     > this for you.
>     >
>     > Let's go through the code:
>     >
>     >
>     > let rec loop steps h n y t =
>     >   if n < steps then loop steps h (n+1) (rk4_step y t h) (t +. h) else
>     >     y
>     >
>     > Here variables `steps` and `h` are loop invariants, so you
>     shouldn't put it
>     > into the loop function parameter list.
>     > Well yes, a compiler can find out loop invariants and remove them
>     for you.
>     > But in our case, to remove them,
>     > it will need to change the type of the function. The compiler will
>     not go
>     > that far for us. It respects a programmer, and if a
>     > programmer decided to make his loop function to depend on `6`
>     arguments,
>     > then it means that the computation is actually depending
>     > on 6 arguments. So, it will allocate 6 registers to hold loop
>     variables,
>     > with all the consequences.
>     >
>     >
>     > Now, let's take a look at the `rk4_step` function:
>     >
>     > let rk4_step y t h =
>     >   let k1 = h *. y' t y in
>     >   let k2 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k1) in
>     >   let k3 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k2) in
>     >   let k4 = h *. y' (t +. h) (y +. k3) in
>     >   y +. (k1+.k4)/.6.0 +. (k2+.k3)/.3.0
>     >
>     >
>     >
>     > This function, is, in fact, a body of the loop, and everything
>     except t is
>     > loop invariant here. Moreover,
>     > function  `y'` is defined as:
>     >
>     > let y' t y = cos t
>     >
>     > I.e., it doesn't really use it the second argument. Probably, a
>     compiler
>     > should inline the call, and eliminate
>     > lots of unecessary computations, and thus free a few registers, but,
>     > apparently, OCaml doesn't do this
>     > (even in 4.03+flambda).
>     >
>     > So we should do this manually:
>     >
>     > let rk4_step y t =
>     >   let k1 = h *. y' t in
>     >   let k2 = h *. y' (t +. 0.5*.h)  in
>     >   let k3 = h *. y' (t +. 0.5*.h)  in
>     >   let k4 = h *. y' (t +. h) (y +. k3) in
>     >   y +. (k1+.k4)/.6.0 +. (k2+.k3)/.3.0
>     >
>     > We can even see, that `k3` and `k2` are equal now, so we can
>     eliminate them:
>     >
>     > let rk4_step y t =
>     >   let k1 = h *. y' t in
>     >   let k2 = h *. y' (t +. 0.5*.h)  in
>     >   let k4 = h *. y' (t +. h) (y +. k3) in
>     >   y +. (k1+.k4)/.6.0 +. k2 *. 1.5
>     >
>     >
>     > Finally, we don't want to pass `y` into the `rk4_step` every time,
>     as we
>     > don't want to require an extra register for it.
>     > After all these manual optimizations, we have the following program:
>     >
>     > let h = 0.1
>     >
>     >
>     > let exact t = sin t
>     >
>     >
>     > let rk4_step t =
>     >
>     >   let k1 = h *. cos t in
>     >
>     >   let k2 = h *. cos (t +. 0.5*.h) in
>     >
>     >   let k4 = h *. cos (t +. h) in
>     >
>     >   (k1+.k4)/.6.0 +. k2*.1.5
>     >
>     >
>     > let compute steps =
>     >
>     >   let rec loop n y t =
>     >
>     >     if n < steps
>     >
>     >     then loop (n+1) (y +. rk4_step t) (t +. h)
>     >
>     >     else y in
>     >
>     >   loop 1 1.0 0.0
>     >
>     >
>     > let () =
>     >
>     >   let y = compute 102 in
>     >
>     >   let err = abs_float (y -. (exact ((float_of_int 102) *. h))) in
>     >
>     >   let large = 50000000 in
>     >
>     >   let y = compute large in
>     >
>     >   Printf.printf "%b\n"
>     >
>     >     (abs_float (y -. (exact (float_of_int large) *. h)) < 2. *. err)
>     >
>     >
>     >
>     > This program has the same performance as the C one... unless I
>     pass really
>     > aggressive optimization options
>     > to the C compiler, that will emit a platform specific code, e.g.,
>     >
>     >     gcc rk.c -lm -O3 -march=corei7-avx -o rksse
>     >
>     >
>     > These options basically double the performance of the C version,
>     leaving
>     > OCaml lagging behind. That is because,
>     > OCaml, obviously, cannot follow the latest developments of intel CPU,
>     > especially in the field of SSE.
>     >
>     > The fun part is that when I've tried to compile the same file with
>     clang,
>     > the resulting program was even slower
>     > than the original non-optimized OCaml. But this is all micro
>     benchmarking of
>     > course, so don't jump to fast conclusions
>     > (although I like to think that OCaml is faster than Clang :))
>     >
>     >
>     > As a final remark, my experience in HPC shows that in general you
>     should not
>     > really rely on compiler optimizations and hope
>     > that the compiler will do the magic for you. Even the GCC
>     compiler. It would
>     > be very easy to accidentally amend the above program
>     > in a such way, that the optimizations will no longer fire in. Of
>     course,
>     > writing in assembly is also not a choice. If you really need
>     > to optimize, then you should find out the performance bottleneck
>     and then
>     > optimize it manually until you get an expected performance.
>     > Alternatively, you can use plain old Fortran to get the reliable
>     > performance. And then call it from C or OCaml.
>     >
>     >
>     > Best wishes,
>     > Ivan
>     >
>     >
>     >
>     > On Mon, Dec 19, 2016 at 6:51 AM, Gerd Stolpmann
>     <info@gerd-stolpmann.de <mailto:info@gerd-stolpmann.de>>
>     > wrote:
>     >>
>     >> Hi Christoph,
>     >>
>     >> the extra code looks very much like an allocation on the minor heap:
>     >>
>     >> sub    $0x10,%r15
>     >> lea    0x25c7b6(%rip),%rax
>     >> cmp    (%rax),%r15
>     >> jb     404a8a <dlerror@plt+0x2d0a>
>     >> lea    0x8(%r15),%rax
>     >> movq   $0x4fd,-0x8(%rax)
>     >>
>     >> r15 points to the used area of the minor heap - by decrementing
>     it you
>     >> get an additional block of memory. It is compared against the
>     beginning
>     >> of the heap to check whether GC is needed. The constant 0x4fd is the
>     >> header of the new block (which must be always initialized).
>     >>
>     >> From the source code, it remains unclear for what this is used.
>     >> Obviously, the compiler runs out of registers, and moves some
>     values to
>     >> the minor heap (temporarily). When you call a C function like cos
>     it is
>     >> likely that this happens because the C calling conventions do not
>     >> preserve the FP registers (xmm*). This could be improved if the OCaml
>     >> compiler tried alternate places for temporarily storing FP values:
>     >>
>     >>  - int registers (which is perfectly possible on 64 bit platforms).
>     >>    A number of int registers survive C calls.
>     >>  - stack
>     >>
>     >> To my knowledge, the OCaml compiler never tries this (but this
>     could be
>     >> out of date). This is a fairly specific optimization that makes
>     mostly
>     >> sense for purely iterating or aggregating functions like yours
>     that do
>     >> not store FP values away.
>     >>
>     >> Gerd
>     >>
>     >> Am Samstag, den 17.12.2016, 14:02 +0100 schrieb Christoph Höger:
>     >> > Ups. Forgot the actual examples.
>     >> >
>     >> > Am 17.12.2016 um 14:01 schrieb Christoph Höger:
>     >> > >
>     >> > > Dear all,
>     >> > >
>     >> > > find attached two simple runge-kutta iteration schemes. One is
>     >> > > written
>     >> > > in C, the other in OCaml. I compared the runtime of both and
>     gcc (-
>     >> > > O2)
>     >> > > produces an executable that is roughly 30% faster (to be more
>     >> > > precise:
>     >> > > 3.52s vs. 2.63s). That is in itself quite pleasing, I think. I do
>     >> > > not
>     >> > > understand however, what causes this difference. Admittedly, the
>     >> > > generated assembly looks completely different, but both compilers
>     >> > > inline
>     >> > > all functions and generate one big loop. Ocaml generates a
>     lot more
>     >> > > scaffolding, but that is to be expected.
>     >> > >
>     >> > > There is however an interesting particularity: OCaml generates 6
>     >> > > calls
>     >> > > to cos, while gcc only needs 3 (and one direct jump).
>     Surprisingly,
>     >> > > there are also calls to cosh, acos and pretty much any other
>     >> > > trigonometric function (initialization of constants, maybe?)
>     >> > >
>     >> > > However, the true culprit seems to be an excess of instructions
>     >> > > between
>     >> > > the different calls to cos. This is what happens between the
>     first
>     >> > > two
>     >> > > calls to cos:
>     >> > >
>     >> > > gcc:
>     >> > > jmpq   400530 <cos@plt>
>     >> > > nop
>     >> > > nopw   %cs:0x0(%rax,%rax,1)
>     >> > >
>     >> > > sub    $0x38,%rsp
>     >> > > movsd  %xmm0,0x10(%rsp)
>     >> > > movapd %xmm1,%xmm0
>     >> > > movsd  %xmm2,0x18(%rsp)
>     >> > > movsd  %xmm1,0x8(%rsp)
>     >> > > callq  400530 <cos@plt>
>     >> > >
>     >> > > ocamlopt:
>     >> > >
>     >> > > callq  401a60 <cos@plt>
>     >> > > mulsd  (%r12),%xmm0
>     >> > > movsd  %xmm0,0x10(%rsp)
>     >> > > sub    $0x10,%r15
>     >> > > lea    0x25c7b6(%rip),%rax
>     >> > > cmp    (%rax),%r15
>     >> > > jb     404a8a <dlerror@plt+0x2d0a>
>     >> > > lea    0x8(%r15),%rax
>     >> > > movq   $0x4fd,-0x8(%rax)
>     >> > >
>     >> > > movsd  0x32319(%rip),%xmm1
>     >> > >
>     >> > > movapd %xmm1,%xmm2
>     >> > > mulsd  %xmm0,%xmm2
>     >> > > addsd  0x0(%r13),%xmm2
>     >> > > movsd  %xmm2,(%rax)
>     >> > > movapd %xmm1,%xmm0
>     >> > > mulsd  (%r12),%xmm0
>     >> > > addsd  (%rbx),%xmm0
>     >> > > callq  401a60 <cos@plt>
>     >> > >
>     >> > >
>     >> > > Is this caused by some underlying difference in the
>     representation
>     >> > > of
>     >> > > numeric values (i.e. tagged ints) or is it reasonable to attack
>     >> > > this
>     >> > > issue as a hobby experiment?
>     >> > >
>     >> > >
>     >> > > thanks for any advice,
>     >> > >
>     >> > > Christoph
>     >> > >
>     >> >
>     >> --
>     >> ------------------------------------------------------------
>     >> Gerd Stolpmann, Darmstadt, Germany    gerd@gerd-stolpmann.de
>     <mailto:gerd@gerd-stolpmann.de>
>     >> My OCaml site:          http://www.camlcity.org
>     >> Contact details:        http://www.camlcity.org/contact.html
>     <http://www.camlcity.org/contact.html>
>     >> Company homepage:       http://www.gerd-stolpmann.de
>     >> ------------------------------------------------------------
>     >>
>     >>
>     >
> 
> 


-- 
Christoph Höger

Technische Universität Berlin
Fakultät IV - Elektrotechnik und Informatik
Übersetzerbau und Programmiersprachen

Sekr. TEL12-2, Ernst-Reuter-Platz 7, 10587 Berlin

Tel.: +49 (30) 314-24890
E-Mail: christoph.hoeger@tu-berlin.de


[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 181 bytes --]