/* simeq_accuracy.c  try big range of values */

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

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

static double B[992];
static double X[31] = {0.0, 1.0, 2.0, 5.0, 10.0, 50.0, 100.0, 200.0, 500.0, 
                       1000.0, 2000.0, 5000.0, 10000.0, 20000.0, 50000.0,
                       100000.0, 200000.0, 500000.0, 1000000.0, 2000000.0,
                       5000000.0, 10000000.0, 20000000.0, 50000000.0,
                       1.0E8, 2.0E8, 5.0E8, 10.0E8, 20.0E8, 50.0E8, 100.0E8};
static double Xc[31];

void test(int n)
{
  int i, j, sz;
  double bmin, bmax, bavg;

  printf("n=%d, xbig=%e \n", n, X[n-1]);
  sz = n+1;
  for(i=0; i<n; i++) /* initialize random matrix */
    for(j=0; j<n; j++)
      B[i*sz+j] = pow(2.0,n+1)*(udrnrt()-0.5);

  bmax = B[0];
  bmin = B[0];
  bavg = 0.0;
  for(i=0; i<n; i++)
  {
    B[i*sz+n] = 0.0;
    for(j=0; j<n; j++)
    {
      B[i*sz+n] += B[i*sz+j]*X[j]; /* compute RHS, Y */
      bmax = max(bmax,B[i*sz+j]);
      bmin = min(bmin,B[i*sz+j]);
      bavg += B[i*sz+j];
    }
  }
  printf("Amin=%f, Amax=%f, Aavg=%f \n", bmin, bmax, bavg/(double)(n*n));

  simeqb(n,B,Xc);
  for(i=0; i<n; i++)
  {
    printf("%d  comp=%f  true=%f err=%e \n", i, Xc[i], X[i], abs(Xc[i]-X[i]));
  }
  printf(" \n");
} /* end one test */

int main(int argc, char * argv[])
{
  int n;

  printf("simeq_accuracy.c running, big range of values \n");
  printf("A*X=Y, given X, generate random A, compute Y \n");
  printf("solve for computed Xc, compare accuracy to X \n");
  n = 6;
  test(n);
  n = 10;
  test(n);
  n = 19;
  test(n);
  n = 25;
  test(n);
  n = 31;
  test(n);
  printf(" \n");
  printf("note that some of the error comes from multiplying A*X \n");
  printf("end simeq_accuracy.c \n");
  return 0;
} /* end simeq_accuracy.c */