/* test_simeq2.c  tough matrix */
#include <stdio.h>
#include "simeq.h"

static void vec_put(int n, double A[]);
static void mat_put(int n, double A[]);
static double AA[100000000]; /* 10,000 by 10,000 */
static double YY[10000];
static double XX[10000];

int main(int argc, char *argv[])
{
  double A[3][3] = {{ 3.0, 2.0, 100000.0},
                    {-1.0, 3.0, 100000.0},
                    { 1.0, 2.0,     -1.0}};
  double Y[3] = {100005.0, 100002.0, 2.0 };
  double X[3];
  double err;
  int i, j, k, n;

  printf( "initial matrix A \n" );
  mat_put( 3, (double *)A );
  printf( "initial vector Y \n" );
  vec_put( 3, Y );
  printf( "solve equations \n" );
  simeq( 3, (double *)A , Y , X );
  printf( "solution X \n" );
  vec_put(3, X );

  for(n=2; n<=1024; n=n*2)
  {
    printf("initializing big matrix, n=%d, n*n=%d \n", n, n*n);
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
      {
        AA[i*n+j] = 1.0/((double)(i+1)+(double)(j+1));
      }
      YY[i] = 1.0/(double)(n*n);
    }
    printf("initialization complete, solving \n");
    if(n<=4) mat_put( n, (double *)AA );
    if(n<=4)vec_put( n, YY );

    simeq( n, (double *)AA , YY , XX );
    err = -YY[0];
    for(j=0; j<n; j++) err = err + AA[j]*XX[j]; 
    printf("solution XX[0]=%g, err=%g \n", XX[0], err);
    err = -YY[n-1];
    for(j=0; j<n; j++) err = err + AA[(n-1)*n+j]*XX[j]; 
    printf("       XX[n-1]=%g, err=%g \n\n", XX[n-1], err);
  }
}


void vec_put(int n, double A[])
{
  int i;
  for(i=0; i<n; i++) printf("VEC[%d] = %lf \n", i, A[i]);
}


void mat_put(int n, double A[])
{
  int i, j;
  for(i=0; i<n; i++) 
    for(j=0; j<n; j++)
      printf("MAT[%d,%d] = %lf \n", i, j, A[i*n+j]);
}