/*********************************************************************************** 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 # include # include # 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 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; }