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*/ 


matrix Matread (int n,int m)
{               
  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");
  A=Matread(n,m);
  printf("\n");
  Matdump(A, n, m);
  A=Gauss(A, n, m);
  printf("\n");
  Matdump(A, n, m);
  A=Reduce(A, n, m);
  printf("\n");
  Matdump(A, n, m);
}