Augmentation & Gauss-Jordan Elimination

[al]

I tried to reduce this to a simple 2 functions for solving a set of linear simultaneous equations in c.  First augment the equations with the results in one ( N X N+1) matrix, then run it through the Gauss-Jordan algorithm to be reduced.  I explained the Gauss Jordan algorithm almost line by line with printing too.  Here is what you need :

Code:

#define MAX_ROWS   100  
#define MAX_COL   100

//Matrix main struct
struct matrix
{
   double m[MAX_ROWS][MAX_COL];
   int r;
   int c;
};

// Macros

////////////////////////////////////////////////////////////////////////////////////////////////
// PRINTMATRIXP Macro can be called to print a matrix of type pointer to "struct matrix" to
// standard output
//
// Macro PRINTMATRIXP take the following arguments
//         1- MAT:   A POINTER Matrix of type struct matrix
#define PRINTMATRIXP(MAT) {int i,j;for(i=0;i<MAT->r;i++){for(j=0;j<MAT->c;j++){printf("%10.5f\t",MAT->m[j]);}printf("\n");}printf("\n");}
////////////////////////////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////////////////////////////
// PRINTMATRIX Macro can be called to print a matrix of type "struct matrix" to standard output
//
// Macro PRINTMATRIX take the following arguments
//         1- MAT:   A NON-POINTER Matrix of type struct matrix
#define PRINTMATRIX(MAT) {int i,j;for(i=0;i<MAT.r;i++){for(j=0;j<MAT.c;j++){printf("%10.5f\t",MAT.m[j]);}printf("\n");}printf("\n");}
////////////////////////////////////////////////////////////////////////////////////////////////



////////////////////////////////////////////////////////////////////////////////////////////////
// INITMATRIX is a Macro for inititializing and priming a NON-POINTER matrix structure with the
// your values from a 2 diminsional array you supply
//
// Macro INITMATRIX take the following arguments
//         1-   TMP: A 2 diminsional array matrix with values set
//         2-   MAT: The matrix you declated in you program (non-pointer)
//         3-   ROW: The number of rows in the matrix  
//         4-   COL: The number of columns in the matrix
#define INITMATRIX(TMP,MAT,ROW,COL) {MAT.r=ROW;MAT.c=COL;COPYARRAY(TMP,MAT);}
////////////////////////////////////////////////////////////////////////////////////////////////

// Called with 2 Matrices and will return a one augmented matrix
bool MatAugment2Mats(matrix *mat1, matrix *augmat, matrix *res)
{
   int i=0, j=0, k=0, l=0;

   if(mat1->r!=augmat->r)
      return false;

   // Initialize result matrix
   res->c=mat1->c+augmat->c;
   res->r = mat1->r;
   for(i=0;i<res->r;i++)
      for(j=0;j<res->c;j++)
         res->m[j] = 0;

   for(i=0;i<res->r;i++)
   {
      for(j=0;j<res->c;j++)
      {
         if(j<mat1->c)
            res->m[j] =mat1->m[j];
         else
            res->m[j] =augmat->m[j-mat1->c];

      }
   }
   return true;
}


// called with one augmented matrix of equations and results and will reduce them to identity matrix and return the result matrix in NX1
bool MatGaussReduce(matrix *Amat,  matrix *res)
{
int i=0, j=0, k=0;
int max =0;
double temp=0.0;

res->c = 1;
res->r = Amat->r;
for (i=0;i<res->r;i++)
             res->m[i][0] = 0.0;

cout<<"For each column:"<<endl;

// for each column
for (i = 0; i < Amat->r; ++i)
{
// find the pivot element; look for the largest element
// and move up the reduction triangle
max = i; // assume it is the current row
for (j = i + 1; j < Amat->r; ++j) // scan the rows below for a larger element
if (Amat->m[j][i] > Amat->m[max][i]) // if found get the element row 'max'
max = j;

cout<<"column "<<i+1<<" Pivot = "<<Amat->m[max][i]<<endl;

// Reorder the matrix so the pivot is at the diagonal
for (j = 0; j < Amat->r + 1; ++j)
{
temp = Amat->m[max][j];
Amat->m[max][j] = Amat->m[i][j];
Amat->m[i][j] = temp;

}
cout<<"After reordering  .."<<endl;
PRINTMATRIXP(Amat);

cout<<"For each column "<<endl;
for (j = Amat->r; j >= i; --j)
{
cout<<"For column "<<j+1<<endl;
for (k = i + 1; k < Amat->r; ++k)
{
cout<<"Substract "<< Amat->m[k][i]/Amat->m[i][i] * Amat->m[i][j]<<" from row "<<k+1<<" column "<<j+1<<endl;
Amat->m[k][j] -= Amat->m[k][i]/Amat->m[i][i] * Amat->m[i][j];
PRINTMATRIXP(Amat);

}
cout <<"Now we have "<<endl;
PRINTMATRIXP(Amat);
}

}

// check that the bottom factor is not Zero (approximate to 0.00001), if it is, then we have no solution
if((Amat->m[Amat->r-1][Amat->c-2] <0.00001) && (Amat->m[Amat->r-1][Amat->c-2] >-0.00001))
return false; // there is NO solution, return false

cout <<"Now back substitution from the bottom "<<endl;
for (i = Amat->r - 1; i >= 0; --i)
{

Amat->m[i][Amat->r] = Amat->m[i][Amat->r] / Amat->m[i][i];
Amat->m[i][i] = 1;

for (j = i - 1; j >= 0; --j)
{

Amat->m[j][Amat->r] -= Amat->m[j][i] * Amat->m[i][Amat->r];
Amat->m[j][i] = 0;

}
cout<<"For row "<<i+1<<endl;
PRINTMATRIXP(Amat);
}

// copy the last column in the augmented matrix to the result matrix
for (i=0;i<res->r;i++)
res->m[i][0] = Amat->m[i][Amat->r];

return true;
}


