// test_crout.java  solve linear system of equations
//                  |A|x = y  solve for x  method: Crout

public class test_crout
{
  public test_crout()
  {
    System.out.println("test_crout.java running");
    for(int n=4; n<=128; n=2*n)   // total number of terms
    {
      double A[][] = new double[n][n];
      double AA[][] = new double[n][n]; // A destroyed
      double y[] = new double[n];
      double yy[] = new double[n]; // c destroyed
      double x[] = new double[n];
      double x_soln[] = new double[n];
      double sum, err, sumerr;
      boolean debug = false;

      // build first test case
      if(n==4) System.out.println("random matrix and solution");
      for(int i=0; i<n; i++) // all subscripts zero based
      {
        for(int j=0; j<n; j++)
        {
	  A[i][j] = Math.random();
          AA[i][j] = A[i][j];
	  if(n==4) System.out.println("A["+i+"]["+j+"]="+A[i][j]);
        }
        x_soln[i] = Math.random();
      }

      for(int i=0; i<n; i++)
      {
        sum = 0.0;
        for(int j=0; j<n; j++)
        {
	  sum = sum + A[i][j]*x_soln[j];
        } // j
        y[i] = sum;
        yy[i] = y[i];
        if(n==4) System.out.println("y["+i+"]="+y[i]+", x="+x_soln[i]);
      } // i
   
      new crout(n, A, y, x);

      if(n==4)
      {
        System.out.println("A destroyed");
        for(int i=0; i<n; i++)
	  for(int j=0; j<n; j++)
            System.out.println("A["+i+"]["+j+"]="+A[i][j]);
      }
      // check accuracy
      sumerr = 0.0;
      err = 0.0;
      for(int i=0; i<n; i++)
      {
        sum = 0.0;
        for(int j=0; j<n; j++)
        {
	  sum = sum + AA[i][j]*x[j];
        } // j
        sumerr = sumerr + Math.abs(yy[i]-sum);
        err = Math.max(err,Math.abs(x[i]-x_soln[i])); 
      } // i
      System.out.println("random n="+n+", sumerr="+sumerr+", maxerr="+err);

      // build Hilbert test case
      if(n==4) System.out.println("Hilbert matrix and solution");
      for(int i=0; i<n; i++) // all subscripts zero based
      {
        for(int j=0; j<n; j++)
        {
	  A[i][j] = 1.0/(double)(i+j+2);
          AA[i][j] = A[i][j];
	  if(n==4) System.out.println("A["+i+"]["+j+"]="+A[i][j]);
        }
        x_soln[i] = (double)(i+1);
      }

      for(int i=0; i<n; i++)
      {
        sum = 0.0;
        for(int j=0; j<n; j++)
        {
	  sum = sum + A[i][j]*x_soln[j];
        } // j
        y[i] = sum;
        yy[i] = y[i];
        if(n==4) System.out.println("y["+i+"]="+y[i]+", x="+x_soln[i]);
      } // i
   
      new crout(n, A, y, x);

      if(n==4)
      {
        System.out.println("A destroyed");
        for(int i=0; i<n; i++)
	  for(int j=0; j<n; j++)
            System.out.println("A["+i+"]["+j+"]="+A[i][j]);
      }
      // check accuracy
      sumerr = 0.0;
      for(int i=0; i<n; i++)
      {
        sum = 0.0;
        for(int j=0; j<n; j++)
        {
	  sum = sum + AA[i][j]*x[j];
        } // j
        sumerr = sumerr + Math.abs(yy[i]-sum);
        err = Math.max(err,Math.abs(x[i]-x_soln[i])); 
      } // i
      System.out.println("Hilbert n="+n+", sumerr="+sumerr+", maxerr="+err);
      if(n==4)
      {
        System.out.println("Returned solution, sumerr="+sumerr);

        for(int i=0; i<n; i++)
        { 
          System.out.println("x["+i+"]="+x[i]+", err="+(x[i]-x_soln[i]));
        } // i
      } // end if;
    } // n
  } // test_crout

  public static void main (String[] args)
  {
      new test_crout();
  }
} // end class test_crout