// matmul_heap.c
// CS4540 Fall 2010
// kapenga
//
// This program is one of a set of three programs that show a
// matrix itteration. The difference between the three is the way
// space for the arrays is allocated:
//   matmul_stack uses local (automatic) variables for the arrays,
//        which use space on the stack for the arrays. The array
//        sizes are dynamic.
//   matmul_static uses global (external) variables for the arrays,
//        which uses space on the stack for the arrays. The maximum
//        array sizes are compiled in.
//   matmul_heap uses dynamic (malloced) space for the arrays,
//        which uses space on the heap for the arrays. The array
//        sizes are dynamic.
//
// The following itteretion can be used to solve linear systems
//   t_{i+1} = A t_i + b
// If the itteration converges to t, then t == t_{i+1} == t_i
// So t = A t + b
//   or  (I-a) t = b
//   where, I is the n*n idenity matrix
// There are several important applied problems where convergence
// will take place. One such case is when for
// each row of A ( rows 0 <= i < n)
//             sum(j=0 ... n-1) abs(a[i][j])  < 1.0
// Then the itteration will converge, assuming no roundoff or overflow.
// Example
// % ./matmul_heap 4 10 5
//
//  a=
//  0.189331   0.147829  -0.009582   0.012830
// -0.020409   0.222627   0.073037   0.042701
//  0.069882   0.228326  -0.001161   0.024936
//  0.116375  -0.100117   0.229832   0.022235
//
//  b=
//  2.411774   9.837874   6.251698   6.576916
//
//  itt  error
//    0   2.878398e+00
//    1   8.266521e-01
//    2   2.688652e-01
//    3   8.817662e-02
//    4   2.832084e-02
//    5   9.015857e-03
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>
#include <fcntl.h>
#include <unistd.h>

// These two function are not ansi C so they do not appear from the
// libstd.h  header if the gcc option -std=c99 option is used.
// I should see if there is a safe way to include them from stdlib.h
// and not place them explicitly here, which is bad style.

void srand48(long int seedval);
double drand48(void);

struct theDatas{
    int tid;
    int tcount;
    int itt;
    int n;
    double *a;
    double *b;
    double *t;
    double *t1;
    double sum;
    double error;
};
void *doTheStuff(void *datas)
{
    double  sum;    // computes the inner products for A * t
    double  error;  // max | t1[i] - t[i] |
    double  errori; // | t1[i] - t[i] |
    int i, j;       // indices into arrays

    int tid;
    int itt;
    int n;
    int tcount;
    double *a;
    double *b;
    double *t;
    double *t1;

    struct theDatas *localDatas;

    localDatas = (struct theDatas *) datas;
    tid = localDatas->tid;
    itt = localDatas->itt;
    tcount = localDatas->tcount;
    n = localDatas->n;
    a = localDatas->a;
    b = localDatas->b;
    t = localDatas->t;



    printf("I'm thread %d!\n",tid);

    error=0.0;
    for(i=0; i< n; i += tcount)
    {
        printf("Thread %d Processing Row %d!\n",tid,i);
        sum = 0.0;
        for(j=0; j< n; j++)
        {
            printf("Thread %d Processing Row %d Col %d!\n",tid,i,j);
            sum += *(a+n*i+j) * t[j];
        }
        t1[i] = sum + b[i];
        errori = fabs(t1[i]-t[i]);
        if(errori > error)
        {
            error=errori;
        }
    }
    localDatas.sum = sum;
    localDatas.error = error;

    pthread_exit(NULL);
}

