## Linear Systems

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

```

### Gauss Elimination

```#include <stdio.h>
#include <math.h>

#define true  1
#define false 0

#define dim 20
#define epsilon 0.000001

typedef struct {
float coeff[dim][dim];
} matrix; /*define type matrix*/

{
int i, j;
matrix A;

printf("Enter matrix:\n");
for ( i=0 ; i< n ; i++ ) {
printf("\nEnter %d th row >> ",i+1);
for ( j=0 ; j < m ; j++ )
scanf("%f",&A.coeff[i][j]);
}
return(A);
}

/*Function Matdump dumps a matrix to the console */

void Matdump (matrix A, int n, int m)
{
int i, j;

for ( i=0 ; i < n ; i++ ) {
for ( j=0 ; j < m ; j++ )
printf("%f\t",A.coeff[i][j]);
printf("\n");
}
}
matrix Change (matrix A, int n, int m, int i, int j)
{
int k;
float x;

if (i != j) {
for ( k=0 ; k < m ; k++ ) {
x = A.coeff[i][k];
A.coeff[i][k] = A.coeff[j][k];
A.coeff[j][k] = x;
}
}
return(A);
}

/* Function MultRow multiplies elements of ith row in the matrix A by mult */

matrix MultRow (matrix A, int n, int m, int i, float mult)
{
int k;

for ( k=0 ; k < m ; k++ )
A.coeff[i][k] = mult * A.coeff[i][k];
return(A);
}

/* Function Rowmul multiplies elements of ith row in the matrix A by mult
and adds them to the jth row */

matrix Rowmul (matrix A, int n, int m, int i, int j, float mult)
{
int k;

for ( k=0 ; k < m ; k++ )
A.coeff[j][k] = A.coeff[j][k] + A.coeff[i][k] * mult;
return(A);
}

/* Function Find finds the element of maximal absolute value in
the jth column below the element A (i,j) */

int Find(matrix A, int n, int m, int i, int j)
{
int k,p;

k = i;
for ( p = i+1 ; p < n ; p++ ) {
if (fabs(A.coeff[k][j]) < fabs(A.coeff[p][j]))
k = p;
}
return(k);
}

matrix Gauss (matrix A, int n, int m)
{
int i, j, counter;
int zero;

counter = 0;
for ( j=0 ; j < m ; j++ ) {
zero = true;
for ( i = j - counter ; i < n ; i++ )
if (fabs(A.coeff[i][j]) > epsilon)
zero = false;
if (zero == true)
counter = counter + 1;
else {
i=Find(A, n, m, j - counter,j);
A=Change(A, n, m, i, j - counter);
A=MultRow(A, n, m, j - counter,1 / A.coeff[j - counter][j]);
for ( i = j - counter + 1 ; i < n ; i++ )
if (fabs(A.coeff[i][j]) > epsilon)
A=Rowmul(A, n, m, j - counter, i, -A.coeff[i][j]);
}
}
return(A);
}

/* This procedure reduces a matrix A in row-echelon for to reduced
row-echelon form */

matrix Reduce (matrix A, int n, int m)
{
int i, j, k, zero;

for ( i = n-1 ; i > 0 ; i-- ) {
zero = true;
for ( j = 0 ; j < n ; j++ )
if (fabs(A.coeff[i][j]) > epsilon)
zero = false;
if (zero == false) {
j = 0;
while (fabs(A.coeff[i][j]) < epsilon)
j++;
for ( k = 0 ; k < i ; k++ )
if (fabs(A.coeff[k][j]) > epsilon)
A=Rowmul(A, n, m, i, k, -A.coeff[k][j]);
}
}
return(A);
}

int main() {
int n,m;
matrix A;

printf("Enter  #rows, #columns >> ");
scanf("%d %d",&n, &m);
printf("\n");