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

#define N 1000

double A[N][N], b[N];

int p;

int mdone, nwdone;
pthread_cond_t ecv = PTHREAD_COND_INITIALIZER;
pthread_mutex_t emtx = PTHREAD_MUTEX_INITIALIZER;
pthread_cond_t scv = PTHREAD_COND_INITIALIZER;
pthread_mutex_t smtx = PTHREAD_MUTEX_INITIALIZER;

typedef struct {
   int startrow;
   int endrow;
} range;

range rows[16];

void geneqns ( )
{
   int i, j, t, x[N], rhs;

   for (i=0; i<N; ++i) x[i] = rand() % 19 - 9;
   for (i=0; i<N; ++i) {
      rhs = 0;
      for (j=0; j<N; ++j) {
         t = rand() % 19 - 9;
         A[i][j] = (double)t;
         rhs += t * x[j];
      }
      b[i] = (double)rhs;
   }
}

void rowredinit ( )
{
   int i, j, k;
   int nr, nrpf, nrpc, nrpr;
   double t;

   for (i=0; i<N; ++i) {
      for (k=i; k<N; ++k) if (A[k][i]) break;
      if (k == N) {
         printf("\n*** Master thread detected underspecified system. Aborting...\n");
         exit(0);
      }
      if (k > i) {
         for (j=i; j<N; ++j) {
            t = A[i][j]; A[i][j] = A[k][j]; A[k][j] = t;
         }
         t = b[i]; b[i] = b[k]; b[k] = t;
      }
      t = A[i][i];
      for (j=i; j<N; ++j) A[i][j] /= t;
      b[i] /= t;

      nr = N-i-1; nrpf = nr / p; nrpc = (nr + (p-1)) / p; nrpr = nr % p;
      j = i+1;
      for (k=0; k<p; ++k) {
         rows[k].startrow = j;
         j += (k < nrpr) ? nrpc : nrpf;
         rows[k].endrow = j-1;
      }

      pthread_mutex_lock(&smtx);
      nwdone = 0; mdone = 1;
      pthread_cond_broadcast(&scv);
      pthread_mutex_unlock(&smtx);

      pthread_mutex_lock(&emtx);
      if (nwdone < p) pthread_cond_wait(&ecv,&emtx);
      pthread_mutex_unlock(&emtx);
   }
}

void *rowred ( void *arg )
{
   int tno, i, j, k;
   double t;

   tno = *((int *)arg);
   printf("+++ Worker thread %d starting\n", tno);

   for (i=0; i<N; ++i) {
      pthread_mutex_lock(&smtx);
      if (!mdone) pthread_cond_wait(&scv,&smtx);
      pthread_mutex_unlock(&smtx);

      /* Do the row reduction in parallel */
      for (k=rows[tno].startrow; k<=rows[tno].endrow; ++k) {
         if ((t = A[k][i])) {
            for (j=i; j<N; ++j) A[k][j] -= t * A[i][j];
            b[k] -= t * b[i];
         }
      }

      pthread_mutex_lock(&emtx);
      ++nwdone;
      if (nwdone == 1) mdone = 0;
      if (nwdone == p) pthread_cond_signal(&ecv);
      pthread_mutex_unlock(&emtx);
   }

   printf("+++ Worker thread %d exiting\n", tno);
   pthread_exit(NULL);
}

void gensol ( double x[], double A[][N], double b[] )
{
   int  i, j;

   for (i=N-1; i>=0; --i) {
      x[i] = b[i];
      for (j=i+1; j<N; ++j) x[i] -= A[i][j] * x[j];
   }
}

double verifysol ( double A[][N], double b[], double x[] )
{
   double maxerr, rowerr;
   int i, j;

   for (i=0; i<N; ++i) {
      rowerr = 0;
      for (j=0; j<N; ++j) rowerr += A[i][j] * x[j];
      rowerr -= b[i];
      if (rowerr > maxerr) maxerr = rowerr;
   }
   return maxerr;
}

int main ( int argc, char *argv[] )
{
   int tno[16], i, j;
   pthread_t tid[16];
   double Acopy[N][N], bcopy[N];
   double x[N], maxerr;

   p = (argc > 1) ? atoi(argv[1]) : 8;
   if ((p < 2) || (p > 16)) p = 8;
   printf("+++ p = %d\n\n", p);

   srand((unsigned int)time(NULL));

   geneqns();
   for (i=0; i<N; ++i) for (j=0; j<N; ++j) Acopy[i][j] = A[i][j];
   for (i=0; i<N; ++i) bcopy[i] = b[i];

   for (i=0; i<p; ++i) tno[i] = i;

   pthread_mutex_trylock(&smtx);
   pthread_mutex_unlock(&smtx);
   pthread_mutex_trylock(&emtx);
   pthread_mutex_unlock(&emtx);

   mdone = nwdone = 0;

   for (i=0; i<p; ++i)
      pthread_create(tid+i,NULL,rowred,(void *)(tno+i));

   /* Wait a little bit for the worker threads to start working */
   usleep(1000);

   printf("\n");

   rowredinit();

   for (i=0; i<p; ++i) pthread_join(tid[i],NULL);

   printf("\n+++ Master thread: Row reduction done\n");

   /* for (i=0; i<N; ++i) {
      for (j=0; j<N; ++j) printf("%lf ", A[i][j]);
      printf("\n");
   }
   for (i=0; i<N; ++i) printf("%lf ", b[i]);
   printf("\n"); */

   gensol(x,A,b);
   printf("\n+++ Master thread: Solution computed\n");

   maxerr = verifysol(Acopy,bcopy,x);
   printf("    Maximum error (initial system) = %16.14lf\n", maxerr);

   maxerr = verifysol(A,b,x);
   printf("    Maximum error (reduced system) = %16.14lf\n", maxerr);

   /* for (i=0; i<N; ++i) printf("%lf ", x[i]);
   printf("\n"); */

   exit(0);
}