/* test_simeq_back.c  compare reduction to back substitution */

#include <stdio.h>
#include <time.h>
#include "udrnrt.h"
#include "simeq_back.h"
#include "simeq_lup.h"
#include "simeq.h"

#define N 2049
static double A[N*N];
static double X[N];
static double XB[N];
static double XP[N];
static double Y[N];

#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();
  }
}

int main(int argc, char * argv[])
{
  int i, j, k;
  int n;
  double t, total;
  double start, now, next;
  
  printf("test_simeq_back.c running\n");
  start = (double)time(NULL);
  for(n=256; n<N; n=n*2)
  {
    init_matrix(n, A, Y);
    printf("matrix initialized at %g\n", now-start);

    now = (double)time(NULL);
    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]);
    }
    printf("n=%d, total error=%g, simeq      took %g seconds\n",
           n, total, next-now);

    now = (double)time(NULL);
    simeq_back(n, A, Y, XB);
    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]*XB[j];
      total = total + abs(t-Y[i]);
    }
    printf("n=%d, total error=%g, simeq back took %g seconds\n",
           n, total, next-now);

    now = (double)time(NULL);
    simeq_lup(n, A, Y, XP);
    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]*XP[j];
      total = total + abs(t-Y[i]);
    }
    printf("n=%d, total error=%g, simeq lup  took %g seconds\n",
           n, total, next-now);
    printf("\n");
  }
  return 0;
} /* end main time_simeq_back.c */