Elimination part of Gaussian Elimination

[vavazoom]

I'm having an issue with the elimination portion of my Gaussian Elimination problem. The program crashes while running, I can't figure out what is causing the wrong rows to switch, and why it is not eliminating correctly.

void Matrix::gaussianElimination()
{

    int i, j, max;
    for(j = 0; j < N; ++j)
    {
        i = j;
        max = i++;
        
        for(max = i++; max < N; max++)
        while(max < N)
        {
            if (fabs(C[i][j]) >= fabs(C[max][j]))
                ++max;
            else
            {
                i = max++;
            }
        }

        if(C[max][j] == 0)
        {
            E = 0;
            return;
        }

        //If p > j: interchange rows p & j
        if(max > j)
        {
            double temp;
            for(int m = 0; m <= N; ++m)
            {
                temp = C[j][m];
                C[j][m] = C[max][m];
                C[max][m] = temp;
                cout.precision(2);
                cout << C[j][m] << " switched with " << C[max][m] << endl;
            }
        }
        //For each i > j: row i - (Cij/Cjj)* row j
        i = j++;
        for(int row_i = i; row_i < N; row_i++)
        {
            for(int row_j = j; row_j <= N; row_j++)
            {
                C[row_i][row_j] = C[row_i][row_j] - ((C[i][j]/C[j][j]) * C[j][row_j]);
                cout.precision(2);
                cout.width(10);
                cout << C[row_i][row_j];
            }
            cout << endl;
        }
    }

    return;
}

        //Elimination to create a diagonal
        for(i = N; i >= j; --i)
        {
            for(k = j++; k < N; ++k)
            {
                C[k][i] -= C[k][j]/C[j][j] * C[j][i];
            }
            for(k = 0; k < N; ++k)
            {
                for(i = 0; i <= N; ++i)
                {
                    cout.width(10);
                    cout.precision(2);
                    cout.setf(ios::fixed);
                    cout << C[k][i];
                }
                cout << endl;
            }
        }
    }


    for(i = N--; i >=0; --i)
    {
        C[i][N] = C[i][N] / C[i][i];
        C[i][i] = 1;
        for(j = i--; j >= 0; --j)
        {
            C[j][N] -= C[j][i] * C[i][N];
            C[j][i] = 0;
        }
    }

    return;
}


void Matrix::solve_x(const Matrix& aMatrix, const Matrix& bMatrix )
{
    int row, col;
    A = aMatrix.get_A();
    b = bMatrix.get_b();
    assert( A != NULL);
    assert( b != NULL);
    E = 1;

    C = new double * [N];
    for(int i = 0; i < N; i++)
    {
        C[i] = new double[N + 1];
    }

    for(row = 0; row < N; row++)
    {
        for(col = 0; col <= N; col++)
        {
            if(col < N)
            {
                C[row][col] = A[row][col];
            }
            else if(col == N)
            {
                C[row][col] = b[row];
            }
            cout.width(10);
            cout.precision(2);
            cout.setf(ios::fixed);
            cout << C[row][col];

        }
        cout << endl;
    }

    cout << "Start Gaussian Elimination" << endl;
    gaussianElimination();

    for(row = 0; row < N; ++row)
    {
        for(col = 0; col <= N; ++col)
        {
            cout.width(10);
            cout.precision(2);
            cout.setf(ios::fixed);
            cout << C[row][col];
        }
    }
 }

Any suggestions, tips, or ideas are appreciated! Thanks...

[arkoenig]

This is more code than I want to bother reading without a clue as to what I'm looking for. However, the following oddity did jump out at me:

The for statement on line 10 has as its subject the while statement on lines 11 through 19. Is that what you had in mind? If so, why is the code not indented that way? If not, what did you have in mind?

And while we're at it, I should note that the code is a little weird: Two nested loops, each of which has max < N as its condition, and each of which increments max each time through. Which means that when the inner loop terminates, so does the outer one--in which case I don't understand why you're bothering to write a loop at all.

So I am wondering whether perhaps line 10 was really intended to be there, or is left over from an earlier version of the program by mistake.

[vavazoom]

Yes, I had a few different versions I was playing with...the for() was supposed to be taken out.

[tesuji]

Hi, glad to see you again

I thought her old but similar problem would have already been solved by master mike (Indeed, I had mistaken decomposition with elimination, I say sorry for that.)

-- tesu

[mike_2000_17]

Well "master mike" is back to quote tesuji...

Your program crashes because of this portion of the code (line 40 to 52):

//For each i > j: row i - (Cij/Cjj)* row j
        i = j++;
        for(int row_i = i; row_i < N; row_i++)
        {
            for(int row_j = j; row_j <= N; row_j++)
            {
                C[row_i][row_j] = C[row_i][row_j] - ((C[i][j]/C[j][j]) * C[j][row_j]);
                cout.precision(2);
                cout.width(10);
                cout << C[row_i][row_j];
            }
            cout << endl;
        }

The problem is that you say, correctly, for each i > j, but the first line "i = j++;" will return the initial value of j and increment j. So you are starting from i = j" not "i = j + 1", this causes the first iteration to eliminate the row that is used to cancel the other rows (row j) and thus, for all other row elimination after that, you get a divide-by-zero and it crashes on an exception, I presume. (I think I actually pointed that out in your previous thread too, this is not a mistake that would go unseen by my C++ eagle eye) So replace "i = j++;" by "i = j + 1;" and it should work (not ++j because you already increment j in the outermost for-statement).