// 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 *error;
    pthread_mutex_t *errorLock;
};
void *doTheStuff(void *datas)
{
    double  sum;    // computes the inner products for A * t
    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;
    double *error;
    pthread_mutex_t *errorLock;

    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;
    t1 = localDatas->t1;
    error = localDatas->error;
    errorLock = localDatas->errorLock;


    //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;
        //printf("n = %d\n",n);
        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];
        }
        //printf("0 ");
        t1[i] = sum + b[i];
        //printf("1 ");
        errori = fabs(t1[i]-t[i]);
        //printf("2 ");
        pthread_mutex_lock(&errorLock);
        if(errori > error)
        {
            error=errori;
        }
        pthread_mutex_unlock(&errorLock);
        //printf("3 ");
    }

    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  error;  // max | t1[i] - t[i] |

// New Variables
int tcount = 4; // 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[1],"%d %[^ /t]", &n, &ch) != 1) || (sscanf(argv[2],"%d %[^ /t]", &seed, &ch) != 1) || (sscanf(argv[3],"%d %[^ /t]", &itt_max, &ch) != 1) || (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[1],"%d %[^ /t]", &n, &ch) != 1) || (sscanf(argv[2],"%d %[^ /t]", &seed, &ch) != 1) || (sscanf(argv[3],"%d %[^ /t]", &itt_max, &ch) != 1) || (sscanf(argv[4],"%d %[^ /t]", &tcount, &ch) != 1) || (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];
}

// Thread stuff!
pthread_t thread[tcount];
struct theDatas someData[tcount];
pthread_attr_t attr;
void *status;
pthread_mutex_t errorLock = PTHREAD_MUTEX_INITIALIZER;


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;
        someData[tc].error = error;
        someData[tc].errorLock = errorLock;

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

        pthread_create(&thread[tc], NULL, doTheStuff, (void *) &someData[tc]);
    }
    for(tc=0; tc<tcount; tc++)
        rc = pthread_join(thread[tc], &status);

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

return(0);
}

// Output
gcc -Wall -pedantic -std=c99 -o matmul_heap matmul_heap.c -lm -lpthread -fstack-protector-all -g
matmul_heap.c: In function ‘doTheStuff’:
matmul_heap.c:127: warning: passing argument 1 of ‘pthread_mutex_lock’ from incompatible pointer type
matmul_heap.c:128: error: invalid operands to binary > (have ‘double’ and ‘double *’)
matmul_heap.c:130: error: incompatible types in assignment
matmul_heap.c:132: warning: passing argument 1 of ‘pthread_mutex_unlock’ from incompatible pointer type
matmul_heap.c: In function ‘main’:
matmul_heap.c:284: error: incompatible types in assignment
matmul_heap.c:285: error: incompatible types in assignment
matmul_heap.c:157: warning: unused variable ‘log’
make: *** [matmul_heap] Error 1