Version française
Home     About     Download     Resources     Contact us    
Browse thread
[Caml-list] Q: Efficient operations on complex numbers
[ Home ] [ Index: by date | by threads ]
[ Search: ]

[ Message by date: previous | next ] [ Message in thread: previous | next ] [ Thread: previous | next ]
Date: -- (:)
From: Jan Kybic <kybic@f...>
Subject: [Caml-list] Q: Efficient operations on complex numbers
Hello,
        I would appreciate if you could give me same help on the
following topic:

Q: How to make efficient computation with complex values in OCaml?

Subquestions:

Q1: Can the boxing of arguments be avoided?

Q2: My ocamlopt (from Linux Mandrake ocaml-3.06-12mdk.i586.rpm) does
    not like the '-ffast-math' switch. Why?

Q2: Is OCaml the right language for me, is there any other 
    you would recommend for me to consider?

Background
==========

I want to develop a fairly complex numerical algorithm
involving hierarchical Boundary Element Method for the electromagnetic
problem, i.e. working with surface meshes, huge data, complicated data
structures, and computational kernel using complex numbers.
 
My first attempt at coding it was in C++ but I did not achieve to
program the full version of the algorithm, since I was simply getting
lost in details of the implementation, handling of the data etc. 

My second attempt was to recode the algorithm in Haskell using ghc. This went
fine and I was quite enjoying it. However, when I started to test the
program on realistic size problems. I found that unwanted laziness
kept creeping in that I could neither detect nor control. The program
consumed a lot of memory and was quite slow, even though the
computational kernel was fast enough. I improved
things a little by adding `seq` all over the place but I gave up in
the end.

My third attempt was to code again in C++, this time using a
third-party library. It went fine while I was doing simple things but
then what I could not make the library do what I needed and with the
limited documentation I had, I was not able to extend it.

So here I am, trying to do the whole thing again, this time in OCaml,
because I hope it would allow me to design the data
structures I will need and at the same time the resulting algorithm
should be sufficiently efficiently compiled, run fast and not used too
much memory.

My problem
==========

So far, I have implemented a part of the computational kernel in
OCaml that calculates approximations using spherical harmonics. It
uses complex values. The Ocaml version is basically is straightforward
rewrite of the C++ version and is almost exactly twice as slow as the
C++ version. This is intriguing because the parts code that only uses real
arithmetic is about as fast as the corresponding C++ version. So why
is working with complex numbers so expensive in OCaml and what can be
done about it?

Experiments
===========

1) I compiled the following trivial function (using ocamlopt on a Linux
box):

    let calc_sqrt _ = 
      let rec calc_sqrt' x i = 
        if i<=0 then x 
        else calc_sqrt' (0.5 *. ( x +. 2.0 /. x )) (i-1)
      in calc_sqrt' 2.0 10000000

It takes 0.4s while the corresponding C++ version 

    double calc_sqrt()
    {
      double x=2.0 ; 
      for (int i=0;i<10000000;i++)  x=0.5 * (x+2.0/x) ;
      return x ;
    }

takes 0.34s. That's fair enough, very result good for OCaml indeed.

2) I did the same thing, only with complex numbers from the Complex
module:

    let calc_sqrt_c _ = 
      let two = { Complex.re=2.0 ; Complex.im=0.0 } and 
          half = { Complex.re=0.5 ; Complex.im=0.0 }
      in
      let rec calc_sqrt_c' x i = 
        if i<=0 then x 
        else calc_sqrt_c' (Complex.mul half 
               ( Complex.add x (Complex.div two x ))) 
                             (i-1)
      in calc_sqrt_c' two 10000000

and in C++:

        
    typedef std::complex<double> complex_type ;

    complex_type calc_sqrt_c()
    {
      complex_type two(2.0,0), half(0.5,0) ; 
      complex_type x=two ;
      for (int i=0;i<10000000;i++)  x=half * (x+two/x) ;
      return x ;
    }


Here however, the OCaml version takes 1.25s against 0.56s in C++,
i.e. a 2-fold slowdown.

3) I could improve the OCaml time to 0.98s by including the code from
   the Complex module in my source file. (Is there a way of achieving
   the same result without actually copying it?) On the other hand, it
   did not matter whether I used tuples or records to represent the
   complex numbers, whether I defined the 'add' and 'mul' functions
   locally or at the top level, and whether I asked the compiler to
   inline or not. 

Speculation
===========

I deduce, that the main factor contributing to the slowdown, is the
boxing and unboxing of values. Can this be optimized away? Would MLton
do it? 

Conclusions
===========

I can live with 2-fold slowdown with respect to C++ if this permits me
to actually write the program. However, it seems to be a shame that
OCaml can generate code on par with C++ in some cases (real
arithmetic) but not in others. So if there is a way to do it, I would
like to know.

Thank you,
                Jan


-- 
-------------------------------------------------------------------------
Jan Kybic <kybic@ieee.org>                  tel. +420 2 2435 7264
       or <kybic@fel.cvut.cz>,     http://cmp.felk.cvut.cz/~kybic

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners