# Thirteenth Week Examples --   Lovely Theory  --  Lousy Computation

Following D. Milicic Thirteenth Week Examples and Fourteenth Week Examples.

## Gaussian Elimination method to find the Determinant

Here is one of the standard ways to compute determinants for general matrices. The other is using LU decompositions. In doing Gaussian Elimination, the determinant is the product of row factors and swaps. This is a modification of our Gaussian Elimination code.
```/***********************************************************************************
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;
}

```

## Row Expansion for the Determinant

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));
}

```

## Cramer's Rule

```
/* 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");
}

```