/* ------------------------------------------------------------------------ 
//  Use Gaussian elimination to solve M equations in M unknowns.  
//  Figures 19.13 Solving a system of linear equations. 
//          19.14  Reading the equations.
//          19.15  Print the equations.
//          19.16  Gaussian elimination.    
//          19.17  Pivoting.    
//          19.18  Scaling the kth equation.
//          19.19  Using the pivot equation to clear coefficients.  
//          19.20  Printing the answers.                 
// ------------------------------------------------------------------------ 
*/
#include "tools.h"

#define EPS .0001       /* Comparison tolerance for zero test. */

typedef double* equation;
typedef equation* linear_system;

linear_system read_system( int* unknowns );
void print_matrix( linear_system mat, int unknowns );
bool solve( linear_system mat, int unknowns );
bool swap_nonzero_row( linear_system mat, int unknowns, int k );
int  find_pivot( linear_system mat, int unknowns, int k );
void scale( equation eqn, int unknowns, int k );
void wipe( linear_system mat, int unknowns, int piv_row, int row );
void print_answers( linear_system mat, int unknowns );

/* --------------------------------------------------------------------- */ 
void main ( void )
{   
    int m;                          /* Number of equations. */  
    linear_system matrix;           /* Array of equation pointers. */
    bool solvable;

    banner();
    puts( " Linear Equations to Solve:" );
    matrix = read_system( &m );
    print_matrix( matrix, m );
    solvable = solve( matrix, m );

    puts( " Equations after Elimination:" );
    print_matrix( matrix, m );
    if (!solvable) puts( " Equations inconsistent or not independent." );
    else print_answers( matrix, m );
    free( matrix );
    bye();
}

/* ------------------------------------------------------------------------ 
**  Allocate and initialize memory for the system of equations. 
**      Allocate one array for each equation in the system.   
**      Store the pointer to each equation-area in one slot of mat.   
**      Read coefficients into each equation-array.   
*/ 
linear_system 
read_system( int* unknowns ) 
{   
    stream eq_in;
    linear_system mat;      /* Matrix is allocated and initialized here. */ 
    equation current;       /* Row currently receiving input. */
    int ok;                 /* Return code from fscanf(). */ 
    int col, k;             /* Loop counters. */     
    int m;                  /* Number of unknowns. */

    /* Open input file and read number of equations. ------------------- */ 
    eq_in = fopen( "gauss.in", "r" );
    if (eq_in == NULL) fatal( " Cannot open gauss.in" );
    fscanf( eq_in, "%i", &m );
    
    /* Allocate equation array. ---------------------------------------- */ 
    mat = (linear_system) malloc( m * sizeof(equation) );
    if (mat == NULL) fatal( " Cannot allocate enough memory." ); 
    
    /* Read in all the equations. -------------------------------------- */ 
    for (k = 0; k < m; ++k) {
        /* Allocate space for equation. -------------------------------- */  
        current = (equation) malloc( (m+1) * sizeof(double) );
        if (current == NULL) fatal( " Cannot allocate enough memory." ); 
        mat[k] = current;                     /* Attach array to matrix. */

        /* Read coefficients for one equation. ------------------------- */  
        for (col = 0; col <= m; ++col) { 
            ok = fscanf( eq_in, "%lg", &current[col] );
            if (ok < 1) fatal( " Read error or unexpected eof." ); 
        }
    }
    /* Return linear system and its size to main. ---------------------- */  
    fclose( eq_in );    /* Done with input; close stream. */ 
    *unknowns = m;      /* Return number of equations. */ 
    return mat;         /* Return matrix. */
}

/* ------------------------------------------------------------------------ 
**  Print the matrix, formatted so that each row is an equation.
*/
void 
print_matrix( linear_system mat, int unknowns )
{   
    int col, row;
    char op;            /* Char to print between terms in output equation. */ 
    equation current;   /* Pointer to new allocation. */

    printf( " ------------------------------------------------------------\n" );
    for (row = 0; row < unknowns; ++row) {  /* Print all equations. */
        current = mat[row];                 /* Next equation pointer. */
        op = '+';
        for (col = 0; col < unknowns; ++col ) {
            if (col == unknowns - 1) op = '='; 
            printf( "%8.3f * %c %c ", current[col], 'a'+ col, op );  
        }
        printf( "%8.3f\n", current[unknowns] ); 
    }   
    printf( " ------------------------------------------------------------\n" );
}

