/* time_tsimeq.c  threaded version */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "udrnrt.h"
#include "tsimeq.h"

#define N 4097
static double A[N*N];
static double X[N];
static double Y[N];
static int np=11;

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))

static void init_matrix(int n, double A[], double Y[])
{
  int i, j;
  for(i=0; i<n; i++)
  {
    Y[i] = udrnrt();
    for(j=0; j<n; j++)
      A[i*n+j] = udrnrt();
  }
  if(n<9)
  {
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++) printf("%5.2f ", A[i*n+j]);
      printf("\n");
      printf("Y[%d]=%5.2f\n", i, Y[i]);
    }
  }
}

int main(int argc, char * argv[])
{
  int i, j, k;
  int n;
  double t, total;
  double wall_start, wall_now, wall_next, cpu_start, cpu_now, cpu_next;
  
  printf("time_tsimeq.c running\n");
  cpu_start = (double)clock()/(double)CLOCKS_PER_SEC;
  wall_start = (double)time(NULL);
  for(n=8; n<=N; n=n*2)
  {
    init_matrix(n, A, Y);
    cpu_now = (double)clock()/(double)CLOCKS_PER_SEC;
    wall_now = (double)time(NULL);
    printf("matrix initialized at wall %g, cpu %g\n",
           wall_now-wall_start, cpu_now-cpu_start);
    tsimeq(np, n, A, Y, X);
    cpu_next = (double)clock()/(double)CLOCKS_PER_SEC;
    wall_next = (double)time(NULL);
    total = 0.0;
    for(i=0; i<n; i++)
    {
      t = 0.0;
      for(j=0; j<n; j++) t = t + A[i*n+j]*X[j];
      total = total + abs(t-Y[i]);
      if(n<9) printf("X[%d]=%g, t=%g, Y[%d]=%g \n",
                      i, X[i], t, i, Y[i]);
    }
    printf("n=%d, total error=%g, tsimeq took wall %g, cpu %g seconds\n",
           n, total, wall_next-wall_now, cpu_next-cpu_now);
  }
  printf("average number of processors used = %f \n",
         (cpu_next-cpu_now)/(wall_next-wall_now));
  printf("time_tsimeq.c exiting \n");
  return 0;
} /* end main time_tsimeq.c */