Gaussian Elimination Linear Equation in C++

I have a fairly basic maths question but the trick is I need it in C++. I'm following the pseudocode given on Wikipedia at the moment. Here is my attempt:

createMatrixForAllSolutions(*this);
std::cout << equationMatrix.to_string() << endl;
bool solved = false;
int rows = equationMatrix.getRows();
int cols = equationMatrix.getCols();
int i = 0;
int j = 0;
int maxi = 0;
double current = 0;
double eqnValue = 0;
double solValue = 0;
std::vector<char> reversedVars;
int sum = 0;
int tempValue;
int tempRHS;
int newValue;
int neRHS;

while (i < rows && j < cols) {
    maxi = i;
    for (int k = i + 1; k < rows; k++) {
        if (abs(equationMatrix.get_element(k, j)) > abs(equationMatrix.get_element(maxi, j)))
            maxi = k;
    }
    if (equationMatrix.get_element(maxi, j) != 0) {
        current = equationMatrix.get_element(i, j);
        for (int x = 0; x < cols; x++) {
            tempValue = equationMatrix.get_element(i, x);
            newValue = equationMatrix.get_element(maxi, x);
            equationMatrix.set_element(i, x, newValue/current);
            equationMatrix.set_element(maxi, x, tempValue);
        }
        tempRHS = solutionMatrix.get_element(i, 0);
        neRHS = solutionMatrix.get_element(maxi, 0);
        solutionMatrix.set_element(i, 0, neRHS/current);
        solutionMatrix.set_element(maxi, 0, tempRHS);
        //SWAP rows i and maxi
        //SWAP RHS i and maxi
        //DIVIDE each entry in row i by开发者_如何学运维 current
        //DIVIDE RHS i by current
        for (int u = i + 1; u < rows; u++) {
            eqnValue = equationMatrix.get_element(u, j) - equationMatrix.get_element(i, j) * equationMatrix.get_element(u, j);
            std::cout << "Equation Value: " << eqnValue << endl;
            equationMatrix.set_element(u, j, eqnValue);
            solValue = solutionMatrix.get_element(u, 0) - solutionMatrix.get_element(i, 0) * solutionMatrix.get_element(u, 0);
            std::cout << "Solution Value: " << solValue << endl;
            solutionMatrix.set_element(u, 0, solValue);
        }
        i++;
    }
    j++;
}

And the pseudo code I am following is from Wikipedia:

i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:
  maxi := i
  for k := i+1 to m do
    if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi := k
    end if
  end for
  if A[maxi,j] ≠ 0 then
    swap rows i and maxi, but do not change the value of i
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry in row i by A[i,j]
    Now A[i,j] will have the value 1.
    for u := i+1 to m do
      subtract A[u,j] * row i from row u
      Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end for
    i := i + 1
  end if
  j := j + 1
 end while

I've done my best match it so far but if anyone is able to figure out why my doe isn't working, it would be lovely. Thanks!

Here's one problem: you are declaring tempValue, tempRHS, newValue, and neRHS all as ints. Even if your matrix starts off with all integer values, they won't stay that way long once you get into the elimination. These should all be declared double - as ints, you will be constantly throwing away the fractional parts.

I didn't go through your code line by line, but the most likely problem lies with the the difference in array indexing conventions between your C++ implementation and the Wikipedia algorithm. C++ uses 0-based arrays (i.e. the index of the first array element is 0) while the Wikipedia algorithm is 1-based (i.e. the index of the first array element is 1). It's likely you missed some conversion somewhere.

Assuming you understand what the algorithm is trying to do, your best bet will be to scrap your code and start over based on your understanding of the algorithm and C++.

If you are having some difficulty with understanding the algorithm, you might want to look at a copy of Numerical Recipes in C (Chapter 2 will be most useful in this regard). Since C and C++ are both 0-based array languages, the changes you need should be relatively minor compared with using the Wikipedia version as a basis for your implementation

You are not dividing by the correct pivot element. At each step of the algorithm the pivot element is A(maxi,j). However, in your code you try to merge the two steps of swapping and dividing by the pivot element. As a result you say

  current = equationMatrix.get_element(i,j) 

This should read

  current = equationMatrix.get_element(maxi,j) 

I hope that helps. There may be more bugs. This is just the first one I saw. It might be helpful when debugging to print out the current matrix at each step of the algorithm. When your Gaussian elimination is working correctly your algorithm will produce an upper triangular matrix with 1's on the diagonal. (For more details see the Wikipedia page.)

I hope this code is for educational purposes. There are many fine linear algebra libraries out there (such as LAPACK). If you need to solve linear systems I highly recommend using one of these nice libraries rather than trying to roll your own.