/* ------------------------------------------------------------------------ 
/*  This function performs the Gaussian elimination algorithm.
**      It uses row operations to reduce the first M columns of the matrix
**      to an identity matrix.  At that point, the last column contains the 
**      solution to the problem.
*/
bool 
solve( linear_system mat, int unknowns )
{   
    int row, k;
    bool solvable;
    
    for (k = 0; k < unknowns; ++k) {    
        solvable = swap_nonzero_row( mat, unknowns, k );
        if (!solvable) break;           /* Do all rows have 0 in column k? */
        scale( mat[k], unknowns, k );            /* Make a 1 in mat[k][k]. */
        /* print_matrix( mat, unknowns );                // For debugging. */
        
        /* Use the 1 in mat[k][k] to zero out the rest of column k. ------ */
        for (row = 0; row < unknowns; ++row)   
            if (k != row)  wipe( mat, unknowns, k, row );
    }
    return solvable;
}

/* ------------------------------------------------------------------------ 
**  Find pivot equation for the kth column. 
**      Pivot is the equation with the largest coefficient in column k.
**      Swap positions of pivot and the kth equation. 
**      If no non-zero coefficient exists, equations cannot be solved. 
*/
bool
swap_nonzero_row( linear_system mat, int unknowns, int k )
{   
    equation swap;                          /* Used to interchange rows. */
    int row = find_pivot( mat, unknowns, k );

    /* Coefficient k of pivot equation must be non-zero. --------------- */
    if (fabs( mat[row][k] ) < EPS ) return false;

    if (row != k) { swap = mat[row];  mat[row] = mat[k]; mat[k] = swap; }
    /* print_matrix( mat, unknowns );                  // For degugging. */
    return true;
}

/* ------------------------------------------------------------------------ 
**  Find the equation with the largest coefficient in column k. 
*/
int 
find_pivot( linear_system mat, int unknowns, int k )
{   
    double coefficient = fabs( mat[k][k] );       /* Largest value so far. */
    int big = k;                          /* Index of current coefficient. */
    int row;                                      /* Index of search loop. */

    for (row = k + 1; row < unknowns; ++row) {
        if (fabs( mat[row][k] ) > coefficient){ 
            coefficient = fabs( mat[row][k] );     /* A bigger coefficient. */
            big=row;                        /* Remember where it was found. */  
        }
    }
    return big;        /* Line number of equation with biggest coefficient. */
}

/* ------------------------------------------------------------------------ 
**  Scale kth equation in matrix to have 1 in kth place 
**      Assumes places 1..k-1 are already zero. 
*/
void 
scale( equation eqn, int unknowns, int k )
{   
    int col;                    /* loop index */
    double factor = eqn[k];     /* scale factor */
    for (col = k; col <= unknowns; ++col) 
        eqn[col] /= factor;
}

/* ------------------------------------------------------------------------ 
**  Clear the k1_th coefficient of k2_th equation by subtracting  
**      mat[k2][k1] * k1_th equation from k2_th equation.
**      Assumes that mat[k1][k1] == 1.   
*/
void 
wipe( linear_system mat, int unknowns, int piv_row, int row )
{   
    int col;
    double coefficient = mat[row][piv_row];

    for(col = piv_row; col <= unknowns; ++col) 
        mat[row][col] -= coefficient * mat[piv_row][col];
}

/* ------------------------------------------------------------------------ 
** Print each variable with its value (from last column of matrix).  
*/
void                        /* Print answers from last column of matrix. */
print_answers( linear_system mat, int unknowns )
{   
    int row;
    for (row = 0; row < unknowns; ++row)  
        printf( "    %c = %8.3f\n",  'a' + row, mat[row][unknowns] ); 
}