#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct { int height; double * entry; }
  VECTOR;

typedef struct { int height, width ; double ** entry; }
  MATRIX;

	/*
	 * zero_matrix returns m x n matrix all zero
	 */

MATRIX * zero_matrix ( int m, int n )
{
  MATRIX * mat = calloc ( 1, sizeof ( MATRIX ) );

  printf("zero_matrix TO BE COMPLETED\n");

	/***********************
	 *** TO BE COMPLETED ***
	 ***********************/

  return mat;
}

	/*
	 * zero_vector returns vector of height n, all zero
	 */

VECTOR * zero_vector ( int n )
{
  VECTOR * v = ( VECTOR * ) calloc(1, sizeof ( VECTOR ) );

  printf("zero_vector TO BE COMPLETED\n");

	/***********************
	 *** TO BE COMPLETED ***
	 ***********************/

  return v;
}

	/*
	 * read_aug_matrix interprets the
	 * data in standard input as an
	 * m x (m+1) augmented matrix.
	 * If the wrong size, it aborts the program.
	 * 
 	 * It creates an n x n matrix and an n-vector
	 * returning them through a and b.
	 *
	 * Usage:
	 * MATRIX * a; VECTOR * b;
	 * read_aug_matrix ( &a, &b );
	 */

void read_aug_matrix ( MATRIX **a, VECTOR **b )
{
  int m, n, i, j;
  MATRIX * mat;
  VECTOR * vec;
  double x;

  scanf ( "%d %d", &m, &n );
  if ( n != m+1 )
  {
    printf("aug matrix %d rows %d columns, wrong shape, abort\n",
		m, n);
    exit(-1);
  }
    
  mat = zero_matrix ( m, n-1 );
  vec = zero_vector ( m );

  for (i=0; i<m; ++i)
  for (j=0; j<n; ++j)
  {
    scanf ( "%lf", & ( x ) );
    if ( j < n-1 )
      mat->entry[i][j] = x;
    else
      vec->entry[i] = x;
  }

  *a = mat;
  *b = vec;
}

void print( MATRIX * a, VECTOR * b )
                /* print the augmented array */
{
  int i,j;
  for ( i=0; i<a->height; ++i )
  {
    for  (j=0; j<a->width; ++j )
      printf ( " %10.5f", a->entry[i][j] );
    printf ( "   %10.5f", b->entry[i]    );
    printf ("\n");
  }
}

	/*
	 * Three EROs: subtract, scale, swap.
	 */

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

void scale ( MATRIX *a, VECTOR *b, int i, double s )
	/* multiply row i by s */
{
  int c;
  for ( c=0; c<a->width; ++c )
    a->entry[i][c] = a->entry[i][c] * s;

  b->entry[i] *= s;
}

void swap ( MATRIX * a, VECTOR * b, int i, int j )
 	/* swap rows i and j */
{
   double temp;
    int c;
  for ( c= 0; c<a->width; ++c )
  {
    temp = a->entry[i][c];
    a->entry[i][c] = a->entry[j][c];
    a->entry[j][c] = temp;
  }

  temp = b->entry[i];
  b->entry[i] = b->entry[j];
  b->entry[j] = temp;
}

	/*
	 * reduce ( a, b, uu, yy )
	 * a is an n x n matrix, b an n-vector.
	 * Creates a copy u of a and y of b,
	 * and reduces them together, until
	 * u is upper triangular.
	 *
	 * Usage:
	 * MATRIX *a, *u; VECTOR * b, * y;
	 * .... set up a, then ....
	 * reduce ( a, b, &u, &y )
	 */
	 

reduce ( MATRIX * a, VECTOR * b, MATRIX ** uu, VECTOR ** yy )
	/* perform Gaussian elimination with partial pivoting */
{
  int i, j, k, ell, r, p;
  double s; /* scale */
  double runningmax, temp;
  MATRIX * u = zero_matrix ( a->height, a->width );
  VECTOR * y = zero_vector ( b->height );

  for (i=0; i<a->height; ++i)
  {
    for (j=0; j<a->width; ++j)
    {
      u->entry[i][j] = a->entry[i][j];
    }
    y->entry[i] = b->entry[i];
  }

  k = 0;
		/* k: no. of columns already cleared */
  for ( j=0; j<a->width; ++j )
		/* attempt to clear column j */
  {
    p = -1;
    runningmax = 0;
    for ( r = k; r < a->width; ++r )
		/* search column j, row k and below, 
		 * for entry max in absolute value */
    {
      temp = fabs ( u->entry[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 ( u, y, p, k );
       s = 1/u->entry[k][j];
	
        for ( ell =k+1; ell < u->height; ++ ell )
		/* subtract u[ell][j] times row k ... */
	  subtract (u, y,k, u->entry[ell][j] * s, ell );

        ++k;
     }
  }
  * uu = u;
  * yy = y;
}

VECTOR * back_substitute ( MATRIX * u, VECTOR * y )
{
  int i,j;
  double temp;
  VECTOR * x = zero_vector ( u->width );

  printf("back_substitute TO BE COMPLETED\n");

	/***********************
	 *** TO BE COMPLETED ***
	 ***********************/
  return x;
}

void print_vector( VECTOR * x )
{
  int i;

  for ( i=0; i<x->height; ++i )
    printf (" %f", x->entry[i] );
  printf ("\n");
}

main()
{
  MATRIX *a, *u;
  VECTOR *b, *y, *x;
  int n, w;
  int i,j;

  read_aug_matrix ( &a, &b );

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

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

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