title: Gauss Elimination for NxM matrix
body: <pre><code>/* Program to demonstrate gaussian &lt;strong class="highlight"&gt;elimination&lt;/strong&gt;
     on a set of linear simultaneous equations
   */
  
  #include &lt;iostream&gt;
  #include &lt;cmath&gt;
  #include &lt;vector&gt;
  
  using namespace std;
  
  const double eps = 1.e-15;
  
  /*Preliminary pivoting strategy
    Pivoting function
   */
  double pivot(vector&lt;vector&lt;double&gt; &gt; &amp;a, vector&lt;double&gt; &amp;b, int i)
   {
       int n = a.size();
       int j=i;
       double t=0;
  
       for(int k=i; k&lt;n; k+=1)
       {
           double aki = fabs(a[k][i]);
           if(aki&gt;t)
           {
               t=aki;
               j=k;
           }
       }
       if(j&gt;i)
       {
           double dummy;
           for(int L=0; L&lt;n; L+=1)
           {
               dummy = a[i][L];
               a[i][L]= a[j][L];
               a[j][L]= dummy;
           }
           double temp = b[j];
           b[i]=b[j];
           b[j]=temp;
       }
       return a[i][i];
   }        
  
  
  
  /* Forward &lt;strong class="highlight"&gt;elimination&lt;/strong&gt; */
  void triang(vector&lt;vector&lt;double&gt; &gt; &amp;a, vector&lt;double&gt; &amp;b) 
  {
      int n = a.size();
      for(int i=0; i&lt;n-1; i+=1)
      {
          double diag = pivot(a,b,i);
          if(fabs(diag)&lt;eps)
          {
              cout&lt;&lt;"zero det"&lt;&lt;endl;
              return;
          }
          for(int j=i+1; j&lt;n; j+=1)
          {
              double mult = a[j][i]/diag;
              for(int k = i+1; k&lt;n; k+=1)
              {
                  a[j][k]-=mult*a[i][k];
              }
              b[j]-=mult*b[i];
          }
      }
  }
  
  /*
  DOT PRODUCT OF TWO VECTORS
  */
  double dotProd(vector&lt;double&gt; &amp;u, vector&lt;double&gt; &amp;v, int k1,int k2)
  {
    double sum = 0;
    for(int i = k1; i &lt;= k2; i += 1)
    {
       sum += u[i] * v[i];
    }
    return sum;
  }
  
  /*
  BACK SUBSTITUTION STEP
  */
  void backSubst(vector&lt;vector&lt;double&gt; &gt; &amp;a, vector&lt;double&gt; &amp;b, vector&lt;double&gt; &amp;x)
  {
      int n = a.size();
    for(int i =  n-1; i &gt;= 0; i -= 1)
    {
      x[i] = (b[i] - dotProd(a[i], x, i + 1,  n-1))/ a[i][i];
    }
  }
  /*
  REFINED GAUSSIAN &lt;strong class="highlight"&gt;ELIMINATION&lt;/strong&gt; PROCEDURE
  */
  void gauss(vector&lt;vector&lt;double&gt; &gt; &amp;a, vector&lt;double&gt; &amp;b, vector&lt;double&gt; &amp;x)
  {
     triang(a, b);
     backSubst(a, b, x);
  }                               
  
  // EXAMPLE MAIN PROGRAM
  int main()
  {
      int n;
      cin &gt;&gt; n;
      vector&lt;vector&lt;double&gt; &gt; a;
      vector&lt;double&gt; x;
      vector&lt;double&gt; b;
      for (int i = 0; i &lt; n; i++) {
          vector&lt;double&gt; temp;
          for (int j = 0; j &lt; n; j++) {
              int no;
              cin &gt;&gt; no;
              temp.push_back(no);
          }
          a.push_back(temp);
          b.push_back(0);
          x.push_back(0);
      }
      /* 
      for (int i = 0; i &lt; n; i++) {
          int no;
          cin &gt;&gt; no;
          b.push_back(no);
          x.push_back(0);
      }
      */
  
    gauss(a, b, x);
    for (size_t i = 0; i &lt; x.size(); i++) {
        cout &lt;&lt; x[i] &lt;&lt; endl;
    }
     return 0;
  }
  </code></pre>
  
  <p>The above gaussian eleimination algorithm works fine on NxN matrices. But I need it to work on NxM matrix. Can anyone help me to do it? I am not very good at maths. I got this code on some website and i am stuck at it.</p>