/**************************************************************************************
                            C-DAC Tech Workshop : hyPACK-2013
                             October 15-18,2013
 Example    : pthread-jacobi.c
 
 Objective      : Jacobi method to solve AX = b matrix system of linear equations.
 
 Input          : Class Size
          Number of Threads
 
 Output         : The solution of  Ax = b or
                  The status of convergence for given bumber of iterations
 
 Created    : MAY-2013  
 
 E-mail         : hpcfte@cdac.in                                         
  
*************************************************************************************/
 
#include <stdio.h>
#include<pthread.h>
#include<sys/time.h>
#include<stdlib.h>
#include<string.h>
 
#define  MAX_ITERATIONS 100
#define MAXTHREADS 8
 
double   Distance(double *X_Old, double *X_New, int matrix_size);
pthread_mutex_t  mutex1 = PTHREAD_MUTEX_INITIALIZER;
 
double   **Matrix_A, *RHS_Vector;
double   *X_New, *X_Old, *Bloc_X, rno,sum;
 
int Number;
void jacobi(int);
 
int main(int argc, char **argv)
{
 
        double diag_dominant_factor  = 4.0000;
        double tolerance  = 1.0E-5;
        /* .......Variables Initialisation ...... */
        int matrix_size,  NoofRows, NoofCols,CLASS_SIZE,THREADS;
        int NumThreads,ithread;
        double rowsum;
        double  sum;
        int irow, icol, index, Iteration,iteration,ret_count;
        double time_start, time_end,memoryused;
        struct timeval tv;
        struct timezone tz;
        char * CLASS;
        FILE *fp;
 
        pthread_attr_t pta;
        pthread_t *threads;
      
        memoryused =0.0;
 
 
        printf("\n\t\t---------------------------------------------------------------------------");
        printf("\n\t\t Centre for Development of Advanced Computing (C-DAC)");
        printf("\n\t\t Email : hpcfte@cdac.in");
        printf("\n\t\t---------------------------------------------------------------------------");
        printf("\n\t\t Objective : To Solve AX=B Linear Equation (Jacobi Method)\n ");
        printf("\n\t\t Performance for solving AX=B Linear Equation using JACOBI METHOD");
        printf("\n\t\t..........................................................................\n");
 
        if( argc != 3 ){
 
            printf("\t\t Very Few Arguments\n ");
            printf("\t\t Syntax : exec <Class-Size (Give A/B/C)> <Threads>\n");
            exit(-1);
        }
        else {
            CLASS = "A";
            THREADS = 2;
        }
    if (THREADS > MAXTHREADS ){
        printf("\n Number of Threads must be less than or equal to 8. Aborting ...\n");
        return 0;
    }
        if( strcmp(CLASS, "A" )==0){
            CLASS_SIZE = 1024;
        }
        if( strcmp(CLASS, "B" )==0){
            CLASS_SIZE = 2048;
        }
        if( strcmp(CLASS, "C" )==0){
            CLASS_SIZE = 4096;
        }
 
 
         matrix_size = CLASS_SIZE;
         NumThreads = THREADS;
         printf("\n\t\t Matrix Size :  %d",matrix_size);
         printf("\n\t\t Threads     :  %d",NumThreads);
           
        NoofRows = matrix_size;
        NoofCols = matrix_size;
         
        
        /* Allocate The Memory For Matrix_A and RHS_Vector */
        Matrix_A = (double **) malloc(matrix_size * sizeof(double *));
        RHS_Vector = (double *) malloc(matrix_size * sizeof(double));
 
 
        /* Populating the Matrix_A and RHS_Vector */
        rowsum = (double) (matrix_size *(matrix_size+1)/2.0);
        for (irow = 0; irow < matrix_size; irow++)
        {
             Matrix_A[irow] = (double *) malloc(matrix_size * sizeof(double));
             sum = 0.0;
             for (icol = 0; icol < matrix_size; icol++)
             {
                 Matrix_A[irow][icol]= (double) (icol+1);
                 if(irow == icol )  Matrix_A[irow][icol] = (double)(rowsum);
                 sum = sum + Matrix_A[irow][icol];
             }
             RHS_Vector[irow] = (double)(2*rowsum) - (double)(irow+1);
         }
           
         memoryused+=(NoofRows * NoofCols * sizeof(double));
         memoryused+=(NoofRows * sizeof(double)); 
          
         printf("\n");
 
         if (NoofRows != NoofCols)
         {
                printf("Input Matrix Should Be Square Matrix ..... \n");
                exit(-1);
         }
 
        /* Dynamic Memory Allocation */
        X_New = (double *) malloc(matrix_size * sizeof(double));
        memoryused+=(NoofRows * sizeof(double));
        X_Old = (double *) malloc(matrix_size * sizeof(double));
        memoryused+=(NoofRows * sizeof(double));
        Bloc_X = (double *) malloc(matrix_size * sizeof(double));
        memoryused+=(NoofRows * sizeof(double));
 
        /* Calculating the time of Operation Start */
        gettimeofday(&tv, &tz);
        time_start= (double)tv.tv_sec + (double)tv.tv_usec / 1000000.0;
 
        /* Initailize X[i] = B[i] */
 
        for (irow = 0; irow < matrix_size; irow++)
        {
                Bloc_X[irow] = RHS_Vector[irow];
                X_New[irow] =  RHS_Vector[irow];
        }
   
        /* Allocating the memory for user specified number of threads */
        threads = (pthread_t *) malloc(sizeof(pthread_t) * NumThreads); 
 
        /* Initializating the thread attribute */
        pthread_attr_init(&pta);
 
        Iteration = 0;
        do {
              for(index = 0; index < matrix_size; index++)
              X_Old[index] = X_New[index];
              for(ithread=0;ithread<NumThreads;ithread++)
              {
 
                 /* Creating The Threads */                
                 ret_count=pthread_create(&threads[ithread],&pta,(void *(*) (void *))jacobi, (void *) (matrix_size/NumThreads));
             if(ret_count)
             {
            printf("\n ERROR : Return code from pthread_create() is %d ",ret_count);
            exit(-1);
             }
 
               }
 
                  Iteration++;
                  for (ithread=0; ithread<NumThreads; ithread++)
                  {
                    ret_count=pthread_join(threads[ithread], NULL);
            if(ret_count)
                {
                printf("\n ERROR : Return code from pthread_join() is %d ",ret_count);
                exit(-1);
                }
                  }                 
                  ret_count=pthread_attr_destroy(&pta);
              if(ret_count)
              {
            printf("\n ERROR : Return code from pthread_attr_destroy() is %d ",ret_count);
            exit(-1);
              }
             } while ((Iteration < MAX_ITERATIONS) && (Distance(X_Old, X_New, matrix_size) >= tolerance));
        
        /* Calculating the time at the end of Operation */
        gettimeofday(&tv, &tz);
        time_end= (double)tv.tv_sec + (double)tv.tv_usec / 1000000.0;
 
        printf("\n\t\t The Jacobi Method For AX=B .........DONE");
        printf("\n\t\t Total Number Of Iterations   :  %d",Iteration);
        printf("\n\t\t Memory Utilized              :  %lf MB",(memoryused/(1024*1024)));
        printf("\n\t\t Time in  Seconds (T)         :  %lf",(time_end - time_start));
        printf("\n\t\t..........................................................................\n");
 
        /* Freeing Allocated Memory */
        free(X_New);
        free(X_Old);
        free(Matrix_A);
        free(RHS_Vector);
        free(Bloc_X);
        free(threads);
 
}
 
double Distance(double *X_Old, double *X_New, int matrix_size)
{
        int             index;
        double          Sum;
 
        Sum = 0.0;
        for (index = 0; index < matrix_size; index++)
        Sum += (X_New[index] - X_Old[index]) * (X_New[index] - X_Old[index]);
        return (Sum);
}
 
void jacobi(int Number)
{
   int i,j;
 
   for(i = 0; i < Number; i++)
   {
        Bloc_X[i] = RHS_Vector[i];
 
        for (j = 0;j<i;j++)
        {
            Bloc_X[i] -= X_Old[j] * Matrix_A[i][j];
        }
        for (j = i+1;j<Number;j++)
        {
            Bloc_X[i] -= X_Old[j] * Matrix_A[i][j];
        }
        Bloc_X[i] = Bloc_X[i] / Matrix_A[i][i];
  }
  for(i = 0; i < Number; i++)
 
  {
    X_New[i] = Bloc_X[i];
  }
}