/* test_simeq_small.c  small matrix, little input error */
/*                     big difference in output, not seen in error */
#include <stdio.h>
#include "simeq.h"

static double A[2][2];
static double Y[2];
static double X[2];

int main(int argc, char *argv[])
{
  double err;
  int i, j;
  int n = 2;

  printf("test_simeq_small.c results should be exactly X[0]=1.0, X[1]=-1.0 \n");
  printf("err is the error multiplying the solution, X, times matrix A*X=Y \n");
  printf("a small change in data gives a big change in solution. \n");
  printf("the big change is not indicated by the computed error. \n\n");

  printf("decimal .835*X[0] + .667*X[1] = .168 \n");
  printf("decimal .333*X[0] + .266*X[1] = .067 \n");
  A[0][0] = 0.835; /* row 0 */
  A[0][1] = 0.667;
  Y[0] = 0.168;
  A[1][0] = 0.333; /* row 1 */
  A[1][1] = 0.266;
  Y[1] = 0.067;
  printf("flt-pt %f*X[0] + %f*X[1] = %f \n", A[0][0], A[0][1], Y[0]);
  printf("flt-pt %f*X[0] + %f*X[1] = %f \n", A[1][0], A[1][1], Y[1]);
  printf("initialization complete, solving \n");
  simeq2( n, A , Y , X );
  for(i=0; i<n; i++)
  {
    err = -Y[i];
    for(j=0; j<n; j++)
      err = err + A[i][j]*X[j]; 
    printf("solution X[%1d]=%g, err=%g \n", i, X[i], err);
  }
  printf(" \n");
  printf("now perturb .067 to .066, .001 change, and run again \n");
  Y[1] = 0.066;
  printf("flt-pt %f*X[0] + %f*X[1] = %f \n", A[0][0], A[0][1], Y[0]);
  printf("flt-pt %f*X[0] + %f*X[1] = %f \n", A[1][0], A[1][1], Y[1]);
  simeq2( n, A , Y , X );
  for(i=0; i<n; i++)
  {
    err = -Y[i];
    for(j=0; j<n; j++)
      err = err + A[i][j]*X[j]; 
    printf("solution X[%1d]=%g, err=%g \n", i, X[i], err);
  }
  printf("\n notice computed errors are still small, \n");
  printf(" yet, the values of X[0] and X[1] are very different.\n");
  printf(" Typical problem with people and computers: \n");
  printf(" Often called garbage in, garbage out! \n");
  printf("\n test_simeq_small.c finished");

  return 0;
} /* end test_simeq_small.c */