Gaussian elimination with scaled partial pivoting

[Crzyrio]

Hi

As part of an assigment i am needed to write a C++ Program to solve a system of equations using Gaussian elimination with scaled partial pivoting method.
Now our prof has told us to simple use the pseudocode found in the book. I did my best to finish it however, the answer the program is outputting is not the same as the one i have calculated by hand. Any help would be greatly appreciated

/* 
 * Developer : Rioch D'lyma
 * 
 * Class : 3MN3
 * Due
 * Purpose :
 *
 * Created on January 20, 2011, 12:23 AM
 */
#include <stdlib.h>
#include <cstdlib>
#include <iostream>
#include <math.h>
using namespace std;

/*
 * 
 */
int main(int argc, char** argv) {

    float x[4];
    float sum;
    float temp;
    int i,j,k,n;
    float r,rmax,smax,xmult;
    float a[5][5] = {
                     {0,-1,1,0,-3},
                       {0,1,0,3,1},
                        {0,0,1,-1,-1},
                        {0,3,4,1,2}};
    float b[5] = {4,0,3,1};
    int l[5] = {1,2,3,4};
    float s[5] = {0,0,0,0};
    n = 3;

    for(i=1; i<=n;i++){
        l[i] = i;
        s[i] = 0;
        for(j=1;j<=n;j++){

            if(a[i][j] > s[i]){
                s[i] =a[i][j];
                
            }

        }
        s[i] = smax;
       
    }

 //   cout << s[0] << endl;
   // cout << s[1] << endl;
  //  cout << s[2] << endl;
   // cout << s[3] << endl;
    
    

    for(k=1;k<=n-1;k++){
        rmax = 0;
        for (i=k;i<=n;i++){
           // cout << a[l[i]][k] << " / " << s[l[i]] << endl;
            r = fabs(a[l[i]][k]/s[l[i]]);
            if(r>rmax){
                rmax = r;
                j=i;
               
            }
\
        }
        temp = l[j];
        l[j] = l[k];
        l[k] = temp;


        for(i=k+1;i<=n;i++){
           
           // cout << "K = " << k << endl;
            //cout << "I = " << i << endl;
            //cout << a[l[i]][k] << " / " << a[l[k]][k] << endl;
            xmult = a[l[i]][k]/a[l[k]][k];
            a[l[i]][k] = xmult;
            //cout << xmult << endl;
            for(j=k;j<=n;j++){
                a[l[i]][j] = a[l[i]][j] - (xmult*a[l[k]][k]);
            }
        }

        
    }
    

    for(k=1;k<=n-1;k++){

        for(i=k; i<=n;i++){
            b[l[i]] = b[l[i]] - (a[l[i]][k]*b[l[k]]);
           
        }

    }
    
    x[n] = b[l[n]]/a[l[n]][n];
   

    for(i=n-1;i>=1;i--){
        sum = b[l[i]];
        for (j=i+1;j<=n;j++){
            sum = sum - (a[l[i]][j]*x[j]);
        }
        x[i] = sum/a[l[i]][j];
    }
    for(int r = 0; r <= n;r++){
        cout << "x" << r+1 << " = " << x[r] << endl;
     }

    return 0;
}

Pseudocode

procedure Gauss(n, (ai j ), (li )) integer i, j, k, n;	real r, rmax, smax, xmult real array (ai j )1:n×1:n , (li )1:n ;	real array allocate (si )1:n for i = 1 to n do
li ←i smax ← 0 for j = 1 to n do
smax ← max(smax, |ai j |) end for
si ← smax end for
for k = 1 to n − 1 do rmax ← 0
for i = k to n do r ← |ali ,k /sli |
if (r > rmax) then rmax ← r
j←i end if
end for
lj ↔lk for i = k + 1 to n do
xmult ← ali ,k /alk ,k ali,k ←xmult for j = k + 1 to n do
ali,j ← ali,j −(xmult)alk,j end for
end for end for
deallocate array (si ) end procedure Gauss


procedure Solve(n, (ai j ), (li ), (bi ), (xi )) integer i, k, n;	real sum real array (ai j )1:n×1:n , (li )1:n , (bi )1:n , (xi )1:n for k = 1 to n − 1 do
for i = k + 1 to n do
bli ←bli −ali,kblk end for
end for
xn ←bln/aln,n fori =n−1to1 step−1do
sum ← bli for j = i + 1 to n do
sum ← sum − ali , j x j end for
xi ←sum/ali,i

I have included both my code and the Pseudocode given to us. Array a is the coefficient's and b is the constants.

THanks
Rioch

[daviddoria]

I suggest you come up with a simple example problem, and then do the problem step by step by hand. Then step through the code with a debugger and ensure that it is doing every step the same as you did by hand. This will let you find the bug, then you just need to fix it!

[ravenous]

Then step through the code with a debugger and ensure that it is doing every step the same as you did by hand. This will let you find the bug, then you just need to fix it!

I second that, definitely use a debugger to check the values that are being calculated at each step.

Also, you have done some strange things in your code. You have used #include <stdlib.h> and #include <cstdlib> , which (I think) are the same thing, except that the first is deprecated in C++. You should also use <cmath> instead of math.h .

[Crzyrio]

Thanks for the reply guys!

I have not learnt how to use a debugger but will look it up. However if you notice i have a number of outputs commented out, i have been using that to keep track of what the calculations are doing.

And the problem i have been using has been solved by hand already, as we had to do it by hand and then code.

Ravenous about the

#include

statements. I primarily use netbeans on my mac to write the programs and test it on windows using visual express before handing it in, so i just put both in and clean it up later based on what is required on the windows. In theory that shouldnt happen but netbeans ask i put diff statements than visual.

As for the

math.h

i shall start using

<cmath>

Thanks :)

I ended up figuring it out. My input was in wrong and a few '+1' and "<' were not used correctly :)

[daviddoria]

You should really really learn to use a debugger. Almost every IDE has one built in.

Also, if your problem is solved please mark the thread as 'solved'.