Hey...I'm working on a Gaussian elimination problem...my pivot doesn't seem to work correctly...the logic must be wrong. Also, does anyone know why "INFINITY" comes up as some of my random numbers chosen?

Thanks in advance for any help! :)

Matrix::Matrix():N(0)
{
    //Empty Matrix by default
}

Matrix::Matrix(int n)
{
    N = n;
    A = new double * [N];
    for(int i = 0; i < N; i++)
    {
        A[i] = new double [N];
    }
    s = new double [N];
    b = new double [N];
    for( int i = 0 ; i < N; ++i){
     b[i] = 0.00;
    }
}

void Matrix::generate_A()
{
    for(int i = 0; i < N; i++)
    {
        for(int j = 0; j < N; j++)
        {
            A[i][j] = random();
            cout.width(15);
            cout.precision(2);
            cout.setf(ios::fixed);
            cout << A[i][j];
        }
        cout << endl;
    }
}

void Matrix::generate_s()
{
    for(int i = 0; i < N; i++)
    {
        s[i] = random();
        cout.precision(2);
        cout.setf(ios::fixed);
        cout << s[i] << endl;
    }
}


void Matrix::solve_b(const Matrix& aMatrix, const Matrix& sMatrix )
{
    A = aMatrix.get_A();
    s = sMatrix.get_s();
    assert( A != NULL);
    assert( s != NULL);
    for(int i = 0; i < N; ++i)
    {
        for(int j = 0; j < N; ++j)
        {
            b[i] = b[i] + (A[i][j] * s[i]);
        }
        cout.precision(2);
        cout.setf(ios::fixed);
        cout << b[i] << endl;
    }
}

void Matrix::pivot()
{

    int start, i, k;
    for(int j = 0; j < N; ++j)
    {
        start = i = j;
        k = i + 1;
        while(k < N)
        {
            if (C[i][j] >= C[k][j])
                ++k;
            else
            {
                i = k;
                ++k;
            }
        }
        int p = i;
        if(C[p][j] == 0)
        {
            E = 0;
            return;
        }
        if(p > start)
        {
            double temp;
            for(int m = 0; m < N; ++m)
            {
                temp = C[start][m];
                C[start][m] = C[p][m];
                C[p][m] = temp;
                cout.precision(2);
                cout << C[start][m] << " switched with " << C[p][m] << endl;
            }
        }
    }

    return;
}


