#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <pthread.h>
#ifdef XCCMEM
#include "xccmem.h"
#endif
#include <sched.h>

#define SIZE 1600
#define LINELENGTH 128
#define MAXPROC 128

/* 行列に型名を付ける */
typedef double mat_t[SIZE][SIZE];

/* LU分解関数 */
typedef void (*lu_func) (int P, int myid, int n, mat_t a);

/* アルゴリズム整理用 */
struct lu_algorithm{
  int id;
  lu_func func;
  char *desc;
};

/* 並列処理へのパラメータ */
struct workload {
  lu_func f;
  int P;
  int myid;
  int n;
  void *a;
};

// algorithms
void lu_rightlooking(int P, int myid, int n, mat_t a);
void lu_rightlooking_p(int P, int myid, int n, mat_t a);

// misc
double elapsed_time(struct timeval tp[2]);
char *readline(FILE *fp);
void p_helper(lu_func f, int P, int n, mat_t a);
void genmat(int n, mat_t a);
void printmat(int n, mat_t a);
int inputmat(mat_t a, FILE *fp);
void outputmat(int n, mat_t a, FILE *fp);
void copy_l(int n, mat_t a, mat_t l);
void copy_u(int n, mat_t a, mat_t u);
void transpose(int n, mat_t a);
void matmul(int n, mat_t a, mat_t b, mat_t c);
int comparemat(int n, mat_t a, mat_t b);
void usage(char *filename);

/* 同期処理用 */
#ifdef XCCMEM
int ba[MAXPROC];
#else
pthread_mutex_t mut = PTHREAD_MUTEX_INITIALIZER;
pthread_cond_t cond = PTHREAD_COND_INITIALIZER;
int progress = 0;
#endif

/* グローバル変数: 行列 */
mat_t a;
mat_t l, u, a_origin;

/* アルゴリズムリスト */
struct lu_algorithm algorithm[] =
{{1, lu_rightlooking,      "Right-looking"},
 {101, lu_rightlooking_p,  "Right-looking Parallel"},
};

