Multiplication of matrix in C and OCaml

["Frédéric Gava" <gava@univ-paris12.fr>]

Dear All,

I compared multiplication of matrix in C and OCaml and I was a little 
surprise to see that the following C code (using -O2) is 8 time faster 
than the OCaml one (even with -unsafe).

Anybody have an idea to optimize my OCaml code or know why is there this 
"big" difference ?

many thanks

Frédéric Gava

ps: the difference is the same even when I use a non-polymorphic 
multiplication

*************** C Code ******************

typedef struct Complexe {double re; double img;} complexe;

__inline__ void* add(void *x, void *y)
{
complexe a=((complexe*) x)[0];
complexe b=((complexe*) y)[0];
return (void *) &((complexe){a.re+b.re, a.img+b.img});
}

__inline__ void* mult(void *x, void *y)
{
complexe a=((complexe*) x)[0];
complexe b=((complexe*) y)[0];
return (void *) &((complexe){a.re*b.re-a.img*b.img, a.re*b.img+a.img*b.re});
}

void multiply_generic(int n, void *zero, void* (*add)(void *a, void *b), 
void* (*mult)(void *a, void *b), void *a, void *b, void *c)
{
  register int i,j,k;
  void *tmp;
  void *cc;

   for (i=0; i<n; i++)
    {
     for (j=0;j<n;j++)
      {
       tmp=zero;
       for (k=0;k<n;k++) tmp=add(tmp,(mult(b+(k*n+j),a+(i*n+k))));
       cc=c+(i*n+j);
       cc=tmp;
      }
   }
}

complexe random_complexe(void)
{
  complexe c = {c.re=(double)(rand()%1000), c.img=(double)(rand()%1000)};
  return c;
}


Call this with (and with random matrix).

complexe zero;
zero.re=(double)0;
zero.img=(double)0;
multiply_generic(i,&zero,add,mult,a,b,c);


********************* OCaml Code *******************************

let random_complex n = {Complex.re=Random.float n; im=Random.float n}

let matrixRandom n = Array.init (n*n) (fun _ -> random_complex 1000.)

let multiplication n zero add mul a b c =
   let tmp=ref zero in
   for i=0 to n-1 do
    for j=0 to n-1 do
     tmp:=zero;
     for k=0 to n-1 do
      tmp:= add (!tmp) (mul b.(k*n+j) a.(i*n+k))
     done;
     c.(i*n+j)<-(!tmp)
    done
   done