Re: caml (special) light and numerics

[Thorsten Ohl <ohl@crunch.ikp.physik.th-darmstadt.de>]

>>>>> "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