/* main */
int main(int argc, char *argv[])
{
  struct timeval tp[2];
  int n = 10, al = 0, debug = 0, mode_in = 0, mode_out = 0, p = 0;
  int i, j;
  int silent = 0;
  double elapsed;
  FILE *fin, *fout;
  char *input = NULL, *output = NULL;
  struct lu_algorithm *lu;

  /* コマンドライン引数の解釈 */
  for(i = 1; i < argc; i++){
    if(strcmp(argv[i], "-i") == 0){
      /* 行列の入力元が指定される */
      i++;
      if(i == argc){
        printf("Specify the filename to read: %s\n", argv[i - 1]);
        return 0;
      }
      input = argv[i];
      if(strcmp(input, "-") == 0){
        /* 標準入力から */
        mode_in = 1;
      }else{
        /* ファイルから */
        mode_in = 2;
      }
    }else if(strcmp(argv[i], "-o") == 0){
      /* 行列の出力先が指定される */
      i++;
      if(i == argc){
        printf("Specify the filename to write: %s\n", argv[i - 1]);
        return 0;
      }
      output = argv[i];
      if(strcmp(output, "-") == 0){
        /* 標準出力へ */
        mode_out = 1;
      }else{
        /* ファイルへ */
        mode_out = 2;
      }
    }else if(strcmp(argv[i], "-d") == 0){
      /* 検算結果を表示する */
      debug = 1;
    }else if(strcmp(argv[i], "-a") == 0){
      /* 分解アルゴリズムが指定される */
      i++;
      if(i == argc){
        printf("Specify the algorithm to use: %s\n", argv[i - 1]);
        return 0;
      }
      sscanf(argv[i], "%d", &al);
    }else if(strcmp(argv[i], "-p") == 0){
      /* スレッド数が指定される */
      i++;
      if(i == argc){
        printf("Specify the number of threads to use: %s\n", argv[i - 1]);
        return 0;
      }
      sscanf(argv[i], "%d", &p);
      if(p < 1){
        printf("The number of threads must be 1 - %d: regarded as 1\n", MAXPROC);
        p = 1;
      }else if(p > MAXPROC){
        printf("The number of threads must be 1 - %d: regarded as %d\n", MAXPROC, MAXPROC);
        p = MAXPROC;
      }
    }else if(argv[i][0] >= '0' && argv[i][0] <= '9'){
      /* 行列のサイズが指定される */
      sscanf(argv[i], "%d", &n);
      if(n < 1 || n > SIZE - 4){
        printf("The size must be 1 - %d\n", SIZE - 4);
        return 0;
      }
    }else if(strcmp(argv[i], "-h") == 0){
      /* ヘルプ */
      usage(argv[0]);
      return 0;
    }else if(strcmp(argv[i], "-s") == 0){
      /* 経過時間のみ表示する */
      silent = 1;
    }else{
      /* 不明なパラメータ */
      printf("Invalid Parameter: %s\n", argv[i]);
      return 0;
    }
  }
  
  /* 既定のプロセッサ数とアルゴリズム */
  if(al == 0){
    /* アルゴリズム未指定 */
    if(p <= 1){
      /* 逐次版の中で速いものを選んでおく */
      al = 1;
    }else{
      /* 並列版の中で速いものを選んでおく */
      al = 101;
    }
  }else{
    /* アルゴリズム指定済み */
    if(p == 0){
      /* スレッド数未指定 */
      if(al >= 101){
        /* 並列版アルゴリズム指定済み */
        /* デフォルトは P = 2 */
        p = 2;
      }
    }
  }

  /* LU分解アルゴリズムを指定する */
  lu = NULL;
  for(i = 0; i < sizeof(algorithm) / sizeof(struct lu_algorithm); i++){
    if(algorithm[i].id == al){
      lu = &algorithm[i];
      if(!silent) printf("*** %s ***\n", lu->desc);
      break;
    }
  }
  if(lu == NULL){
    printf("Invalid Parameter: -a %d\n", al);
    return 0;
  }

  /* データを入力する */
  switch(mode_in){
  case 0:
    /* 生成関数を用いる */
    genmat(n, a);
    break;
  case 1:
    /* 標準入力から読み込む */
    n = inputmat(a, stdin);
    break;
  case 2:
    /* ファイルから読み込む */
    fin = fopen(input, "r");
    n = inputmat(a, fin);
    fclose(fin);
    break;
  }
  
  /* テスト行列の内容確認 */
  if(debug){
    /* 後で比較するために退避する */
    for(i = 0; i < n; i++){
      for(j = 0; j < n; j++){
        a_origin[i][j] = a[i][j];
      }
    }
  }

  /* 時間測定開始 */
  gettimeofday(tp, 0);

  /* LU分解 */
  if(al < 100){
    /* 逐次(最初の2引数はダミー) */
    (lu->func)(1, 0, n, a);
  }else{
    /* 並列  */
    p_helper(lu->func, p, n, a);
  }

  /* 時間測定終了, 経過時間表示 */
  gettimeofday(tp+1, 0);
  elapsed = elapsed_time(tp);

  if(!silent){
    printf("Elapsed Time: %f sec.\n", elapsed);
  }else{
    printf("%f\n", elapsed);
  }

  /* LUを取り出しておく */
  copy_l(n, a, l);
  copy_u(n, a, u);

  /* 結果を出力する */
  switch(mode_out){
  case 0:
    /* 何もしない */
    break;
  case 1:
    /* 標準出力へ書き込む */
    puts("\nOutput L:");
    outputmat(n, l, stdout);
    puts("\nOutput U:");
    outputmat(n, u, stdout);
    break;
  case 2:
    /* ファイルへ書き込む */
    fout = fopen(output, "w");
    fputs("Output L:\n", fout);
    outputmat(n, l, fout);
    fputs("\nOutput U:\n", fout);
    outputmat(n, u, fout);
    fclose(fout);
    break;
  }

  /* 検算処理 */
  if(debug){
    matmul(n, l, u, a);
    /* 元の行列と等しいかどうかを調べる */
    if(comparemat(n, a, a_origin)){
      /* 等しい */
      puts("Valid.");
    }else{
      /* 等しくない */
      puts("Invalid.");        
    }
  }

  return 0;
}

