#include <upc_relaxed.h>
#include <stdio.h>
#include <upc.h>
#include<upc_strict.h>

#define N 4
#define P 4
#define M 4
#define D N/THREADS

shared double a[1*THREADS][N];
shared double x[1*THREADS];
shared double buf[1*THREADS];
shared double b[1*THREADS];
const double tau = .5;
const double epsilon = 1e-6;
double normBsquared;

int main(void) {

  int i,j;
  double normXSquared = 0.;


// Arrays initialization   by thread 0

    if(MYTHREAD==0) {
            normBsquared = 0;
            printf("Matrix A\n");
                for(i=0;i<N;i++)
            {
                        for(j=0;j<P;j++)
                {
                    a[i][j] = i * N +j ;
                    printf("%f\t",a[i][j]);
                    if(j==(N-1))
                        b[j]= 3*(j+1)-1;
                    else
                        b[j]= 4*(j+1);
                    normBsquared += b[j] * b[j];
                }
                printf("\n");
            }
            printf("VECTOR B\n");
            for(i<0;i<N;i++)
            {
                printf("%f\t",b[i]);
            }
            printf("\n");
            for(i=0;i<N;i++)
            {
                x[i]=0;
                buf[i]=0;
            }
            printf("Start find!\n");
            printf("Bsquared %f\n", normBsquared);
    }
    for(;;)
    {
        upc_barrier;
        upc_forall(i = 0; i < N; i++;i)
            for (int j = 0; j < N; j++)
                buf[i] += a[i][j] * x[j];
        upc_forall(i = 0; i < N; i++;i)
        {
            buf[i] -= b[i];
            x[i] -= buf[i] * tau;
        }
        upc_barrier;
        if(MYTHREAD==0)
        {
            for(i = 0; i < N; i++)
            {
                normXSquared += buf[i] * buf[i];
            }
            //printf("current Xsquared %f\n", normXSquared);
        }
        upc_barrier;
        if (normXSquared / normBsquared < epsilon * epsilon)
            break;
        if(MYTHREAD==0)
        {
            for(i=0;i<N;i++)
            {
                buf[i]=0;
                normXSquared = 0;
            }
        }
    }
    if(MYTHREAD==0)
        for(i=0;i<N;i++)
        {
            printf("%f\t",x[i]);
        }

    /*upc_barrier;

    // all threads  perform matrix multiplication

    upc_forall(i=0;i<N;i++;i)
    {
      // &a[i][0] specifies that this iteration will be executed by the thread that has affinity to
    //  element a[i][0]
        //printf("Matrix element %d THREAD (%d)",i,MYTHREAD);
         for (j=0; j<M; j++) {
            c[i][j] = 0;
            for(l=0; l< P; l++) c[i][j] +=a[i][l]*b[l][j];
         }
    }*/

    //upc_barrier;

}

/*void solve1()
{

    int i,j;
    //initialize x
    if(MYTHREAD==0)
        for(i=0;i<N;i++)
        {
            x[i]=0;
            buf[i]=0;
        }
    for(;;)
    {
        upc_barrier;
        upc_forall(i = 0; i < N; i++;i)
            for (int j = 0; j < N; j++)
                buf[i] += a[i][j] * x[j];
        upc_forall(i = 0; i < N; i++;i)
        {
            buf[i] -= b[i];
            x[i] -= buf[i] * tau;
        }
        upc_barrier;
        if(MYTHREAD==0)
        {
            for(i = 0; i < N; i++)
            {
                normXSquared += buf[i] * buf[i];
            }
        }
        if (normXSquared / normBsquared < epsilon * epsilon)
            break;
        if(MYTHREAD==0)
        {
            for(i=0;i<N;i++)
            {
                buf[i]=0;
                normXSquared = 0;
            }
        }
    }
}*/