caml (special) light and numerics

caml (special) light and numerics

[Thorsten Ohl]

A suggestion for caml special light v1.07:

   -  on Linux (and other i386) systems, the runtime should be linked
      with `-lieee' to enable IEEE arithmetic (`-mieee-fp' does the
      same thing).  This way the runtime will not abort on division by
      0.0, but give the result `Infinity' instead.  Predefined
      constants for `Infinity' and `NaN' would be handy, but can be
      prepared in the program.

Does anybody have numbers on how csl (in particular the new native
compiler) compares in performance to C++ or FORTRAN in numerical
applications (quadrature, ODEs, linear algebra, ...)?

Cheers,
-Thorsten

/// Thorsten Ohl, TH Darmstadt, Schlossgartenstr. 9, D-64289 Darmstadt, Germany
//////////////// net: ohl@crunch.ikp.physik.th-darmstadt.de, ohl@gnu.ai.mit.edu
/// voice: +49-6151-16-3116, secretary: +49-6151-16-2072, fax: +49-6151-16-2421

Re: caml (special) light and numerics

[Xavier Leroy]

> A suggestion for caml special light v1.07:
>    -  on Linux (and other i386) systems, the runtime should be linked
>       with `-lieee' to enable IEEE arithmetic (`-mieee-fp' does the
>       same thing).

Yes, the next release will run in IEEE mode under Linux and FreeBSD.

> Does anybody have numbers on how csl (in particular the new native
> compiler) compares in performance to C++ or FORTRAN in numerical
> applications (quadrature, ODEs, linear algebra, ...)?

I have two numerical benchmarks that are available both in CSL and in C:
a simple FFT implementation and the Pseudoknot benchmark used by
Feeley, Hartel et al (a nucleotide placement problem). The figures on
a DecStation 3000/300X are as follows:

                        C       CSL
        Pseudoknot      0.73s   1.27s
        FFT             3.07s   3.74s

Your mileage may vary greatly, in particular depending on how
carefully the CSL programs are written. Also, the I386 port of CSL is
not as good as the other ports on floating-point operations, because
CSL assumes a standard, direct-access set of floating-point registers,
while the I387 imposes a stack organization on its floating-point
registers.

- Xavier Leroy

Re: caml (special) light and numerics

[Thorsten Ohl]

>>>>> "Xavier" == Xavier Leroy <Xavier.Leroy@inria.fr> writes:

Xavier> Your mileage may vary greatly, in particular depending on how
Xavier> carefully the CSL programs are written.

This may be asking too much, but are there any (written) guidelines
for efficient numerical CSL coding (besides the one on the WWW site
that gave me some hints)?  I'm trying to find out if a CSL like
language can be useful in physics simulations.  There will be *some*
speed penalty, but the important question is whether the elegance of a
modern language is worth this penalty.

I guess what I'm looking for are also some examples, because the
typical tutorials and textbooks stay away from numerical applications.

An example: the FORTRAN or C programmers in me sometimes want an
array local to a function, which would be cheap, because it's
allocated on the stack.  But CSL will allocate it on the heap if I use
`Array.new' and I assume that this is expensive.  Is this true and are
there any standard procedures around it?  (Typical example)

  val f : float array -> float array

will certainly be expensive, but a `subroutine'

  val g : float array -> unit

modifying the vector in place will also (in general) need a temporary
array to store either input or output.

Another example: how can I find out (without deciphering assembler
code) if a particular recursive function has been recognized as tail
recursive?

Xavier> Also, the I386 port of CSL is not as good as the other ports
Xavier> on floating-point operations

I can imaginee, but I'm still wating for these cheap Alpha-clones to
take off.  Unfortunately most fast machines I have an account on are
either HP-PA or RS6000 ...

Thanks for your attention,
-Thorsten

/// Thorsten Ohl, TH Darmstadt, Schlossgartenstr. 9, D-64289 Darmstadt, Germany
//////////////// net: ohl@crunch.ikp.physik.th-darmstadt.de, ohl@gnu.ai.mit.edu
/// voice: +49-6151-16-3116, secretary: +49-6151-16-2072, fax: +49-6151-16-2421

Re: caml (special) light and numerics

[Pierre Weis]

