/* time_simeq.c */

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

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

#define N 4097
static double A[N*N];
static double X[N];
static double Y[N];

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<5)
  {
    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 start, now, next;
  
  printf("time_simeq.c running\n");
  start = (double)time(NULL);
  for(n=4; n<N; n=n*2)
  {
    init_matrix(n, A, Y);
    now = (double)time(NULL);
    printf("matrix initialized at %g\n", now-start);
    simeq(n, A, Y, X);
    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, simeq took %g seconds\n",
           n, total, next-now);
  }
  return 0;
} /* end main time_simeq.c */