#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.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);


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 next itteraton of solution vector, space to be malloced
double  *ttemp; // used to swap t1 and t at each itteration
int     itt_max=5;      // number of itterations to preform
int     itt;            // current itteration
int     i, j;           // indices into arrays
double  sum;    // computes the inner products for A * t
double  error;  // max | t1[i] - t[i] |
double  errori; // | t1[i] - t[i] |
char ch;                // for error checking on command line args.

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


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, &ch) != 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
if(logfile != null)
        log = fopen(logfile,"w+");
else
        log = fopen(STDOUT,"w+");

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

fprintf(log,"\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;
        fprintf(log,"%10.6f ", *(a+n*i+j) );
        }
        fprintf(log,"\n");
}


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

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

// Creat our threads
pthread_t thread[tcount];

// Do stuff!
fprintf(log,"\n  itt  error\n");
for(itt=0; itt<=itt_max; itt++)
{
        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;
                }
        }
        ttemp = t1;
        t1 = t;
        t = ttemp;
        fprintf(log,"%5d %14.6e\n", itt, error);
}

return(0);
}