int main(int argc, char *argv[])
{
int n=4;        // problenm size
double  *a;     // pointer to transfromation array, space to be malloced
double  *b;     // pointer to transfromation vector, space to be malloced
int seed=10;    // seed for srand48() / drand48()
double *t;     // pointer to solution vector, space to be malloced
double *t1;    // pointer to solution vector, space to be malloced
int itt_max=5;  // number of itterations to preform
int itt;        // current itteration
char ch;        // for error checking on command line args.
int i, j;       // indices into arrays
double  *ttemp; // used to swap t1 and t at each itteration
double  sum;    // computes the inner products for A * t
double  error;  // max | t1[i] - t[i] |
double  errori; // | t1[i] - t[i] |

// New Variables
int tcount = 2; // Number of threads, default 2
char* logfile=NULL; // File to log too, default will be STDOUT
FILE * log;     // File descriptor for our log file
int tc;         // thread counter
int rc;

if( argc == 4 )
{
    if( (sscanf(argv[1],"%d %[^ /t]", &n, &ch) != 1) || (sscanf(argv[2],"%d %[^ /t]", &seed, &ch) != 1) || (sscanf(argv[3],"%d %[^ /t]", &itt_max, &ch) != 1) )
    {
        fprintf(stderr," ERROR : useage: %s [ <n> <seed> <itt_max>]\n", argv[0]);
        return(1);
    }
}
else if(argc == 5)
{
    if(sscanf(argv[4],"%d %[^ /t]", &tcount, &ch) != 1)
    {
        fprintf(stderr," ERROR : useage: %s [ <n> <seed> <itt_max>]\n", argv[0]);
        return(1);
    }
}
else if(argc == 6)
{
    if(sscanf(argv[5],"%s", logfile) != 1)
    {
        fprintf(stderr," ERROR : useage: %s [ <n> <seed> <itt_max>]\n", argv[0]);
        return(1);
    }
}
else if(argc != 1 )
{
    fprintf(stderr," ERROR : useage: %s [ <n> <seed> <itt_max>]\n", argv[0]);
    return(1);
}

if( n<1 )
{
    fprintf(stderr," ERROR : n must be positive\n");
    return(1);
}

if( (a=(double *)malloc(sizeof(double)*n*n)) == NULL )
{
    fprintf(stderr," ERROR : malloc for a failed\n");
    return(1);
}

if( (b=(double *)malloc(sizeof(double)*n)) == NULL )
{
    fprintf(stderr," ERROR : malloc for b failed\n");
    return(1);
}

if( (t=(double *)malloc(sizeof(double)*n)) == NULL )
{
    fprintf(stderr," ERROR : malloc for t failed\n");
    return(1);
}
if( (t1=(double *)malloc(sizeof(double)*n)) == NULL )
{
    fprintf(stderr," ERROR : malloc for t1 failed\n");
    return(1);
}

// Open Log File
/*printf("%s\n",logfile);
if(logfile != NULL)
{
    log = fopen(logfile,O_APPEND | O_WRONLY);
    dup2(log,stdout);
}*/

// Generate matrix a with | eigenvalues | < 1
srand48((long int)seed);

printf("\n  a=\n");
for(i=0; i< n; i++)
{
    for(j=0; j< n; j++)
    {
        *(a+n*i+j) = 1.999 * (drand48() - 0.5) / n;
        printf("%10.6f ", *(a+n*i+j) );
    }
    printf("\n");
}


// Generate vector b
printf("\n  b=\n");
for(i=0; i< n; i++)
{
    b[i] = 10.0 * drand48();
    printf("%10.6f ", b[i]);
}
printf("\n");

// Initialize t
for(i=0; i< n; i++)
{
    t[i] = b[i];
}

// Create our thread array
pthread_t thread[tcount];
struct theDatas someData[tcount];
pthread_attr_t attr;
void *status;

pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

// Do stuff!
printf("\n  itt  error\n");

for(itt=0; itt<=itt_max; itt++)
{
    for(tc=0; tc<tcount; tc++)
    {
        someData[tc].tid = tc;
        someData[tc].tcount = tcount;
        someData[tc].itt = itt;
        someData[tc].n = n;
        someData[tc].a = a;
        someData[tc].b = b;
        someData[tc].t = t;
        someData[tc].t1 = t1;

        printf("Spawning thread %d of %d\n",tc,tcount-1);

        rc = pthread_create(&thread[tc], NULL, doTheStuff, (void *) &someData[tc]);
        if(rc)
        {
            fprintf(stderr,"Error spawning thread %d (%d)\n",tc,rc);
            return(-1);
        }

        /*
        error=0.0;
        for(i=0; i< n; i++)
        {
            sum = 0.0;
            for(j=0; j< n; j++)
            {
                sum += *(a+n*i+j) * t[j];
            }
            t1[i] = sum + b[i];
            errori = fabs(t1[i]-t[i]);
            if(errori > error)
            {
                error=errori;
            }
        }
        */

    }
    for(tc=0; tc<tcount; tc++)
    {
        rc = pthread_join(thread[tc], &status);
        if(rc)
        {
            fprintf(stderr,"Error joining thread %d (%d)\n",tc,rc);
            return(-1);
        }
         printf("Main: completed join with thread %d having a status of %ld\n",tc,(long)status);
    }

    ttemp = t1;
    t1 = t;
    t = ttemp;
    printf("%5d %14.6e\n", itt, error);
}

return(0);
}

// output
travis@batman:~/cs454/ass2/3[139] 2864 % ./matmul_heap 10 10 10

  a=
  0.075732   0.059132  -0.003833   0.005132  -0.008163   0.089051   0.029215   0.017080   0.027953   0.091330
 -0.000464   0.009974   0.046550  -0.040047   0.091933   0.008894  -0.051739   0.096709   0.025021   0.031523
 -0.076330  -0.076081  -0.068221  -0.056610  -0.058096   0.054306  -0.079693  -0.005993  -0.010118   0.029982
  0.026604  -0.027442   0.003887   0.094350  -0.026043  -0.024950   0.087433  -0.075697   0.005936  -0.085100
  0.000359  -0.076407  -0.045923   0.098488   0.055573  -0.081033   0.056287  -0.061312   0.070243   0.071067
  0.087992   0.060836   0.084124   0.019586  -0.049597  -0.007719   0.064988  -0.049935  -0.059120   0.079362
 -0.053771   0.057492  -0.009375  -0.010235   0.046699   0.073086   0.077445  -0.005320  -0.024860  -0.075497
 -0.068074  -0.067880  -0.047793  -0.032344   0.058348  -0.093961   0.054170  -0.085748  -0.051241  -0.079509
  0.095479  -0.007781  -0.090568   0.075227   0.067584  -0.062678  -0.052532   0.001142   0.030670  -0.025270
 -0.067617   0.094996   0.026292   0.094573  -0.012125  -0.029150   0.023878  -0.037532  -0.022051   0.083907

  b=
  6.298272   3.537073   8.809945   8.331442   4.834493   5.205376   9.472396   2.323569   1.059944   6.950048

  itt  error
Spawning thread 0 of 1
I'm thread 0!
Thread 0 Processing Row 0!
Thread 0 Processing Row 0 Col 0!
Thread 0 Processing Row 0 Col 1!
Thread 0 Processing Row 0 Col 2!
Thread 0 Processing Row 0 Col 3!
Thread 0 Processing Row 0 Col 4!
Thread 0 Processing Row 0 Col 5!
Thread 0 Processing Row 0 Col 6!
Thread 0 Processing Row 0 Col 7!
Thread 0 Processing Row 0 Col 8!
Thread 0 Processing Row 0 Col 9!
[2]    19581 segmentation fault  ./matmul_heap 10 10 10