#include <stdio.h>
#include <stdlib.h>
#include <math.h>
	/*
	 * MUST BE COMPILED with -lm
	 */

void print( int n, double a[10][10], double b[10] )
                /* print the augmented array */
{
  int i,j;
  for ( i=0; i<n; ++i )
  {
    for  (j=0; j<n; ++j )
      printf ( " %10.5f", a[i][j] );
    printf ( "   %10.5f", b[i] );
    printf ("\n");
  }
}

void print_answer( int n, double x[10] )
{
  int i;

  printf("Answer (transposed) ");
  for ( i=0; i<n; ++i )
    printf (" %f", x[i] );
  printf ("\n");
}


void subtract ( int n, double a[10][10], double b[10], int i, double s, int j )
	/* subtract s times row i from row j */
{
  int c;
  for ( c=0; c<n; ++c )
    a[j][c] = a[j][c] - s * a[i][c];
  b[j] = b[j] - s * b[i];
}

void scale ( int n, double a[10][10], double b[10], int i, double s )
	/* multiply row i by s */
{
  int c;
  for ( c=0; c<n; ++c )
    a[i][c] = a[i][c] * s;

  b[i] *= s;
}

void swap ( int n, double a[10][10], double b[10], int i, int j )
 	/* swap rows i and j */
{
   double temp;
    int c;
  for ( c= 0; c<n; ++c )
  {
    temp = a[i][c];
    a[i][c] = a[j][c];
    a[j][c] = temp;
  }

  temp = b[i];
  b[i] = b[j];
  b[j] = temp;
}

reduce ( int n, double a[10][10], double b[10], double u[10][10],
		double y[10] )
	/* perform Gaussian elimination with partial pivoting */
{
  int i, j, k, ell, r, p;
  double s; /* scale */
  double runningmax, temp;

  for (i=0; i<n; ++i)
  {
    for (j=0; j<n; ++j)
    {
      u[i][j] = a[i][j];
    }
    y[i] = b[i];
  }

  k = 0;
		/* k: no. of columns already cleared */
  for ( j=0; j<n; ++j )
		/* attempt to clear column j */
  {
    p = -1;
    runningmax = 0;
    for ( r = k; r < n; ++r )
		/* search column j, row k and below, 
		 * for entry max in absolute value */
    {
      temp = fabs ( u[r][j] );
      if ( temp > runningmax )
      {
        p = r;
	runningmax = temp;
      }
    }
		/* either p == -1 or p>=k indexes a
		 * row such that u[p][j] is a pivot element */
    if ( p >= 0 )
    {
       swap (n, u, y, p, k );
       s = 1/u[k][j];
	
        for ( ell =k+1; ell < n; ++ ell )
		/* subtract u[ell][j] times row k ... */
	  subtract (n, u, y,k, u[ell][j] * s, ell );

        ++k;
     }
  }
}

void back_substitute ( int n, double u[10][10], double y[10], double x[10])
{
	/* TO BE FILLED IN */
}

main()
{
  double a[10][10], u[10][10];
  double b[10], y[10], x[10];
  int n, w;

	/* Fill in code here to scan in a and b together
	 * from input in the form of an augmented matrix such as
	 * 2 3
	 * 1 2 3
	 * 4 5 6
	 * 
	 * n=2 w=3 (n rows and w columns.  Expect w = n+1)
	 * a[0][0]=1 a[0][1]=2  b[0]=3 etcetera
	 */
	

  printf ("Input matrix\n");
  print( n, a, b );

  reduce ( n, a, b, u, y );

  printf ("Reduced matrix\n");
  print( n, u, y);
  back_substitute(n, u, y, x);
  print_answer(n, x);
}