> I guess what I'm looking for are also some examples, because the
> typical tutorials and textbooks stay away from numerical applications.
> 
> An example: the FORTRAN or C programmers in me sometimes want an
> array local to a function, which would be cheap, because it's
> allocated on the stack.  But CSL will allocate it on the heap if I use
> `Array.new' and I assume that this is expensive.  Is this true and are
> there any standard procedures around it?

There is no way to explicitely allocate on the stack in any Caml
implementation including CSL. This is just too low level and must be a
compiler optimization, not a user's concern. Unfortunately no Caml
compiler (except Bigloo in some cases) works as to allocate data on
the stack. Anyway, allocating in the heap is not too expansive,
especially since we have a very fast GC in CSL.

>  (Typical example)
> 
>   val f : float array -> float array
> 
> will certainly be expensive, but a `subroutine'
> 
>   val g : float array -> unit
> 
> modifying the vector in place will also (in general) need a temporary
> array to store either input or output.

If your code does not explicitely allocate a new array in the body of
the procedure, CAML won't allocate any temporary array. There is no
semantic feature in the language that assumes that some arrays have to
be copied (as in Pascal or Ada for instance): arrays are always
efficiently handled as pointers to some memory locations.

> Another example: how can I find out (without deciphering assembler
> code) if a particular recursive function has been recognized as tail
> recursive?

You cannot. Once more this is too low level to be of programmer's
concern (in the Caml philosophy). However, our compilers are rather
good at optimizing these tail calls (not only recursive ones). So if
it is reasonably evident, your Caml compiler must do it.

Remember that Caml is a functionnal programming language: so function
definitions and function calls are very heavily used, and the compiler
optimizes them as much as it can. For instance, some Caml compilers
may systematically turn imperative loops into recursive definitions of
functions, since they are able to recognize loops into those recursive
definitions, and to compile them as such.

Concerning CSL specifically, Xavier may probably tell you more ...

I would like to end this answer on efficiency considerations by
telling a little true story: somebody complains once that the Caml
Light system failed to execute his code with a fatal "Out of
memory..." error message, after having allocated more than 200
Megabytes of data. The supposed ``bug'' of Caml Light was that the
size of useful data was supposed to be strictly smaller than 200
Megabytes. The conclusion of the user was: ``Well, there's no problem,
I'm going to rewrite all this stuff in C, and use a huge malloc at the
beginning'' :(

After a look at the source code of his example, it appears that the
program got an exponential complexity anyway, and would not answer
before years regardless of the size of the available memory (in fact
the program needed hundred of Gigabytes of memory to complete, so that
a ``huge'' malloc was hopeless).
After a bit of brain storming, changing the algorithm leads to a
program that runs within 100 words of memory, runs exponentially
faster, and gives the result after a few minutes :)

So be sure that you really need to know this kind of very low level
details (such as the ``tail calls'' treatment for instance). My
feeling is that, if efficiency is a problem, you have to change the
complexity of the algorithm, not to try to change the code generated by
the compiler!
And if you need to test and change your programs rapidly, the
readibility and expressiveness of Caml programs may help you...

Pierre Weis
----------------------------------------------------------------------------
WWW Home Page: http://pauillac.inria.fr/~weis
Projet Cristal
INRIA, BP 105, F-78153 Le Chesnay Cedex (France)
E-mail: Pierre.Weis@inria.fr
Telephone: +33 1 39 63 55 98
Fax: +33 1 39 63 53 30
----------------------------------------------------------------------------

Re: caml (special) light and numerics

[Thorsten Ohl]

>>>>> "Pierre" == Pierre Weis <Pierre.Weis@inria.fr> writes:

Pierre> If your code does not explicitely allocate a new array in the
Pierre> body of the procedure, CAML won't allocate any temporary
Pierre> array.

How can I get around Array.new in 

    val f : float array -> float array
    
    let f x =
      let n = Array.length x in
      let y = Array.new n 0.0 in
	for i = 0 to n do
	  for j = 0 to n do
	    y.(i) <- y.(i) +. a.(i).(j) *. x.(j)
	  done
	done