/* 経過時間の計算 */
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);
}

/* バリア同期 */
#ifdef XCCMEM
/* 共有メモリ用 */
void init_barrier_sync(int P){
  int i;
  for(i = 0; i < P; i++) ba[i] = 1;
}
 
void barrier_sync(int P, int myid){
  if (myid == 0){
    int i, b = -ba[0];
    for(i = 1; i < P; i++)
      do{
        while(xread_int(ba[i]) != b);
      }while(atomic_read_int(ba[i]) != b);
    start_access_after_read();
    atomic_write_int_to_finish_access(ba[0],b);
  }else{
    int b = -ba[myid];
    atomic_write_int_to_finish_access(ba[myid], b);
    do{
      while(xread_int(ba[0]) != b);
    }while(atomic_read_int_to_start_access(ba[0]) != b);
  }
}
#else
/* pthreadの同期制御 */
void barrier_sync(int P, int myid)
{
  int prg;
  pthread_mutex_lock(&mut);
  if((prg = progress + 1) < P)
  {
    progress = prg;
    pthread_cond_wait(&cond, &mut);
  }
  else
  {
    progress = 0;
    pthread_cond_broadcast(&cond);
  }
  pthread_mutex_unlock(&mut);
}
#endif

/* スレッド生成 */
int systhr_create(void * (*start_func)(void *), void *arg){
  int status = 0;
  pthread_t tid;
  pthread_attr_t attr;
  pthread_attr_init(&attr);
  status = pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
  if(status == 0) status = pthread_create(&tid, &attr, start_func, arg);
  if(status != 0) status = pthread_create(&tid, 0,     start_func, arg);
  return status;
}

/* 使用するCPUを指定する */
void assign_cpu(int cpu){
  cpu_set_t cpu_set;
  /* 指定したCPUのみを集合に加える */
  CPU_ZERO(&cpu_set);
  CPU_SET(cpu, &cpu_set);
  /* 割り当て */
  sched_setaffinity(0, sizeof(cpu_set_t), &cpu_set);
}

/* 1行を読み込んだメモリ領域を返す */
char *readline(FILE *fp){
  char *buf;
  int len = LINELENGTH;
  buf = (char *)malloc(sizeof(char) * len);
  buf[0] = '\0';
  while(fgets(buf + strlen(buf), len - strlen(buf), fp)){
    /* 改行文字以外で終わっていれば終了する */
    if(buf[strlen(buf) - 1] == '\n') break;
    /* 続行する */
    len += LINELENGTH;
    buf = (char *)realloc(buf, len);
  }
  return buf;
}

/* 並列処理の起動 */
void *exec_workload(void *w0)
{
  struct workload *w = (struct workload *)w0;
  /* CPUの指定 */
  assign_cpu(w->myid);
  /* メイン処理 */
  (w->f)(w->P, w->myid, w->n, w->a);
  return NULL;
}

/* 並列処理の準備 */
void
p_helper(lu_func f, int P, int n, mat_t a)
{
  int i;
  struct workload *w = (struct workload *)malloc(sizeof(struct workload)*P);

  /* 同期処理の初期化 */
#ifdef XCCMEM
  init_barrier_sync(P);
#else
  progress = 0;
#endif
  
  /* 各並列処理に与えるパラメータ */
  for(i=0;i<P;i++){
    w[i].f = f;
    w[i].P = P;
    w[i].myid = i;
    w[i].n = n;
    w[i].a = a;
  }
  /* 処理を開始する */
  for(i=1;i<P;i++) systhr_create(exec_workload, w+i);
  exec_workload(w);
}

/* 行列生成関数 */
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++){
      fprintf(stdout, "%6.3f%c", a[i][j], (j==n-1) ? '\n' : ' ');
    }
  }
  fputc('\n', stdout);
}

