Gaussian-Jordan Elimination

[hbcpp]

I have created a program that solves the linear system by Gaussian-Jordan Elimination method.

My program algorithm follows every step I was taught in school, swap rows, multiply rows by constant, (like in this video: https://www.youtube.com/watch?v=0fTSBIBD7Cs ), in the end my matrix has changed like in the video too.

I did the same steps as shown in this site too, https://matrix.reshish.com/

I would post the code here but I think it is too verbose and I haven't yet written any comments.

But I felt it was doing too many computations so I went online and found this https://www.tutorialspoint.com/cplusplus-program-to-implement-gauss-jordan-elimination

#include<iostream>
using namespace std;
int main() {
   int i,j,k,n; // declare variables and matrixes as
   input
   float a[10][10],b,x[10];
   printf("\nEnter the size of matrix: ");
   scanf("%d",&n);
   printf("\nEnter the elements of augmented matrix (rowwise):\ n");
   for(i=1; i<=n; i++) {
      for(j=1; j<=(n+1); j++) {
         cout << "A[" << i << ", " << j << " ]=";
         cin >> a[i][j];
      }
   }
   //to find the elements of diagonal matrix
   for(j=1; j<=n; j++) {
      for(i=1; i<=n; i++) {
         if(i!=j) {
            b=a[i][j]/a[j][j];
            for(k=1; k<=n+1; k++) { 
               a[i][k]=a[i][k]-b*a[j][k];
            }
         }
      }
   }
   cout<<"\nThe solution is:\n";
   for(i=1; i<=n; i++) {
      x[i]=a[i][n+1]/a[i][i];
      cout<<"x"<<i << "="<<x[i]<<" ";
   }
   return(0);
}

This one works greatly with extremely less computations but it fails when the linear system is one with an infinity of solutions.

I feel that the solution that I created would be a half-done work to determine the matrix inverse and determinant, even though it envolves a lot computations. But there must be a better solution than this.


So which one should I choose?? Is there a better way to solve it a linear system?

Thanks

[lastchance]

Try your code on the (perfectly invertible) system

( 0 1 ) (x) = (3)
( 1 0 ) (y) = (4)

You may need to read up on "partial pivoting".

You should also put in lines to detect singular matrices (which this one isn't).


BTW, arrays start at 0 in c/c++

[hbcpp]

arrays start at 0 in c/c++

yeah noticed that too, the code is the from the link

[hbcpp]

@lastchance

what do you mean by

Try your code on the (perfectly invertible) system

[lastchance]

Hello, @hbcpp.

I'm referring to the code in your post. If you try that particular matrix in the code then it will fail because it divides by a zero diagonal element. You need to swap some rows first (partial pivoting).

My point was that the code would fail for reasons other than the matrix being singular. This matrix isn't singular: it has determinant -1.

There are better methods than Gauss-Jordan for large matrices. However, the latter is a sound and (with partial pivoting) a relatively stable approach, which is good for checking more advanced methods.

Another advantage of Gauss-Jordan is that, even if the matrix is singular, it can, by reducing the system to row-echelon form, tell you the rank of the matrix (number of linearly independent rows).

Your earlier comment was also right - the technique can be used to find the determinant (product of the diagonal elements after reducing to upper triangular) and invert a matrix (doing the same row operations simultaneously on the identity matrix). Long ago these were some of the first programs that I ever wrote - but I was coding in ZX Spectrum Basic, not in the still-to-be-invented C++!

[hbcpp]

I have finished up the inverse and determinant code, and I did the same

...find the determinant (product of the diagonal elements after reducing to upper triangular) and invert a matrix (doing the same row operations simultaneously on the identity matrix...

The code, however, got super verbose I am afraid it will be a little challenging for others to read it even though it is simple to understand if you know the Gauss-Jordan method.

Here is the inverse code...

Matrix is a 2D vector containing objects of Fraction.

The functions zeroOutTheColumn, changePivotTo_one, pivotEqualTo_one_Found, pivotNot_zero_Found, swapRows are contained in a header file.

The code posted below is only to give an idea of what I meant, not to be completely understood as it is only a small part of a much larger program..

Just take a look at the code and have an idea how it works, and see if you can give me a feedback

Matrix Matrix::inverse()
{
    assert(is_square());

    Matrix mx = *this;
    Matrix inverse = Matrix::IDENTITY(rows);

    bool alternative_pivot_1_found;

    bool pivot_not_zero_found;

    int row_with_alternative_pivot;

    int row_with_pivot_not_zero;

    int pivot_row = 0;
    int pivot_col = 0;

// Gauss Elimination
    while (pivot_row < (rows - 1))
    {
        alternative_pivot_1_found = pivotEqualTo_one_Found (mx.data, pivot_row, pivot_col, rows, row_with_alternative_pivot);

        pivot_not_zero_found = pivotNot_zero_Found(mx.data, pivot_row, pivot_col, rows, row_with_pivot_not_zero);

        if (mx.data[ pivot_row ] [ pivot_col ] != 1 && alternative_pivot_1_found )
        {
            swapRows(inverse.data, pivot_row, row_with_pivot_not_zero, columns);
            swapRows(mx.data, pivot_row, row_with_pivot_not_zero, columns);
        }
        else if (mx.data[ pivot_row ] [ pivot_col ] == 0 && pivot_not_zero_found )
        {
            swapRows(inverse.data, pivot_row, row_with_pivot_not_zero, columns);
            swapRows(mx.data, pivot_row, row_with_pivot_not_zero, columns );
        }

        int col_dif_zero;

        firstNumberNot_zero(mx.data, pivot_row, columns, col_dif_zero);

        if (( mx.data[pivot_row] [col_dif_zero] ) != 1)
        {
            changePivotTo_one(inverse.data, pivot_row, columns, mx.data[ pivot_row ][ col_dif_zero ]);

            changePivotTo_one(mx.data, pivot_row, columns,
                              mx.data[ pivot_row ][ col_dif_zero ]);
        }

        int n = pivot_row + 1;

        while (n < rows)
        {
            Fraction constant = mx.data[ n ][ col_dif_zero ];

            zeroOutTheColumn(inverse.data, n, pivot_row, columns, constant);
            zeroOutTheColumn(mx.data, n, pivot_row, columns, constant);
            

            ++n;
        }

        ++pivot_row;
        ++pivot_col;
    }

//Jordan Elimination
    while(pivot_row > 0)
    {
        int col_dif_zero;

        firstNumberNot_zero(mx.data, pivot_row, mx.columns, col_dif_zero);

        if (( mx.data[pivot_row] [col_dif_zero] ) != 1)
        {
            changePivotTo_one(inverse.data, pivot_row, mx.columns, mx.data[ pivot_row ][ col_dif_zero ]);
            changePivotTo_one(mx.data, pivot_row, mx.columns, mx.data[ pivot_row ][ col_dif_zero ]);
        }

        int n = pivot_row - 1;

        while (n >= 0)
        {
            Fraction constant = mx.data[ n ][ col_dif_zero ];

         
            zeroOutTheColumn(inverse.data, n, pivot_row, mx.columns, constant);
            zeroOutTheColumn(mx.data, n, pivot_row, mx.columns, constant);
            

            --n;
        }
        --pivot_row;
    }

    return inverse;
}

[lastchance]

That looks like code for the "inverse", not the determinant.

Difficult to say much about the code without being able to test it. However, you appear to be using a "Fraction" class. Unless the matrix is small (and fairly "nice") I think the numerator and denominator of this would overflow quite rapidly - the method has a lot of multiplies and divides.

[hbcpp]

Sorry, I meant to write "inverse".

Difficult to say much about the code without being able to test it.

When I finish I will be uploading the whole project(codeblocks) to github, I will post the link here.

Unless the matrix is small (and fairly "nice") I think the numerator and denominator of this would overflow quite rapidly

I actually had a problem with overflow a while back.
In the Fraction, the numerator and denominator type is boost::multiprecision::cpp_int.

the method has a lot of multiplies and divides

Yeah that is the reason I am tryna find a better way

[jonnin]

inverse via append identity and solve?
you can also look at the algorithm for the pseudoinverse, which gives the true inverse if it has one with an interesting approach.
dets are best computed from the eigens, rather than the brute force approach. Getting the eigens, I used to use over-relax approach. But my work was in controls & stable matrices, didn't even have to normalize them they were well conditioned etc.

[hbcpp]

@jonnin thanks. I will look them up.

@ lastchance, I posted the full code for review, https://codereview.stackexchange.com/questions/236728/a-matrix-library-in-c

or you can download here, https://github.com/hbtalha/Matrix-Library