Found some of my earlier college C++ work on an old hard drive... man this was bad.

gistfile1.cpp

template <class T>
Matrix<T> gaussian_solver<T>::operator()(const Matrix<T>& A, const Vector<T>& b) {
   
   if(!A.is_square()) {
      throw CSignal("Gaussian solver can only be applied to square matrices");
   }
   
   Matrix<T> system = A;       // make a copy because we're going to mess it up
   Vector<T> row_temp;
   system = system.augment(b); // Augment vector b, now we have a system
   
   const size_t m = system.getRows();
   const size_t n = system.getCols();
 
   size_t i=0, j=0, max_ind=0, temp=0;
   T val, max_val;
   
   while (i < m && j < n)
   {
      //Find pivot in column j, starting in row i
      max_val = system(i,j);
      max_ind = i;
  
      for(size_t k=i; k<m; k++)
      {
         val = system(k, j);
         if(tabs(val) > tabs(max_val))
         {
            max_val = val;
            max_ind = k;
         }
      }
      if(tabs(max_val) > epsilon)
      {
         system.switchRow(i,max_ind);
         
         // Now system(i, j) will contain the same value as max_val
         
         // divide row i by max_val
         for(size_t k=0; k < n; k++)
         {
            system(i,k) = system(i,k) / max_val;
         }
         // Now sytem(i, j) will have the value 1
         
         for(size_t u=0; u < m; u++)
         {
            if (u != i)
            {
               // subtract A[u, j] * row i from row u
               for(size_t k=i+1; k<n; k++)
               {  
                  system(u,k) -= (system(u,j) * system(i,k));
               }
            }
         }
         i++;
      }
      j++;
   }
   return system;
}

template <class T>
void gaussian_solver<T>::setprecision(double eps)
{
   if(eps < 1 && eps > 0)
      epsilon = eps;
}

template <class T>
T gaussian_solver<T>::tabs(const T number)
{
   T temp = number;
   
   if(temp < 0) {
      temp = -temp;
   }
   return temp;
}