/***********************************************************************************
Treibergs Mar. 19, 2006
Gaussian Elimination Determinant
Enter a square matrix. Finds the row-echelon form of the matrix. Then gets det.
Idea is that each of the basic row operations
R1 = divide a row by c. det(R1)=1/c.
R2 = swap two rows. det(R2)=-1.
R3 = add a multiple of ith row to jth row. det(R3)=1.
Then the Row-echelon matrix, U, is upper triangular with 1's or 0's on the diagonal.
Hence det(U) = U[n][n] = the product of its diagonals.
Doing a sequence or row ops yields
R * R * ... R * A = U.
where R is one of the three operations. Thus, taking det,
det(R) det(R) ... det(R) det(A) = det(U)
or
det(A) = det(U) . (1/det(R)) . (1/det(R) . ... . (1/det(R))
We collect the product on the right in the variable detfac.
gauelidet.c
*************************************************************************************/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define MAXDIM 20
# define TOL 0.5e-7
typedef struct
{
double m[MAXDIM][MAXDIM];
int size;
} matrix;
matrix
enterm ( void );
matrix
addmurom ( matrix, int i, double c, int j );
matrix
swaprom (matrix, int i, int j );
matrix
mulrom ( matrix a, int i, double c );
void
printm ( matrix );
int
main ( void )
{
double pivot, detfac=1.0;
int i, j, k, ipivot, lasti=0, lastj=0;
matrix a;
a = enterm ();
printf("\n\nInitial Matrix");
printm ( a );
for( j=lastj; j < 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);
detfac = -detfac;
}
if ( pivot != 1.0 )
{
a = mulrom ( a, lasti, 1.0/pivot );
detfac *= pivot;
}
for ( i = lasti + 1; i < a.size; i++ )
{
if ( fabs ( a.m[i][j] ) < TOL )
a.m[i][j] = 0.0;
else
a = addmurom ( a, i, -a.m[i][j], lasti );
}
lasti++;
}
}
printf("\nRow-Echelon Form of the Matrix");
printm ( a );
j = a.size - 1;
printf ( "\nThe determinant = %f\n", a.m[j][j]*detfac );
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
enterm ( void )
{
int i=0, j, ir, ic;
matrix g;
while ( i++ <= 10 )
{
printf ( "\nEnter the size = number of rows = number of columns : " );
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++ )
scanf ( "%lf", &g.m[ir][ic] );
}
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 < 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 < 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 < a.size; k++ )
{
temp = a.m[i][k];
a.m[i][k] = a.m[j][k];
a.m[j][k] = temp;
}
return a;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm ( 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;
}
There are more efficient ways than row expansion to compute the determinant. One of these is Gaussian Elimination with partial pivoting. By keeping track of row multiplications and swaps, the determinant of the matrix may be computed in the Gaussian Elimination program by just a couple of extra lines of code.
It is impressive to see the recursions do row expansions. However, the number of operations grows like O(n!) in the dimension for expansion by cofactors. Gaussian elimination takes O(n3).
/***********************************************************************************
Treibergs Mar. 19, 2006
Determinants via Column Expansion
This is an extremely inefficient way to compute determinants. However it illustrates
structures, arrays and recursive functions.
The determinant may be expanded on any row or column. let a be a nxn matrix.
If we let strike(a,i,j) be the (n-1)x(n-1) matrix gotten from a by crossing out
the ith row and jth column. Then, e.g., the determinat equals its expansion along
the jth column, given by
det(a) = (-1)^(1+j) a[1][j] det( strike(a,1,j) ) + . . . +
(-1)^(i+j) a[i][j] det( strike(a,i,j) ) + . . . +
(-1)^(n+j) a[n][j] det( strike(a,n,j) )
Note that the complexity increases like n! In this program, we have used the last
column of a to expand the det.
determ.c
*************************************************************************************/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define MAXDIM 10
typedef struct
{
double m[MAXDIM][MAXDIM];
int size;
} matrix;
double
determ( matrix );
void
printm ( matrix );
matrix
enterm ( void );
matrix
strike ( matrix, int );
int
main ( void )
{
int i,j,k;
matrix a;
a = enterm ();
printf("\n\nInitial Matrix");
printm ( a );
printf(" Determinant = %f\n", determ(a) );
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
strike( matrix a, int i) /* strike ith row and last column */
{
int j, k;
for ( j = i + 1; j < a.size; j++ )
{
for ( k = 0; k < a.size - 1; k++ )
a.m[j-1][k] = a.m[j][k];
}
a.size = a.size - 1;
return a;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
double
determ ( matrix a )
{
double sum=0.0, sign = 1.0;
int i,n;
n = a.size;
if( n == 1 )
return a.m[0][0];
else
{
for ( i = n - 1 ; i >= 0; i-- )
{
sum = sum + sign * a.m[i][n-1] * determ( strike( a, i) );
sign = -sign;
}
return sum;
}
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
enterm ( void )
{
int i=0, j, k, ir, ic;
matrix g;
while ( i++ <= 10 )
{
printf ( "\nEnter the dimension of a square matrix, n : " );
scanf ( "%d", &j );
if ( (1 <=j) && (j <= MAXDIM) )
{
for ( ir=0; ir < j; ir++ )
{
for ( ic = 0; ic < j; ic++ )
{
printf ( "Enter a(%d,%d) :", ir+1, ic+1 );
scanf ( "%lf", &g.m[ir][ic] );
}
}
g.size = j;
return g;
}
else
printf ( "Bad input: size must be between %d and %d.\n", 1, MAXDIM );
}
exit ( 0 );
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm ( 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;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/* Milicic
Determinants using cofactor expansion
*/
#include <stdio.h>
#define dim 20
typedef struct {
float coeff[dim][dim];
int size;
} matrix;
matrix cofactor( matrix A, int k) {
matrix B;
int i, j ;
B.size = A.size - 1;
for (j = 0 ; j < B.size ; j++ )
for (i = 0 ; i < B.size ; i++ ){
if ( i < k )
B.coeff[j][i] = A.coeff[j+1][i];
else
B.coeff[j][i] = A.coeff[j+1][i+1];
}
return(B);
}
float Det( matrix A ){
matrix B;
float det;
int i;
det = 0;
if (A.size == 1)
det = A.coeff[0][0];
else
for (i = 0 ; i < A.size ; i++ ){
B = cofactor(A,i);
if (i%2 == 0)
det = det + A.coeff[0][i] * Det(B);
else
det = det - A.coeff[0][i] * Det(B);
}
return(det);
}
int main(){
matrix A;
int i, j;
for ( i = 0 ; i < dim ; i++ ){
for ( j =0 ; j < dim ; j++ )
A.coeff[i][j] = 0;
}
printf("The size of the determinant is: ");
scanf("%d",&A.size);
for ( i = 0 ; i < A.size ; i++ ){
printf("Enter the %dth row of the matrix: ",i+1);
for ( j =0 ; j < A.size ; j++ )
scanf("%f",&A.coeff[i][j]);
printf("\n");
}
printf("The determinant is %f\n",Det(A));
}
/* Milicic
Cramer's Rule
cramer.c */
#include <stdio.h>
#include <math.h>
#define dim 20
#define epsilon 0.000001
typedef struct {
double coeff[dim][dim];
int size;
} matrix;
matrix cofactor( matrix A, int k) {
matrix B;
int i, j ;
B.size = A.size - 1;
for (j = 0 ; j < B.size ; j++ )
for (i = 0 ; i < B.size ; i++ ){
if ( i < k )
B.coeff[j][i] = A.coeff[j+1][i];
else
B.coeff[j][i] = A.coeff[j+1][i+1];
}
return(B);
}
double Det( matrix A ){
matrix B;
double det;
int i;
det = 0;
if (A.size == 1)
det = A.coeff[0][0];
else
for (i = 0 ; i < A.size ; i++ ){
B = cofactor(A,i);
if (i%2 == 0)
det = det + A.coeff[0][i] * Det(B);
else
det = det - A.coeff[0][i] * Det(B);
}
return(det);
}
int main(){
matrix A, C;
double B[dim], X[dim];
int i, j, k, n;
double detsys;
for ( i = 0 ; i < dim ; i++ ){
for ( j =0 ; j < dim ; j++ )
A.coeff[i][j] = 0;
}
printf("The size of the system is: ");
scanf("%d",&n);
A.size = n;
for ( i = 0 ; i < n ; i++ ){
printf("Enter the %dth row of the matrix: ",i+1);
for ( j =0 ; j < n ; j++ )
scanf("%lf",&A.coeff[i][j]);
printf("\n");
}
detsys = Det(A);
if (fabs(detsys) < epsilon)
printf("The system seems to be singular.");
else {
printf("Enter the last column of the extended matrix: ");
for ( i =0 ; i < n ; i++ )
scanf("%lf",&B[i]);
printf("\n");
C.size = n;
for(k = 0; k < n; k++ ){
for (i = 0; i < dim; i++ )
for ( j = 0 ; j < dim ; j++ )
C.coeff[i][j] = A.coeff[i][j];
for (i = 0; i < n; i++ )
C.coeff[i][k] = B[i];
X[k] = Det(C)/detsys;
}
printf("The solution of the system is: \n");
for(k = 0; k < n; k++ )
printf("%f \t",X[k]);
};
printf("\n");
}