/* param 1 = size of the matrices at the beginning
   param 2 = incrementation
   param 3 = size of the matrices at the end
   param 4 = number of iteration 

with 
gcc -O3 -o cmult mult.c

*/

#include <stdlib.h> 
#include <stdio.h> 
#include <math.h>

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

complexe zero;

__inline__ complexe complexe_add(complexe a, complexe b)
{
 return (complexe){a.re+b.re, a.img+b.img};
}


__inline__ complexe complexe_mult(complexe x, complexe y)
{
 return (complexe) {x.re*y.re-x.img*y.img, x.re*y.img+x.img*y.re};
}


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


void multiply_complex(int n, complexe *a, complexe *b, complexe *c)
{
 register int i,j,k;
 complexe tmp;

  for (i=0; i<n; i++) 
   {
    for (j=0;j<n;j++)
     {
      tmp=zero;
      for (k=0;k<n;k++) 
       tmp=complexe_add(tmp,(complexe_mult(b[k*n+j],a[i*n+k])));
      c[i*n+j]=tmp;
     }
  }
 return;
}

void multiply_complex_generic(int n, complexe (*add)(complexe x, complexe y) , complexe (*mult)(complexe x, complexe y),complexe *a, complexe *b, complexe *c)
{
 register int i,j,k;
 complexe tmp;

  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])));
      c[i*n+j]=tmp;
     }
  }
 return;
}


int main(int argc, char* argv[])
{
 complexe *a, *b, *c;
 int i,j,k,t1,pas,maxi,iter;

 zero.re=(double)0;
 zero.img=(double)0;
 
 i=atoi(argv[1]);
 pas=atoi(argv[2]);
 maxi=atoi(argv[3]);
 iter=atoi(argv[4]);

 srand(1000);
 while (i<maxi)
  {  
   /* initialisation des matrices */
   a=(complexe *)malloc(i*i*sizeof(complexe));
   b=(complexe *)malloc(i*i*sizeof(complexe));
   c=(complexe *)malloc(i*i*sizeof(complexe));
      
   for (j=0;j<i;j++)
    for (k=0,t1=j*i;k<i;k++)
     {
      a[t1+k]=random_complexe();
      b[t1+k]=random_complexe();
     }

  for (k=1;k<=iter;k++) 
   multiply_complex(i,a,b,c);
   /* multiply_complex_generic(i,complexe_add,complexe_mult,a,b,c); */

      

  /* suppression des vieilles matrices */
  free((complexe *)a);
  free((complexe *)b);
  free((complexe *)c);

  i=i+pas;
  }
 return 0;
}