The problem is that I cannot modify `x' in place (even if I know that
I won't need it later).

The solution

    val g : float array -> float array -> unit
    
    let f x y =
      let n = Array.length x in
	for i = 0 to n do
	  for j = 0 to n do
	    y.(i) <- y.(i) +. a.(i).(j) *. x.(j)
	  done
	done

comes to mind, but it forces me to think of the operation in terms of
side efects, not in a functional way.

I was hoping to build a library in which I could use functors to build
different incarnations of an algorithm: a slow multi-dimentional one
using arrays and a fast one for one dimensional problems.

functor (In : Vectorspace, Out : Vectorspace) ->
  sig
    val f : In.type -> Out.type
    ...
  end

But it seems that I will have to use references for the
one-dimensional case too.

BTW: If this sound like sensless babble to you, please warn me that
     I'm wasting my time and should better buy a Fortran90 compiler ...

Pierre> However, our compilers are rather good at optimizing these
Pierre> tail calls (not only recursive ones). So if it is reasonably
Pierre> evident, your Caml compiler must do it.

I'll take your word for it.

Pierre> I would like to end this answer on efficiency considerations
Pierre> by telling a little true story: [...]  After a bit of brain
Pierre> storming, changing the algorithm leads to a program that runs
Pierre> within 100 words of memory, runs exponentially faster, and
Pierre> gives the result after a few minutes :)

This is the effect I'm looking for -- and I won't be satified with
anything less :-).

Pierre> My feeling is that, if efficiency is a problem, you
Pierre> have to change the complexity of the algorithm, not to try to
Pierre> change the code generated by the compiler!

Sure.  However: if we assume for simplicity that the compiler induces
a multiplicative overhead (wrt. low level languages like FORTRAN), and
I manage to reduce the complexity by one power

    FORTRAN:  C_F * O(n^2)
    CSL:      C_C * O(n)

then I have to estimate (roughly) C_C/C_F to find out where the break
even point is.  Otherwise I will not be able to convince anyone.

Pierre> And if you need to test and change your programs rapidly, the
Pierre> readibility and expressiveness of Caml programs may help
Pierre> you...

No doubt about that, but I'm trying to figure out if it is _likely_
that also production level code can be produced.  Cslopt is certainly
an example, but from a very different field.  I'm thinking about
implementing a non-trivial system and I need _some_ idea if it has a
chance to fly.  Flashes of genius that reduce an exponential solution
to a power law are nice, but rare in numerics.  OTOH, a factor of ten
can be painful if it means that you have to wait a week, instead of of
a night.

Thanks for your patience,
-Thorsten

/// Thorsten Ohl, TH Darmstadt, Schlossgartenstr. 9, D-64289 Darmstadt, Germany
//////////////// net: ohl@crunch.ikp.physik.th-darmstadt.de, ohl@gnu.ai.mit.edu
/// voice: +49-6151-16-3116, secretary: +49-6151-16-2072, fax: +49-6151-16-2421

Re: caml (special) light and numerics

[Pierre Weis]

> How can I get around Array.new in 
> 
>     val f : float array -> float array
>     
>     let f x =
>       let n = Array.length x in
>       let y = Array.new n 0.0 in
> 	for i = 0 to n do
> 	  for j = 0 to n do
> 	    y.(i) <- y.(i) +. a.(i).(j) *. x.(j)
> 	  done
> 	done
> 
> The problem is that I cannot modify `x' in place (even if I know that
> I won't need it later).

You're right, you must allocate it (but don't forget to return it at the
end of the loop).

> The solution
> 
>     val g : float array -> float array -> unit
>     
>     let f x y =
>       let n = Array.length x in
> 	for i = 0 to n do
> 	  for j = 0 to n do
> 	    y.(i) <- y.(i) +. a.(i).(j) *. x.(j)
> 	  done
> 	done
> 
> comes to mind, but it forces me to think of the operation in terms of
> side efects, not in a functional way.

That's true. Write it that way only if this computation is a
side-effect, that is if the y parameter is in fact an accumulator.

> No doubt about that, but I'm trying to figure out if it is _likely_
> that also production level code can be produced.  Cslopt is certainly
> an example, but from a very different field.  I'm thinking about
> implementing a non-trivial system and I need _some_ idea if it has a
> chance to fly.
> OTOH, a factor of ten can be painful if it means that you have to
> wait a week, instead of a night.

Well it's very difficult to answer. I think that it is worth
trying. Especially because Caml is easy to interface to C: if you get
a bottleneck somewhere, and if the corresponding C code is much more
efficient than the Caml code, you may delegate this computation to
C. On the other hand, the rest of your system may stay in Caml and be
written in a clean and easy to maintain style.