/* 行列入力関数 返却値は行列のサイズ */
int
inputmat(mat_t a, FILE *fp)
{
  char *line, *p;
  int n;
  int i, j;
  /* サイズの入力 */
  line = readline(fp);
  sscanf(line, "%d", &n);
  free(line);

  /* 行列の成分の入力 */
  for(i = 0; i < n; i++){
    /* 行の入力 */
    line = readline(fp);
    /* 成分ごとに分割する */
    p = strtok(line, " ");
    for(j = 0;j < n; j++){
      if(p == NULL){
        /* 成分が不足している */
        fputs("Too few elements given. Retry!\n", stderr);
        break;
      }
      sscanf(p, "%lf", &a[i][j]);
      p = strtok(NULL, " ");
    }
    free(line);
    /* やり直し */
    if(j < n) i--;
  }
  return n;
}

/* 行列出力関数 返却値は行列のサイズ */
void
outputmat(int n, mat_t a, FILE *fp)
{
  int i,j;
  /* サイズの出力 */
  fprintf(fp, "%d\n", n);

  /* 行列の成分の出力 */
  for(i=0;i<n;i++){
    for(j=0;j<n;j++){
      fprintf(fp, "%6.3f%c", a[i][j], (j==n-1) ? '\n' : ' ');
    }
  }
}

/* Lの抽出(Crout型) */
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];
    }
  }
}

/* Uの抽出(Crout型) */
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=32;
  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;
        }
      }
    }
  }
}

/* 行列の比較 */
int comparemat(int n, mat_t a, mat_t b){
  int i, j;
  double err;
  for(i = 0; i < n; i++){
    for(j = 0; j < n; j++){
      /* 誤差を許容して各成分を比較する */
      err = a[i][j] - b[i][j];
      if(!(err > -1e-8 && err < 1e-8)){
        /* 等しくないとみなす */
        return 0;
      }
    }
  }
  return 1;
}

/* 利用方法 */
void usage(char *filename){
  int i;

  printf("Usage: %s [matrix_size] [options...]\n", filename);
  puts("Options:");
  puts("  -a <number>: choose the algorithm to use");
  for(i = 0; i < sizeof(algorithm) / sizeof(struct lu_algorithm); i++){
    printf("     %3d: %s\n", algorithm[i].id, algorithm[i].desc);
  }
  puts("  -h: show this help");
  puts("  -i -: input matrix from standard input ([matrix_size] is ignored)");
  puts("  -i <filename>: input matrix from <filename> ([matrix_size] is ignored)");
  puts("  -o -: output matrix to standard input");
  puts("  -o <filename>: output matrix to <filename>");
  puts("  -d: check the validity of the result");
  puts("  -s: only show elapsed time (if no error occured)");
  puts("  -p <number>: the number of threads to use");
}


/********** Implementations of algorithms *******************/

/* right-looking [-a 1] : 基本アルゴリズム */
void
lu_rightlooking(int P, int myid, 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];
      }
    }
  }
}


/* LU分解の実行(right-looking) 並列版 [-a 101] */
void
lu_rightlooking_p(int P, int myid, int n, mat_t a)
{
  int i,j,k;
  for(k=0;k<n-1;k++){
    /* 分担 */
    {
      /* 通常入口では同期しなくてよい。ここでは同期2で十分 */
      /* barrier_sync(P, myid); */ 
      if (myid == 0) {
        /* メインスレッドだけが実行 */
        for(j=k+1;j<n;j++){
          a[k][j] /= a[k][k];
        }
      } else {
        /* 他のスレッドは何もしないが同期1はする */
      }
      /* 同期1 */
      barrier_sync(P, myid);
    }
    /* 分担 */
    {
      /* 担当範囲 */
      int i0 = (k+1) + (n - (k + 1)) * myid / P;
      int i1 = (k+1) + (n - (k + 1)) * (myid+1) / P;
      /* 通常入口では同期しなくてよい。ここでは同期1で十分 */
      /* barrier_sync(P, myid); */
      for(i=i0;i<i1;i++){
        for(j=k+1;j<n;j++){
          a[i][j] -= a[i][k] * a[k][j];
        }
      }
      /* 同期2 */
      barrier_sync(P, myid);
    }
  }
  /* ここでは同期2で十分 */
  /* barrier_sync(P, myid); */
}