int main()
{
   matrix   res, res2,Gmat, Gaug;

   //--------------------------------------------------------------------------------------------------
   // GAUSS Elimination
   // You can either create 2 matrices and augment them then solve the augmented
   // matrix as in the example below **OR** build a matrix already augmented with the results (not shown)
   // First: you need to create the linear equations matrix
   double tmp4[][3]=
   {   {1,1,2},
      {-1,-2,3},
      {3,-7,4}
   };
   // Second: Call the Init Matrix Macro to prime it with the temporary matrix
   INITMATRIX(tmp4,Gmat,3,3);

   // Create the result matrix
   double tmp5[][1]=
   {   {8},
      {1},
      {10}
   };
   INITMATRIX(tmp5,Gaug,3,1);

   // Augment the 2 matrices into 1 using the MatAugment2Mats() function
   if(!MatAugment2Mats(&Gmat, &Gaug, &res))
   {
      printf("Error: The 2 matrices cannot be augmented due to row numbers mismatch!! \n");
   }
   else
   {
      // Print the full aumented matrix before the elimination call
      printf("Gauss Elimination:: Printing ORIGINAL Matrix (%dX%d)\n",res.r,res.c);
      printf("Original set of equations:\n");
      PRINTMATRIX(res);  

      // if there is a solution find it and print it
      if(MatGaussReduce(&res,&res2))
      {
         printf("Gauss Elimination:: Printing RESULT Matrix (%dX%d)\n",res.r,res.c);
         PRINTMATRIX(res);  
      }
      else
      {
         // MatGaussReduce() can return a false if there is inconsistencies OR the set of equations
         // cannot be solved (has infinate solutions/indeterministic )
         printf("Error: Either bad matrix or the matrix is indeterministic!! \n");
         printf("Reduced Matrix:\n");
         PRINTMATRIX(res);
      }
   }



}

[Seedhom]

Here is a common use of the of this elimination algorithm I used for line intersection in 2D plane.  Note that the hardest part of this is populating the array with the equations coefficients for each line as in Ax + By = C.

Code:

/**********************************************************************
* Function get_intersect_pt2
* IN : 2 lines end points --> Line 1: (l1x1, l1y1) to (l1x2, l1y2)
* Line 2: (l2x1, l2y1) to (l2x2, l2y2)
* OUT : Intersection point (x, y) by reference
*
**********************************************************************/
bool get_intersect_pt2(double l1x1, double l1y1, double l1x2, double l1y2,
 double l2x1, double l2y1, double l2x2, double l2y2,
 double *x, double *y)
{
// First: you need to create the linear equations matrix
matrix Gmat,  // The matrix to be filled by the macro from the 2D array
res;   // The matrix to hold the results returned to us
double tmp[2][3];   // our tmp 2D array with the lines equations in the format:
    // line1: A1x + B1y = C1
    // line2: A2x + B2y = C2
    // filled temp array will be
    // tmp[][] = [A1,B1,C1
    //   A2,B2,C2 ]


tmp[0][0] = l1y2 - l1y1; //A1
tmp[0][1] = l1x1 - l1x2; //B1
tmp[0][2] = tmp[0][0] * l1x1  + tmp[0][1] * l1y1; //C1

tmp[1][0] = l2y2 - l2y1; //A2
tmp[1][1] = l2x1 - l2x2; //B2
tmp[1][2] = tmp[1][0] * l2x1  + tmp[1][1] * l2y1; //C2

// Call the Init Matrix Macro to prime it with the temporary matrix
// and fill in Gmat
INITMATRIX(tmp,Gmat,2,3);

// Reduce it
if(MatGaussReduce(&Gmat,&res))
{
*x = res.m[0][0];
*y = res.m[1][0];
return true;  // True if we find intersection
}
else
{
return false; // Lines are parallel
}

}

----------------------------
Al Sabawi