We got a similar architecture for the bignum arithmetic facility of
Caml: the very low level layer of the arithmetic that performs in place
computation is written in assembly code (or in C, depending on the
hardware). Higher level layers of the arithmetic were written in Caml:
for these parts, the algorithmic was not so well established and once
more we found it comfortable to experiment in Caml.

Pierre Weis
----------------------------------------------------------------------------
WWW Home Page: http://pauillac.inria.fr/~weis
Projet Cristal
INRIA, BP 105, F-78153 Le Chesnay Cedex (France)
E-mail: Pierre.Weis@inria.fr
Telephone: +33 1 39 63 55 98
Fax: +33 1 39 63 53 30
----------------------------------------------------------------------------

Re: caml (special) light and numerics

[Stefan Monnier]

> How can I get around Array.new in 
> 
>     val f : float array -> float array
>     
>     let f x =
>       let n = Array.length x in
>       let y = Array.new n 0.0 in
> 	for i = 0 to n do
> 	  for j = 0 to n do
> 	    y.(i) <- y.(i) +. a.(i).(j) *. x.(j)
> 	  done
> 	done
> 
> The problem is that I cannot modify `x' in place (even if I know that
> I won't need it later).

>From what I understand, you're concerned more about the fact that "new" might
be slow than about the fact that the code will require two arrays in memory at
the same time (hence having a bigger working set).

You should be aware of the fact that heap allocation is very frequent in caml
and is hence made fairly efficient. In the code above, "new" is very unlikely
to represent more than a few pathetic percents of the time spent in the function
(unless "n" is *very* small (like 0 or 1)).


	Stefan

Re: caml (special) light and numerics

[Xavier Leroy]

> This may be asking too much, but are there any (written) guidelines
> for efficient numerical CSL coding (besides the one on the WWW site
> that gave me some hints)?

As Pierre Weis said, don't worry too much about tail calls or heap
allocations: these are pretty efficient. Other factors like boxing of
floats (see below) will kill you much earlier if not taken into
account.

Here are what I think are the main points to remember for writing
efficient numerical code in Caml Special Light. I don't have much
experience with numerical code, in CSL or otherwise, so these are only
hints.

0- The cslopt compiler is poorly named because it basically does not
optimize the user's code. That is, everything will be executed pretty
much as written in the source. Function applications will always
perform a function call (no inlining). Function definitions fun x -> ...
will always heap-allocate a closure (no lambda-lifting, no closure
hoisting). Arithmetic expressions will be computed at the point where
they occur in the source (no loop-invariant code motion). So, the
programmer is very much in control. (Inlining might be added some day,
but probably no other optimization.)

1- Matrices in Caml are really arrays of (pointers to) arrays. So, a
matrix access a.(i).(j) reads the i-th element of a, which is a
pointer to an array, then read the j-th element of that array.
It's about as efficient as the standard layout for matrices when the
dimensions are not known at compile-time (we're trading a dereference
for a multiplication), but can be a big loss if the dimensions are
statically known and one happens to be, say, a power of 2.

A first consequence is that rows of a matrix are really individual
arrays: they need not be the same length (great for triangular matrices)
and a whole row of a matrix can be replaced in constant time: just do
a.(i) <- some_array. There might also be some creative things you can
do by sharing a row between several matrices.

Another consequence is that matrices are row-major (hope that's the
right term): the first coordinate should vary slower than the second.
E.g. to sum two N*N matrices, don't write

        for j = 0 to N-1 do
          for i = 0 to N-1 do
            c.(i).(j) <- a.(i).(j) + b.(i).(j)
          done
        done

but write:

        for i = 0 to N-1 do
          let row_a = a.(i) and row_b = b.(i) and row_c = c.(i) in
          for j = 0 to N-1 do
            row_c.(j) <- row_a.(j) + row_b.(j)
          done
        done

(Notice the application of point 0: three loop-invariant expressions
have been lifted manually outside of the inner loop.)

2- Floating-point boxing. Because of the combination of garbage
collection, polymorphism/type abstraction, and separate compilation,
floating-point numbers are often boxed, that is, stored in a
heap-allocated block and handled through a pointer. Most Lisp and ML
implementations box (or tag) floats, and that's the main reason why
they deliver poor floating-point performance: any arithmetic operation
on boxed floats has to load from memory its arguments, compute the
result, heap-allocate a block, and store the result in it.

One of the ancestors of Caml Special Light, the Gallium experimental
compiler, implemented a fairly aggressive type-based unboxing strategy
that would completely eliminate boxing from Fortran-style
code. i.e. code without polymorphic functions. Unfortunately, it
turned out to complicate greatly the runtime system and make the
garbage collector less efficient, hence hampering the performance of
symbolic code. Since symbolic computation is still the main
application of ML, I dropped this unboxing strategy.

However, some of the unboxing tricks developed for Gallium are still
used in Caml Special Light, at least on a local scale (where they
don't interact with the garbage collector and everything). Here are
some hints on which floats are boxed and which are not:

a. Function arguments are always boxed.

b. Free variables of functions are always boxed.

c. Intermediate results of arithmetic expressions are not boxed.
   By "intermediate result", I mean the result of an arithmetic
   operation (+. -. *. /. ** sin cos exp log ...) that is immediately
   used as argument of an arithmetic operation.

   Example: fun x -> 3.14 *. x *. x
   The parameter x and the result are boxed, but not x *. x.

d. "let"-bound results of arithmetic operations are not boxed if they
   are only used later as arguments to arithmetic operations.

   Example: let x = y +. z in x *. x               x not boxed
            let x = y +. z in x *. x *. f x        x boxed (passed as argument)

e. Floats in data structures are generally boxed. E.g. a list of
   floats is really a list of pointers to floats. But there are two
   exceptions:

   - arrays of floats (the floats are laid out contiguously, as in C;
                       this is not an array of pointers to floats)
   - record types whose fields are all floats (same as arrays)

   An access in a float array or in such a record avoids boxing in the
   same way as for arithmetic operations.

   Example: a.(i) <- 2.0 *. a.(i) performs only one float read and one
   float store, instead of two reads, a boxing and a store as in most ML
   implementations.

A few consequences of these hints:

a => inline small functions over floats

b + d => iterate with "for" and "while" loops instead of recursive functions
    Example:
        let x = y +. z in
        for i = 0 to N-1 do a.(i) <- x *. a.(i) done
    does not box x, but
        let x = y +. z in
        let rec iter i =
          if i >= N then () else (a.(i) <- x *. a.(i); iter(i+1))
        in iter 0
    is not only unreadable, but causes x to be boxed because it is
    free in iter.

e => use records in preference to tuples for representing points,
    complex numbers, etc. Example:
        type complex = {re: float; im: float}
    saves 4 loads and 2 allocations for each complex addition compared with
        type complex = float * float

all => Fortran/C style (big functions, loops all over the place) is
    compiled better than "classic" functional style (functions,
    functionals, iterators, combinators till you scream).

Another consequence of the optimization of float arrays is that there
are really two formats of arrays in the system: arrays of floats and
arrays of pointers/integers.

The illusion of a parametric array type 'a array is maintained at some
cost: when accessing an array whose static type is 'a array (e.g. in a
polymorphic function), a run-time test is generated to distinguish the
two array formats, and some additional boxing/unboxing may take place
if it's a float array. When more is known on the static type of the
array, no test is generated and no extra boxing takes place.

So, for maximal performance, don't operate over float arrays with
polymorphic array functions. For instance, transposing a float matrix with

        let transpose m =
          let m' = Array.new_matrix (Array.length m.(0)) (Array.length m) in
          for i = 0 to Array.length m do
            let row_i = m.(i) in
            for j = 0 to Array.length row_i do
              m'.(j).(i) <- row_i.(j)
            done
           done;
           m'

is much more efficient if you add a type constraint:

        let transpose (m: float array array) = ...

Well, that's all I can say in general. To know more, you'll have to
write some code and test it. Remember: write it cleanly, then profile
it, then optimize the hot spots -- not the other way around. And let us
know if you beat Fortran 90 in development time + running time.

One last thing:

> Another example: how can I find out (without deciphering assembler
> code) if a particular recursive function has been recognized as tail
> recursive?

Reading assembler code is always an option. You can also compile with
the "-dcmm" option. This will print one of the intermediate
representations, which shows pretty clearly tail calls vs. regular calls,
direct calls vs. indirect calls, and boxing/unboxing steps. That is,
if you can guess the semantics of that intermediate language, because
I will not answer any questions about it. Have fun.

- Xavier Leroy