// time_crout.c 
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "crout.h"
#include "udrnrt.h"

#define N  4096
static double A[N][N]; // destroyed
static double AA[N][N]; // for checking
static double X[N];
static double X_soln[N];
static double Y[N];
static double YY[N];

#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

void init_matrix(int n)
{
  double sum;
  int i, j;

  for(i=0; i<n; i++)
  {
    X_soln[i] = udrnrt();
    for(j=0; j<n; j++)
    {
      A[i][j] = udrnrt();
      AA[i][j] = A[i][j];
    } // j
  } // i
  for(i=0; i<n; i++)
  {
    sum = 0.0;
    for(j=0; j<n; j++)
      sum = sum + A[i][j]*X_soln[j];
    Y[i] = sum;
    YY[i] = Y[i];
  } // i
    
  if(n<5)
  {
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
        printf("A[%d][%d]=%f\n", i, j, A[i][j]);
      printf("X_soln[%d]=%f, Y[%d]=%f\n", i, X_soln[i], i, Y[i]);
      printf("\n ");
    } // i
  }
} // end init_matrix

int main(int argc, char *argv[])
{
  double sum, sumerr, maxerr;
  double start, now, next;
  int i, j, n;
  
  printf("time_crout.c running\n");
  start  = (double)time(NULL);
  for(n=4; n<N; n=n*2)
  {
    init_matrix(n);
    now  = (double)time(NULL);
    printf("matrix initialized in %f seconds\n", now-start);
    crout(n, N,  A, Y, X);
    next  = (double)time(NULL);
    sumerr = 0.0;
    maxerr = 0.0;
    for(i=0; i<n; i++)
    {
      sum = 0.0;
      for(j=0; j<n; j++) sum = sum + AA[i][j]*X[j];
      sumerr = sumerr + abs(sum-YY[i]);
      maxerr = max(maxerr, abs(X_soln[i]-X[i]));
      if(n<5) printf("X[%d]=%f, err=%e\n", i, X[i], X_soln[i]-X[i]);
    } // i
    printf("n=%d, sumerr=%e, maxerr=%e\n", n, sumerr, maxerr);
    printf("crout took %f seconds\n", next-now);
  } // n
  return 0;
} // end time_crout.c