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

#define ABS(x) (x < 0 ? -(x) : (x))
#define EPS 0.00001
#define TRUE  1
#define FALSE 0
#define DEBUG 1

shared[1001] double a[1000][1001];
shared double x[1000];
shared int size;
shared int maxrow;

/*
    Solve a system of n equations in n unknowns using Gaussian Elimination
    Solve an equation in matrix form Ax = b
    The 2D array a is the matrix A with an additional column b.
    This is often written (A:b)

    A0,0     A1,0     A2,0     ....  An-1,0      b0
    A0,1     A1,1     A2,1     ....  An-1,1      b1
    A0,2     A1,2     A2,2     ....  An-1,2      b2
    :         :         :                 :             :
    :         :         :                 :             :
    A0,n-1  A1,n-1  A2,n-1  ....  An-1,n-1    bn-1

    The result is returned in x, otherwise the function returns FALSE
    if the system of equations is singular.
*/

void WriteSolution()
{
    int j,k;
    for (j=0;j<size;j++) {
        for (k=0;k<size+1;k++) {
            printf("%10.3f ",a[j][k]);
        }
        printf(" | %10.3f\n",x[j]);
    }
    printf("\n");
}

void WriteSolution2()
{
    int j,k;

    for (j=0;j<size;j++) {
        printf("%10.3f\n",x[j]);
    }
    printf("\n");
}



int GSolve()
{
    int i,j,k;
    double tmpa,tmpb;
    int b = 0;
    for (i=0;i<size;i++) {
        /* Find the row with the largest first value */
        
        if(MYTHREAD == 0){
            maxrow = i;
            for (j=i+1;j<size;j++) {
//                printf("I AM ROOT : j = %d, maxrow : %d\n",j,maxrow); fflush(stdout);
//                printf("a[%d][%d] : %f\n",j,i,a[j][i]); fflush(stdout);
//                printf("a[%d][%d] : %f\n",maxrow,i,a[maxrow][i]); fflush(stdout);
//                printf("ABS(a[%d][%d]) : %f\n",j,i,ABS(a[j][i])); fflush(stdout);
//                printf("ABS(a[%d][%d]) : %f\n",maxrow,i,ABS(a[maxrow][i])); fflush(stdout);
                tmpa = ABS(a[j][i]);
                tmpb = ABS(a[maxrow][i]);
//                printf("tmpa : %f\ntmpb : %f\n",tmpa,tmpb); fflush(stdout);
                if (ABS(a[j][i]) > ABS(a[maxrow][i])){
//                    printf("Seting maxrow to %d\n",maxrow);
                    maxrow = j;
                }
//                printf("GOT HERE blabla\n"); fflush(stdout);
            }
            printf("Found maxrow : %d\n",maxrow); fflush(stdout);
            /* Swap the maxrow and ith row */
            for (k=i;k<size+1;k++) {
                tmpa = a[i][k];
                a[i][k] = a[maxrow][k];
                a[maxrow][k] = tmpa;
            }
            /* Singular matrix? */
            if (ABS(a[i][i]) < EPS)
                return(FALSE);
        }
        
        upc_barrier(b++);
        
        //printf("maxrow : %d\n",maxrow);



        /* Eliminate the ith element of the jth row */
        for (j=i+1;j<size;j++) {
            // change this to upc_forall
            if(MYTHREAD == j % THREADS){
                for (k=size;k>=i;k--) {
                    a[j][k] -= a[i][k] * a[j][i] / a[i][i];
                }
            }
        }
        upc_barrier(b++);
        //WriteSolution(a,n,x);
    }


    /* Do the back substitution */
    if(MYTHREAD == 0){
        for (j=size-1;j>=0;j--) {
            tmpa = 0;
            for (k=j+1;k<size;k++)
                tmpa += a[j][k] * x[k];
            x[j] = (a[j][size] - tmpa) / a[j][j];
        }
    }

    return(TRUE);
}


int main(int argc,char **argv)
{
    int i,j;
    int tmp;
    double dtmp;
    struct timeval tim;
    double t1,t2;
    FILE* fin;

    if(MYTHREAD == 0){
        printf("Opening file %s\n",argv[1]); fflush(stdout);
        if((fin = fopen(argv[1],"r")) == NULL){
            printf("ERROR OPENING FILE\n"); fflush(stdout);
        }
        printf("Opened\n"); fflush(stdout);
        printf("Reading size\n"); fflush(stdout);
        fscanf(fin,"%d",&tmp);
        size = tmp;
        
        printf("Size = %d\n",size); fflush(stdout);

        // get A of (Ax = b)
        printf("Reading...\n"); fflush(stdout);
        for(i = 0;i<size;i++){
            for(j = 0;j<size;j++){
                fscanf(fin,"%lf",&dtmp);
                a[i][j] = dtmp;
            }
        }
        //get b of (Ax = b)
        for(i = 0;i<size;i++){
            fscanf(fin,"%lf",&dtmp);
            a[i][size] = dtmp;
        }
        fclose(fin);
    }
    
    upc_barrier(1);
    printf("Solving...\n"); fflush(stdout);
    gettimeofday(&tim, NULL);
    t1=tim.tv_sec+(tim.tv_usec/1000000.0);
    GSolve();
    gettimeofday(&tim,NULL);
    t2=tim.tv_sec+(tim.tv_usec/1000000.0);
    printf("%.6lf seconds to solve\n",t2-t1);
    upc_barrier(2);
    if(MYTHREAD == 0){
        WriteSolution2();
    }
    
}