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

#include "temps.h"

static double **nouvelleMatrice(int n) {
  double **m =(double **) calloc(n,sizeof(double *));
  for ( int i = 0 ; i < n ; i++ )
    m[i] = (double *) calloc(n,sizeof(double));
  return(m);
}

static void liberationMatrice(double **m,int n) {
  for ( int i = 0 ; i < n ; i++ )
    free(m[i]);
  free(m);
}

static void initMatrice(double **m,int n,int seed) {
  srand(seed);
  for ( int i = 0 ; i < n ; i++ )
    for ( int j = 0 ; j < n ; j++ )
	  m[i][j] = (rand()%1000)/1000.0;
}

static void produitMatricesCarrees(double **m1,double **m2,double **r,int n) {
  for ( int i = 0 ; i < n ; i++ ) {
    for ( int j = 0 ; j < n ; j++ ) {
      double v = 0.0;
      for ( int k = 0 ; k < n ; k++ ) {
        v += m1[i][k] * m2[k][j]; }
      r[i][j] = v; } }
}

static void transposition(double **m,int n) {
  for ( int i = 0 ; i < n-1 ; i++ ) {
    for ( int j = i+1 ; j < n ; j++ ) {
      double aux = m[i][j];
      m[i][j] = m[j][i];
      m[j][i] = aux; } }
}

static void produitMatricesCarreesAvecTransposition(double **m1,double **m2,double **r,int n) {
  transposition(m2,n);
  for ( int i = 0 ; i < n ; i++ ) {
    for ( int j = 0 ; j < n ; j++ ) {
      double v = 0.0;
      for ( int k = 0 ; k < n ; k++ ) {
        v += m1[i][k]*m2[j][k]; }
      r[i][j] = v; } }
  transposition(m2,n);
}

static void testProduit(int sz,int n) {
  double ** m1 = nouvelleMatrice(sz);
  double ** m2 = nouvelleMatrice(sz);
  double ** m = nouvelleMatrice(sz);
  initMatrice(m1,sz,111);
  initMatrice(m2,sz,222);
  { long i = 0;
    debutChrono();
    do {
      if ( i%2 == 0 )
        produitMatricesCarrees(m1,m2,m,sz);
        else
        produitMatricesCarrees(m2,m1,m,sz);
      i++; }
    while ( i < n );
    finChrono();
    printf("%5d : %5d : %12lf %12lf",sz,n,tempsIntermediaire()/1000.0,tempsIntermediaire()*1000000.0/n/((double) sz*sz*sz)); }
  { long i = 0;
    debutChrono();
    do {
      if ( i%2 == 0 )
        produitMatricesCarreesAvecTransposition(m1,m2,m,sz);
        else
        produitMatricesCarreesAvecTransposition(m2,m1,m,sz);
      i++; }
    while ( i < n );
    finChrono();
    printf(" : %12lf %12lf\n",tempsIntermediaire()/1000.0,tempsIntermediaire()*1000000.0/n/((double) sz*sz*sz)); }
  liberationMatrice(m1,sz);
  liberationMatrice(m2,sz);
  liberationMatrice(m,sz);
}
  
static int testEgalite(double **m1,double **m2,int n) {
  for ( int i = 0 ; i < n ; i++ )
    for ( int j = 0 ; j < n ; j++ )
      if ( m1[i][j] != m2[i][j] )
        return false;
  return true;
}

void main(void) {
    //int n = 200;
    //double **m1 = nouvelleMatrice(200);
    //double **m2 = nouvelleMatrice(200);
    //double **p1 = nouvelleMatrice(200);
    //double **p2 = nouvelleMatrice(200);
    //initMatrice(m1,n,111);
    //initMatrice(m2,n,222);
    //produitMatricesCarrees(m1,m2,p1,n);
    //produitMatricesCarreesAvecTransposition(m1,m2,p2,n);
    //printf("%d\n",testEgalite(p1,p2,n));
    //liberationMatrice(m1,n);
    //liberationMatrice(m2,n);
    //liberationMatrice(p1,n);
    //liberationMatrice(p2,n);
  for ( int i = 400 ; i <= 600 ; i += 10 ) {
    double div = i/100.0;
    int sz =(int) (0.9+10000/div/div/div);
    testProduit(i,sz); }

  //testProduit(100,10000);
  //testProduit(200,1200);
  //testProduit(300,400);
  //testProduit(400,150);
  //testProduit(500,80);
  //testProduit(600,45);
  //testProduit(700,30);
  //testProduit(800,20);
  //testProduit(900,14);
  //testProduit(1000,10);
  //testProduit(1100,8);
  //testProduit(1200,6);
  //testProduit(1300,5);
  //testProduit(1400,4);
  //testProduit(1500,3);

  getchar();
}