#include <stdio.h>
#include <sys/time.h>

#define SIZE 1600

typedef double mat_t[SIZE][SIZE];
typedef double vec_t[SIZE];

double
elapsed_time(struct timeval tp[2])
{
  return tp[1].tv_sec-tp[0].tv_sec+1e-6*(tp[1].tv_usec-tp[0].tv_usec);
}

void
lu_f(int n, mat_t a)
{
  int i,j,k;
  for(k=0;k<n-1;k++)
    {
      for(j=k+1;j<n;j++)
        a[k][j] /= a[k][k];
      for(i=k+1;i<n;i++)
        for(j=k+1;j<n;j++) a[i][j] -= a[i][k] * a[k][j];
    }
}

void
genmat(int n, mat_t a)
{
  int i,j;
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      a[i][j] = n - abs(i-j);
}

void
printmat(int n, mat_t a)
{
  int i,j;
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.3lf%c", a[i][j], (j==n-1) ? '\n' : ' ');
  putchar('\n');
}

void
copy_l(int n, mat_t a, mat_t l)
{
  int i,j;
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      l[i][j] = (i<j) ? 0.0 : a[i][j];
}

void
copy_u(int n, mat_t a, mat_t u)
{
  int i,j;
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      u[i][j] = (i==j) ? 1.0 : (i>j) ? 0.0 : a[i][j];
}

void
transpose(int n, mat_t a)
{
  int i,j;
  for(i=0;i<n;i++) 
    for(j=0;j<i;j++)
      {
	double t1 = a[i][j], t2 = a[j][i];
	a[i][j] = t2, a[j][i] = t1;
      }
}

void
matmul(int n, mat_t a, mat_t b, mat_t c)
{
  int bs=12;
  int i0,j0,i,j,k;
  transpose(n, b);
  for(i0=0;i0<n;i0+=bs)
    for(j0=0;j0<n;j0+=bs)
      for(i=i0;i<i0+bs && i<n;i++)
	for(j=j0;j<j0+bs && j<n;j++)
	  {
	    double s = 0.0;
	    for(k=0;k<n;k++) s += a[i][k] * b[j][k];
	    c[i][j] = s;
	  }
}

mat_t a, l, u;

main(int argc, char *argv[])
{
  struct timeval tp[2];
  int n = 10, al = 1, d = 0;
  if(argc > 1) n = atol(argv[1]);
  if(argc > 2) al = atol(argv[2]);
  if(argc > 3) d = 1;
  genmat(n, a);
  if(d) printmat(n, a);
  gettimeofday(tp, 0);
  switch(al){
  case 1: lu_f(n, a); break;
  }
  gettimeofday(tp+1, 0);
  printf("%lf\n", elapsed_time(tp));

  if(d){
    printmat(n, a);
    copy_l(n, a, l);
    copy_u(n, a, u);
    printmat(n, l);
    printmat(n, u);
    matmul(n, l, u, a);
    printmat(n, a);
  }
}