/***********************************************************************************
Treibergs                                                              Apr. 12, 2006 

Solution to HW7. Invert a matrix using Gaussian Elimination

Enter  n  and an  n x n   matrix. Builds an  n x 2n  augmented matrix  (A,I). Then 
it finds the row-echelon and reduced-row-echelon form of the augmented matrix.  If
it is  (I,B)  then  B is the desired inverse. Otherwise the input matrix was 
singular


Programming Notes

This is a modification of the code for doing Gaussian elimination.  enterm
is modified to return an  n x 2n  matrix.  printma  and  printmb  print the first half
or the second half of the matrix. After the reduction to row echelon form, we can tell
if the matrix is nonsingular if the lower right hand corner,  a[n-1][n-1],  is one 
and not zero. Then we call reducem to get the reduced row echelon form.

gausselim.c
*************************************************************************************/
      
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define  MAXDIM  15
# define  MAXDIM2 30
# define  TOL     0.5e-7


typedef struct
{
    double m[MAXDIM][MAXDIM2];
    int size;
} matrix;

matrix
enterm ( void );

matrix
addmurom ( matrix, int, double, int );

matrix
swaprom (matrix, int, int );

matrix
mulrom ( matrix, int, double );

matrix
reducem ( matrix );


void 
printma ( matrix );

void 
printmb ( matrix );


int 
main ( void ) 
{
    double pivot;
    int i, j, k, ipivot, lasti=0;
    matrix a;
    
    a = enterm (); 
    printf("\n\nInitial Matrix =");
    printma ( a );
   
    
    for( j=0 ; j < 2*a.size; j++ )
    {
        ipivot=-1;
        pivot=0.0;
        
/* Partial pivoting step: find the largest in pivotal col below last pivotal row */ 
       
        for ( i = lasti; i< a.size; i++ )
        {
            if( fabs ( a.m[i][j] )  > fabs ( pivot ) + TOL )
            {
                pivot = a.m[i][j];
                ipivot = i;
            }
        }
        
/* Skip to next col if all zeros. Else swap in pivoti row and eliminate other rows */
        
        if ( ipivot != -1 )
        {
            if (ipivot != lasti)
                a = swaprom ( a, ipivot, lasti);
                
        
            if ( pivot != 1.0 )    
                a = mulrom ( a, lasti, 1.0/pivot );
                
        
            for ( i = lasti  +  1; i < a.size; i++ )
            {
                if ( fabs ( a.m[i][j] )  >=  TOL )
                    a = addmurom ( a, i, -a.m[i][j], lasti );

                a.m[i][j] = 0.0;
                
            }
            lasti++;
        }
    }
    

    i = a.size - 1;
    if(  a.m[i][i]  < 0.5 )
    {
            printf("\nThe matrix is singular.\n");
    }
    else
    {
        a = reducem ( a );
        printf("\nInverse Matrix = ");
        printmb ( a );
    }
     return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
enterm ( void )
{
      int i=0, j,  ir, ic;
      matrix g;
      
      while ( i++ <= 10 )
      {
          printf ( "\nEnter the size  : " );
          scanf ( "%d", &j );
          
          if ( (1<=j) && (j <= MAXDIM)  )
          {
               for ( ir=0; ir< j; ir++ )
               {
                   printf ( "Enter row  a(%d, * ) :", ir+1 );
                   for ( ic = 0; ic < j; ic++ )
                   {
                       g.m[ir][ic+j] = 0.0;
                       scanf ( "%lf", &g.m[ir][ic] );
                   }
               g.m[ir][ir+j] = 1.0;
               }
               g.size = j;
               return g;          
          }
          else
                printf ( "Bad input: size must be between %d and %d.\n", 1, MAXDIM );

      }
      exit ( 0 );
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
mulrom ( matrix a, int i, double c )
{
      int j;
      for ( j=0; j < 2*a.size; j++ )
          a.m[i][j] = c * a.m[i][j];
      return a;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
addmurom( matrix a, int i, double c,  int j )
{
      int k;
      for ( k = 0; k < 2*a.size; k++ )
          a.m[i][k] += c * a.m[j][k];
      return a;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
swaprom (matrix a, int i, int j)
{
      int k;
      double temp;
      for ( k=0; k < 2*a.size; k++ )
         {
              temp = a.m[i][k];
              a.m[i][k] = a.m[j][k];
              a.m[j][k] = temp;
          }
      return  a;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

void 
printma ( matrix a )
{
    int i, j;
    printf ( "\n\n" );
    for ( i = 0; i < a.size; i++ )
    {
        for ( j = 0; j < a.size; j++ )
            printf ( "\t%f", a.m[i][j] );
        printf ( "\n\n" );
    }
    return;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

void 
printmb ( matrix a )
{
    int i, j;
    printf ( "\n\n" );
    for ( i = 0; i < a.size; i++ )
    {
        for ( j = a.size; j < 2*a.size; j++ )
            printf ( "\t%f", a.m[i][j] );
        printf ( "\n\n" );
    }
    return;
}


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
reducem ( matrix a )
{
    int i, j, ilast = 1, jlast = 1;
    if ( a.size > 1)
    {
        while ( (ilast < a.size) && (jlast < 2*a.size) )
        {
            if ( fabs ( a.m[ilast][jlast] ) < TOL )
                 jlast++;
            else
            {
                  for ( i=0; i<ilast; i++)
                  {
                      if ( fabs(a.m[i][jlast]) > TOL )
                      {
                           for ( j = jlast + 1; j< 2*a.size; j++ )
                                a.m[i][j] -= a.m[i][jlast] * a.m[ilast][j];
                      }
                      a.m[i][jlast] = 0.0;
                  }
                  ilast++;
                  jlast++;
            }
        } 
    }
    return a;
}