void Matrix::solve_x(const Matrix& aMatrix, const Matrix& bMatrix )
{
    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(int row = 0; row < N; row++)
    {
        for(int 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 << "Does pivot work? " << endl;
    pivot();
    cout << "Looks like it " << endl;
 }


int Matrix::get_N()
{
    return(N);
}


double ** Matrix::get_A() const
{
    return (A);
}

double * Matrix::get_s() const
{
    return (s);
}

double * Matrix::get_b() const
{
    return (b);
}

double ** Matrix::get_C()
{
    return (C);
}

double Matrix::random()
{
    double dividend;
    double divisor;
    double rand;
    do
    {
        dividend = ( std::rand() % 1000);
        divisor = ( std::rand() % 100);
        rand = (static_cast <double>(dividend)) / (static_cast<double>(divisor));
    }
    while( dividend < divisor );
    return (rand);
}







int setN()
{
   int n;
    do
    {
        cout << "\nEnter a positive integer n, the size of the (n x n) matrix : ";
        cin >> n;
    }
    while(n < 0);
    return(n);
}

int main()
{
    srand(time(NULL));
    int n = setN();

    //Create matrix A
    Matrix A(n);
    cout << "\nMatrix A:\n";
    A.generate_A();

    //Create sol'n matrix s
    Matrix s(n);
    cout << "\nMatrix s:\n";
    s.generate_s();

    Matrix b(n);
    cout << "\nMatrix b:\n";
    b.solve_b(A, s);

    Matrix x(n);
    cout << "\nMatrix C:\n";
    x.solve_x(A, b);

    system("pause");

    return (0);

}

First, INFINITY comes up because of a divide-by-zero, so check that you don't have a 0 divisor before dividing. NOTE: your random function is really unnecessarily complicated.

Line 59 should read: "b = b + (A[j] * s[j]);" notice the index "j" for the vector s.

Line 94 should read: "for(int m = 0; m <= N; ++m)" notice the <= instead of < because you need to pivot the solution vector as well.

Other then that, I think it's ok. Now you just have the actual elimination to do, which is not harder than what you did already. :)


BTW: just a little tip, to do a pivot, you don't actually need to swap the two rows, you can just add the proper one (the one with max diagonal) to the improper one (the one with smaller diagonal value), but keep it like that, it's fine.

Edited 6 Years Ago by mike_2000_17: n/a

Hi vavazoom

the only thing what immediately leaped out at me is the line where you determine pivot element:

if (C[i][j] >= C[k][j]) ...

// what should be:

if ([B]fabs[/B](C[i][j]) >= [B]fabs[/B](C[k][j])) ...

More serious is where you are exchanging rows:

for(int m = 0; m < N; ++m)
{
temp = C[start][m];
C[start][m] = C[p][m];
C[p][m] = temp;
}

After this loop right side vector elements r[start] of your LEQ Ax=r MUST also be exchanged, for example:

temp = r [ start ]; r [ start ] = r [ p ]; r [ p ] = temp;
where r is right side vector, see my following code.

Because I didn't find this right side vector in your code you may have stored it together with nxn matrix in an nx(n+1) matrix, where column (n+1) contains this vecor.

Below is a method which is taken from a larger matrix class container I wrote.

This method does exactely what are you want to do: determine pivot element, exchanges pivot row (in both nxn matrix and right side vector r), calculates triangular matrix and finally solves Ax=r.

If you are into your code, the kernel algorithm of my method is easily understood. Only consider that in my matrix class array indexes are denoted by (), e.g. conventional C++ a[i,j] is a(i,j) in my class. Then I decided () instead of [] for strictly discriminating between c++ arrays and my own class arrays (and yes, still Fortran adicted).

// Gauss with pivoting to solve Ax=r
template <class T> bool matrix<T>::gsolver (matrix<T>& a, matrix<T>& r, matrix<T>& x)
/*
 Solution of A*x=r by Gauss algorithm with pivoting
 a - coefficient matrix, will be overwritten
 r - right side vector, will be overwritten
 x - resulting vector

 return
 true  - success
 false - singular matrix

 Further reading:
 Carl Friedrich Gauss: Disquisitiones arithmeticae. Hannover(1801).
*/

{ int i, j, k, p, n=a.ncols; T h, e = matrix<T>::getSingularLimit();
  for (i = 0; i<n; i++)
  {
     p = i;
     for (j = i+1; j<n;j++)                                      /* get pivot */
       if (fabs ( a ( p, i ) ) < fabs ( a ( j , i ) )) p = j;
     if (p != i)
     {                                                       /* exchange rows */
      for (k = i; k<n; k++)
      {
        h= a ( i, k ); a ( i, k )= a ( p, k ); a ( p, k )= h;
      }
      h = r ( i ); r ( i ) = r ( p ); r ( p ) = h;
     }
     for (j = i + 1; j<n;j++)
       if (e * fabs ( a ( j, i ) ) < fabs(a ( i, i )))
       { h = a ( j, i ) / a ( i, i ); r ( j ) = r ( j ) - h * r ( i );
         for (k = n-1; k >= i; k--) a ( j, k ) = a ( j, k ) - h * a ( i, k );
       }
        else {a.rank(i+1); return false;}                  /* singular matrix */
  }
  for (j = n-1; j>=0;j--)
   { h = T(0);
     for (i = j + 1; i < n; i++) h = h + a ( j, i ) * x(i);
     h = r ( j ) - h;
     if (e * fabs ( h ) < fabs( a ( j, j ) )) x(j) = h  / a ( j, j );
       else {a.rank(i+1); return false;}                   /* singular matrix */
   }
  a.rank(n); return true;
}

Remarks: n=a.ncols can be set to number of rows/colums. e = matrix<T>::getSingularLimit(): simply set e = 1.0e-300 for singulatity test. drop a.rank() for it only stores current rank of matrix.


-- tesu

Edited 6 Years Ago by tesuji: n/a

@Tesu, note that he does augment the matrix with the RHS vector, see method solve_x() where he allocated C as N x N+1 and fills it up with values properly.

Also @Tesu, try not to post ready-made code. If the OP wanted to just pick a Gaussian elimination algorithm to copy-paste in his own code, he would have done so instead of posting a thread here, it's not hard to find simple algorithms with a simple Google search, but it takes the fun and learning out of doing it yourself. For instance, someone can download the source code for CLARAty and get all matrix numerical methods you can imagine, and much much more, in one shot, or I could give him my own extensive matrix numerical methods libraries.. But that's no fun.

@mike_2000_17, I already assumed that this offside vector is stored in n+1-th column, see my posting. If so, he should replace

---> for(int m = 0; m < N; ++m)

by

---> for(int m = 0; m <= N; ++m)

to swap RHS vector elements also.

As for posting complete code, I believe you are partly right. He did almost already solve the complete Gaussian elimimation algorithm, so my hidden agenda was to show him some errorless code which he might compare with his almost complete solution to figure out some essential flaws he made. I hadn't believe that he would have replaced his fine code by mine.

-- tesu

Edited 6 Years Ago by tesuji: n/a

hahaha...I'm definitely a "she" :)

Thanks for all the advice and help guys...I did make those 2 line changes, however I'm still having the same issue with my pivot function. It's staying in the first column and not moving on to the 2nd column, which is then switching rows over and over wrongly.

Ah, he is definitely a She. So I say sorry for having mistaken he and she. OK, you did make those two line changes, so you did also add missing fabs() twice? However, there is definitely more to change, for example lines 186&204: potentially endless loops.

Can you explain what solve_b and solve_x are supposed to do?

What I am missing completely in your code is where you calculate the LU decomposition (lines 31to 34 of my code) and where you solve the equation system now being in triangular form where you get the solution vector from X(n), X(n-1), ...X(1) by way of back substitution, see lines 38 to 42.


-- tesu

Sorry too for mistaking he for she. Tesuji is right about the fabs() That is probably why rows get switched wrongly.

You should also understand that pivoting is part of the Gaussian elimination. The algorithm is that for every column, you find the row with the maximum (fabs()) value at the column, pivot with the "start" row (whose index is equal to the column index), eliminate the lower off-diagonal entries by adding that row * <some factor>, and repeat until the last row. So the pivot step is in the middle of the main elimination loop. So if you try to pivot for the first column then the next and so on without doing the elimination step, it is not going to do anything useful besides randomly swapping rows around.

@Tesuji:
1. I don't see any endless loops, just properly used "do { } while();" loops.
2. solve_b computes the RHS vector for a solution vector "s" (for checking if the algorithm computes "s" back again)
3. solve_x finds the actual solution vector for Ax = b.
4. the elimination part is not done, she is still at verifying that the pivot step is working correctly (which is a good way to program, kudos to vavazoom). BTW this is a Gaussian elimination, not a LU-decomposition (LU-decomp. is one, slightly better way to solve Ax=b, but this problem is on Gaussian elim.).

I've changed the fabs()...and I actually took out int start because it looked to me as though it was useless. This is my pivot function which I'm basically doing the elimination in as well now...and it's actually turning position 0,0 to a 0.00 and crashing after running through a few times...any help is appreciated...you guys have already helped me out a lot!! :)

void Matrix::pivot()
{

    /*Compute pivot index j <= p <= n
        such that:
        |Cpj| = max(|Cij|) from i = j to n
    */
    int i, k;
    for(int j = 0; j < N; ++j)
    {
        i = j;
        k = i++;
        while(k < N)
        {
            if (fabs(C[i][j]) >= fabs(C[k][j]))
                ++k;
            else
            {
                i = k;
                ++k;
            }
        }

        //If Cpj = 0: set E = 0 and exit
        int p = i;
        if(C[p][j] == 0)
        {
            E = 0;
            cout << "E = "<< E << " <- Should be 0" << endl;
            return;
        }

        //If p > j: interchange rows p & j
        if(p > j)
        {
            double temp;
            for(int m = 0; m <= N; ++m)
            {
                temp = C[j][m];
                C[j][m] = C[p][m];
                C[p][m] = temp;
                cout.precision(2);
                cout << C[j][m] << " switched with " << C[p][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;
}

Thanks!

First, you are incrementing your index j twice in the loop, once at line 9 and once at line 47. That would skip one line. Second, you are starting the "//For each i > j: row i - (Cij/Cjj)* row j" part with "i = j" (because j++ returns the original value of j) not with "i > j" (or i = j + 1) as you should and said in the comment line. So to fix it all:

//replace this:
i = j++;
for(int row_i = i; row_i < N; row_i++)

//with that:
for(int row_i = j+1; row_i < N; row_i++)

Now, your almost there, just a simple back-substitution left to do!

Nice quality code btw, simple and easy to read, not very usual on this forum. If you want to clean up even more (like you did with the start variable), note that "p" is also useless...

This article has been dead for over six months. Start a new discussion instead.