/* Linear equation solver using Gaussian elimination.

   The coefficients are stored as

   m[0][0],m[0][1],...,m[0][n-1]
   m[1][0],m[1][1],...,m[1][n-1]
                ...
   m[n-1][0],m[n-1][1],...,m[n-1][n-1]

   The constants are stored as

   m[0][n],m[1][n],...,m[n-1][n]
*/

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

#define MAXSIZE 15

void lin_solve( float m[ ][ MAXSIZE + 1 ], float x[ ], int n );

main()
{
   float workspace[ MAXSIZE ][ MAXSIZE + 1 ], 
         x[ MAXSIZE ], factor, temp;
   int i, j, n;

   printf( "Number of equations? " );
   scanf( "%d", &n );
   printf( "\n\nEnter coefficients by rows--\n" );
   for ( i = 0; i < n; i++ )
      for ( j = 0; j < n; j++ )
         scanf( "%f", &workspace[ i ][ j ] );

   printf( "\n\nEnter constants--\n" );
   for ( i = 0; i < n; i++ )
      scanf( "%f", &workspace[ i ][ n ] );

   /* Solve using Gaussian elimination */
   lin_solve( workspace, x, n );
   printf( "\n\nSolution--\n" );
   for ( i = 0; i < n; i++ )
      printf( "\tx[ %d ] = %f\n", i, x[ i ] );

   return EXIT_SUCCESS;
}

/* solve using Gaussian elimination */
void lin_solve( float m[ ][ MAXSIZE + 1 ], float x[ ], int n )
{
   int i, j, k, pivot;
   float factor, temp;

   for ( i = 0; i < n; i++ ) {
      if ( m[ i ][ i ] == 0.0 ) {
         pivot = 0;
         /* find next nonzero entry in col i */
         for ( j = i + 1; j < n; j++ )
            if ( m[ j ][ i ] != 0.0 ) {
               pivot = j;
               break;
            }
         /* if no nonzero entry in col i, system is singular */
         if ( pivot == 0 ) {
            printf( "System is singular\n" );
            exit( EXIT_SUCCESS );
         }
         /* swap so m[i][i] != 0 */
         for ( j = 0; j < n + 1; j++ ) {
            temp = m[ i ][ j ];
            m[ i ][ j ] = m[ pivot ][ j ];
            m[ pivot ][ j ] = temp;
         }
      }

      /* make col i, row j >= i + 1, zero */
      for ( j = i + 1; j < n; j++ ) {
         factor = -m[ j ][ i ] / m[ i ][ i ];
         for ( k = i; k < n + 1; k++ )
            m[ j ][ k ] += factor * m[ i ][ k ];
      }
   }

   /* solve for unknowns */
   x[ n - 1 ] = m[ n - 1][ n ] / m[ n - 1 ][ n - 1];
   for ( j = n - 2; j >= 0; j-- ) {
      x[ j ] = m[ j ][ n ];
      for ( k = j + 1; k < n; k++ )
         x[ j ] -= m[ j ][ k ] * x[ k ];
      x[ j ] /= m[ j ][ j ];
   }
}