/* Compilation line without OpenMP */
/*   $ gcc -Wall -lm -pedantic prog.c -o prog */

/* To compile and run an OpenMP program: */
/*   $ gcc -o program -fopenmp program.c */
/*   $ export OMP_NUM_THREADS=2 */
/*   $ ./program */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* #include <omp.h> */

typedef double** mat;

/* set matrix to zero matrix */
void mat_zero(mat x, int n) {
  int i, j;
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      x[i][j] = 0;
}

/* make a new matrix (look at how the pointers work) */
mat mat_new(int n) {
  int i;
  mat x = malloc(sizeof(double*) * n);
  x[0] = malloc(sizeof(double) * n * n);

  for (i = 0; i < n; i++)
    x[i] = x[0] + n * i;
  mat_zero(x, n);
  
  return x;
}

/* copy (from a linear vector) a matrix */
mat mat_copy(double *v, int n) {
  int i, j;
  mat x = mat_new(n);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      x[i][j] = v[i*n+j];
  return x;
}

/* free a matrix */
void mat_del(mat x) { free(x[0]); free(x); }

/* print a matrix */
void mat_show(mat x, char *fmt, int n) {
  int i, j;
  if (!fmt) fmt = "%2.2g";
  for (i = 0; i < n; i++) {
    printf(i ? "      " : " [ ");
    for (j = 0; j < n; j++) {
      printf(fmt, x[i][j]);
      printf(j < n - 1 ? "  " : i == n - 1 ? " ]\n" : "\n");
    }
  }
}

/* to check your programs */
mat mat_mul(mat a, mat b, int n) {
  int i, j, k;
  mat c = c = mat_new(n);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      for (k = 0; k < n; k++)
	c[i][j] += a[i][k] * b[k][j];
  return c;
}

/* the factorisation */
void mat_LU(mat A, mat L, mat U, int n) {
  int i, j, k;
  mat_zero(L, n);
  mat_zero(U, n);
  
  /* complete here!  */
  
}

/* returns a Walsh matrix of order n (should be of the form 2^k) */
mat walsh_matrix(int n)
{
  mat m, w;
  int i, j;
  if(n <= 1) {
    m = mat_new(1);
    m[0][0] = 1;
  } else {
    w = walsh_matrix(n / 2);
    m = mat_new(n);
    for(i = 0; i < n / 2; ++i) {
      for(j = 0; j < n / 2; ++j) {
        m[i][j]                 =   w[i][j];
        m[i][j + n / 2]         =   w[i][j];
        m[i + n / 2][j]         =   w[i][j];
        m[i + n / 2][j + n / 2] = - w[i][j];
      }
    }
    mat_del(w);
  }
  return m;
}

/* print non-zero values */
void boolprintmat(mat m, int n) {
  int i, j;
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      if(m[i][j] > 0.0001 || m[i][j] < -0.0001)
	printf("█");
      else
	printf(" ");
    }
    printf("\n");
  }
  printf("\n");
}

int main() {
  int n;
  mat A, L, U;
  
  /* should be a nice result */
  n = 1 << 4;
  L = mat_new(n);
  U = mat_new(n);
  A = walsh_matrix(n);
  mat_LU(A, L, U, n);
  boolprintmat(U, n);
  boolprintmat(L, n);
  mat_del(A);
  mat_del(L);
  mat_del(U);
  
  /* with a big matrix */
  if(0) {
    n = 1 << 10;
    L = mat_new(n);
    U = mat_new(n);
    A = walsh_matrix(n);
    mat_LU(A, L, U, n);
    mat_del(A);
    mat_del(L);
    mat_del(U);
  }
  return 0;
}