title: Gauss-Jordan Elimination in C++
body: <p>I was using Gauss-Jordan elimination in C++ to solve a system of linear equations.
  Code works fine. Was wondering why Lines 1,2,3 in <code>void gauss()</code> can't be replaced by Line 4 (getting incorrect output after doing so)?</p>
  
  <pre><code>#include &lt;iostream&gt;
  using namespace std;
  class Gauss
  {
      float a[50][50];
      int n;
  public:
      void accept()
      {
          cout&lt;&lt;"Enter no. of variables: ";
          cin&gt;&gt;n;
          for(int i=0;i&lt;n;i++)
          {
              for(int j=0;j&lt;n+1;j++)
              {
                  if(j==n)
                      cout&lt;&lt;"Constant no."&lt;&lt;i+1&lt;&lt;" = ";
                  else
                      cout&lt;&lt;"a["&lt;&lt;i+1&lt;&lt;"]["&lt;&lt;j+1&lt;&lt;"] = ";
                  cin&gt;&gt;a[i][j];
              }
          }
      }
      void display()
      {
          for(int i=0;i&lt;n;i++)
          {
              cout&lt;&lt;"\n";
              for(int j=0;j&lt;n+1;j++)
              {
                  if(j==n)
                      cout&lt;&lt;" ";
                  cout&lt;&lt;a[i][j]&lt;&lt;"\t";
              }
          }
      }
  
      void gauss()//converting augmented matrix to row echelon form
      {
          float temp;//Line 1
          for(int i=0;i&lt;n;i++)
          {
              for(int j=i+1;j&lt;n;j++)
              {
                  temp=a[j][i]/a[i][i];//Line 2
                  for(int k=i;k&lt;n+1;k++)
                  {
                        a[j][k]-=temp*a[i][k];//Line 3
                      //a[j][k]-=a[j][i]*a[i][k]/a[i][i];//Line 4
                  }
              }
          }
      }
  
      void EnterJordan()//converting to reduced row echelon form
      {
          float temp;
          for(int i=n-1;i&gt;=0;i--)
          {
  
              for(int j=i-1;j&gt;=0;j--)
              {
                  temp=a[j][i]/a[i][i];
                  for(int k=n;k&gt;=i;k--)
                  {
                      a[j][k]-=temp*a[i][k];
                  }
              }
          }
  
          float x[n];
          for(int i=0;i&lt;n;i++)//making leading coefficients zero
              x[i]=0;
          for(int i=0;i&lt;n;i++)
          {
              for(int j=0;j&lt;n+1;j++)
              {
                  if(x[i]==0&amp;&amp;j!=n)
                      x[i]=a[i][j];
                  if(x[i]!=0)
                      a[i][j]/=x[i];
              }
          }
      }
      void credits()
      {
          for(int i=0;i&lt;n;i++)
          {
              cout&lt;&lt;"\nx"&lt;&lt;i+1&lt;&lt;" = "&lt;&lt;a[i][n]&lt;&lt;endl;
          }
      }
  
  };
  
  int main()
  {
      Gauss obj;
      obj.accept();
      cout&lt;&lt;"\n\nAugmented matrix: \n\n\n";
      obj.display();
      obj.gauss();
      cout&lt;&lt;"\n\nRow Echelon form: \n\n\n";
      obj.display();
      obj.EnterJordan();
      cout&lt;&lt;"\n\nReduced row echelon form:\n\n\n";
      obj.display();
      cout&lt;&lt;"\n\nSolution: \n\n\n";
      obj.credits();
      return 0;
  }
  </code></pre>
  
  <blockquote>
    <p>Note: My code doesn't take into consideration the problem of division when the pivot is zero (I'm choosing the diagonal element as the pivot every time). 
    For the particular example I tried however, such a case was not encountered.</p>
  </blockquote>
  
  <p>Augmented matrix is :</p>
  
  <pre><code> 2   1  -1    8 
  -3  -1   2   -11    
  -2   1   2   -3
  </code></pre>
  
  <p>The output matrix is :</p>
  
  <pre><code>1   0   0    2  
  0   1   0    3  
  0   0   1    -1 
  </code></pre>
  
  <p>and the solution  is :</p>
  
  <pre><code>x1 = 2
  
  x2 = 3
  
  x3 = -1
  </code></pre>
  
  <p>Using Line 4, the output matrix is :</p>
  
  <pre><code>1   0   0    -0.75  
  0   1   -0   8  
  0   0   1    -1.5   
  </code></pre>
  
  <p>and the solution is :</p>
  
  <pre><code>x1 = -0.75
  
  x2 = 8
  
  x3 = -1.5
  </code></pre>