/* time_simeq_lup.c */

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

#define N 4096
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_lup.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_lup(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[j]);
      if(n<5) printf("X[%d]=%g, t=%g, Y[%d]=%g \n",
                      i, X[i], t, i, Y[i]);
    }
    printf("n=%d, total error=%g, simeq_lup took %g seconds\n",
           n, total, next-now);
  }
  return 0;
} /* end main time_simeq_lup.c */