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


double random_number(double *rmin, double *rmax){
  return (*rmin) + (float)rand()/(float)RAND_MAX*((*rmax)-(*rmin));
}


int main(){
  double **A, **L, **U;
  double *Af, *Lf, *Uf;
  int i,j,k;
  double r_min,r_max;
  int n_ges;

  clock_t tcpu_1, tcpu_2;

  n_ges = 4;
  r_min = 0.0;
  r_max = 1.0;

  srand(time(0));

  A = (double**)malloc(n_ges*sizeof(double));
  L = (double**)malloc(n_ges*sizeof(double));
  U = (double**)malloc(n_ges*sizeof(double));

  Af = (double*)malloc(n_ges*n_ges*sizeof(double));
  Lf = (double*)malloc(n_ges*n_ges*sizeof(double));
  Uf = (double*)malloc(n_ges*n_ges*sizeof(double));

  for (i=0; i<n_ges; i++){
    A[i] = &Af[i*n_ges];
    L[i] = &Lf[i*n_ges];
    U[i] = &Uf[i*n_ges];
  }

  printf("\nA:\n"); 
  for(i=0; i<n_ges; i++){
    for(j=0; j<n_ges; j++){
      A[i][j] = random_number(&r_min,&r_max);
      L[i][j] = 0.;
      U[i][j] = 0.;
      printf("%6.3f",A[i][j]);
    }
    printf("\n");
  }

  for(k=0; k<n_ges-1; k++){
    for(i=k+1; i<n_ges; i++){
      L[i][k] = A[i][k]/A[k][k];
      for(j=k; j<n_ges; j++){
        A[i][j] = A[i][j] - L[i][k]*A[k][j];
      }
    }
  }

  for (i=0; i<n_ges; i++){
    L[i][i] = 1.;
    for (j=i; j<n_ges; j++){
      U[i][j] = A[i][j];
    }
  }

  printf("\nL:\n"); 
  for(i=0; i<n_ges; i++){
    for(j=0; j<n_ges; j++){
      printf("%6.3f",L[i][j]);
    }
    printf("\n");
  }

  printf("\nU:\n"); 
  for(i=0; i<n_ges; i++){
    for(j=0; j<n_ges; j++){
      printf("%6.3f",U[i][j]);
    }
    printf("\n");
  }

  printf("\nA:\n"); 
  for(i=0; i<n_ges; i++){
    for(j=0; j<n_ges; j++){
      A[i][j] = 0.;
      for (k=0; k<n_ges; k++){
        A[i][j] += L[i][k]*U[k][j];
      }
      printf("%6.3f",A[i][j]);
    }
    printf("\n");
  }

  return 0;
}