Gaussian elimination in pascal

•

What is DaniWeb IT Discussion Community?

You're currently browsing the Pascal and Delphi section within the Software Development category of DaniWeb, a massive community of 455,965 software developers, web developers, Internet marketers, and tech gurus who are all enthusiastic about making contacts, networking, and learning from each other. In fact, there are 3,684 IT professionals currently interacting right now! Registration is free, only takes a minute and lets you enjoy all of the interactive features of the site.

Join our community!

Please support our Pascal and Delphi advertiser: Programming Forums

Thread Tools Display Modes

Views: 1455 | Replies: 2

•

Join Date: Nov 2007

Posts: 3

Reputation: pukepuke is an unknown quantity at this point

Rep Power: 0

Solved Threads: 0

pukepuke

Offline

Newbie Poster

Gaussian elimination in pascal

Nov 17th, 2007

I'm trying to make a program that solvas a linear system Ax =b with a n x n symmetric matriz and a n x 1 vector by gaussian elimination with partial pivotization. The program must search for the pivot element in each column. If the pivot isn't placed in the main diagonal, the program change this line for the line that has the pivot. And it solves the system by gaussian elimination.
I tried this algorithm but it works mainly for 2x2 and 3x3 matrices (it must work for n x n). My program is changing some whole lines in the matrix to zero. So, when it calculates the solution x a division by zero occurs. The higher is n, the higher the chance of getting this error. For some 10 x 10 matrices the program worked too. The algorithm:

Help with Code Tags

(Toggle Plain Text)

Procedure gauss (n: integer;
                 a: matriz;
                 b: vetor);
Var
        x: vetor;
        k,w,e,l :integer;
        c,aux,m,soma :real;
Begin
        w:=1;  {w is a variable for fixed column}
        for k:=1 to n do  {k is the number of iterations}
        begin
                e:=k; {e is the adress of the pivot}
                c:=a[k,w];
                if k<n
                then
                begin
                        for i:=1 to n do {searching for the pivot}
                        begin
                                if k+1 <= n
                                then
                                begin
                                        if abs(a[i+1,w]) > c
                                        then
                                        begin
                                                c:=a[k+i,w];
                                                e:=k+i;
                                        end;
                                end;
                        end;
                        for j:=1 to n do {changing lines}
                        begin
                                aux:=a[k,j];
                                a[k,j]:=a[e,j];
                                a[e,j]:=aux;
                        end;
                        if c <> 0
                        then
                        begin
                                for i:=1 to n do {gaussian elimination}
                                begin
                                        if i > k then
                                        begin
                                                 m:=a[i,w]/a[k,w];
                                                 for j:=1 to n do
                                                 begin
                                                        a[i,j]:=a[i,j]-m*a[k,j];
                                                 end;
                                                 b[i]:=b[i]-m*b[k];
                                        end;
                                end;
                        end;
                end;
                w:=w+1;
        end;
        if a[n,n] <> 0 {calculating x}
        then
        begin
                x[n]:=b[n]/a[n,n];
                for l:=n-1 downto 1 do
                begin
                        soma:=0;
                        for j:=n downto l+1 do
                        begin
                                soma:=soma+x[j]*a[l,j];
                        end;
                        x[l]:=(b[l]-soma)/a[l,l];
                end;
        end
        else
        begin
                writeln('Error');
        end;
end;

Procedure gauss (n: integer;
                 a: matriz;
                 b: vetor);
Var
        x: vetor;
        k,w,e,l :integer;
        c,aux,m,soma :real;
Begin
        w:=1;  {w is a variable for fixed column}
        for k:=1 to n do  {k is the number of iterations}
        begin
                e:=k; {e is the adress of the pivot}
                c:=a[k,w];
                if k<n
                then
                begin
                        for i:=1 to n do {searching for the pivot}
                        begin
                                if k+1 <= n
                                then
                                begin
                                        if abs(a[i+1,w]) > c
                                        then
                                        begin
                                                c:=a[k+i,w];
                                                e:=k+i;
                                        end;
                                end;
                        end;
                        for j:=1 to n do {changing lines}
                        begin
                                aux:=a[k,j];
                                a[k,j]:=a[e,j];
                                a[e,j]:=aux;
                        end;
                        if c <> 0
                        then
                        begin
                                for i:=1 to n do {gaussian elimination}
                                begin
                                        if i > k then
                                        begin
                                                 m:=a[i,w]/a[k,w];
                                                 for j:=1 to n do
                                                 begin
                                                        a[i,j]:=a[i,j]-m*a[k,j];
                                                 end;
                                                 b[i]:=b[i]-m*b[k];
                                        end;
                                end;
                        end;
                end;
                w:=w+1;
        end;
        if a[n,n] <> 0 {calculating x}
        then
        begin
                x[n]:=b[n]/a[n,n];
                for l:=n-1 downto 1 do
                begin
                        soma:=0;
                        for j:=n downto l+1 do
                        begin
                                soma:=soma+x[j]*a[l,j];
                        end;
                        x[l]:=(b[l]-soma)/a[l,l];
                end;
        end
        else
        begin
                writeln('Error');
        end;
end;

If someone could help me with this...

•

Join Date: Oct 2007

Location: Cherry Hill, NJ

Posts: 1,878

Reputation: Duoas is a splendid one to behold

Rep Power: 13

Solved Threads: 193

Duoas

Offline

Posting Virtuoso

Re: Gaussian elimination in pascal

Nov 20th, 2007

Hey there, just a heads-up to let you know that help is forthcoming. I just had to go learn how to solve a linear system of equations using gaussian elimination and partial pivotization first.

Now I'm looking over your code. I've got a few suggestions about commentary and variable names, but that will come when I post again.

•

Join Date: Oct 2007

Location: Cherry Hill, NJ

Posts: 1,878

Reputation: Duoas is a splendid one to behold

Rep Power: 13

Solved Threads: 193

Duoas

Offline

Posting Virtuoso

Re: Gaussian elimination in pascal

Nov 21st, 2007

Eh, OK. Don't take this the wrong way, but you need to read some more on gaussian elimination using partial pivotization. Your code is missing loops where there should be loops, and where there are loops, you are indexing stuff you shouldn't.

BTW. I'm sure you know that the code you posted doesn't compile. I haven't tried, but right at the start I noticed two fatal errors.

First you assign c (an integer) an element of A (a real). Without explicit conversion in Pascal this will fail as a type error. Even so, the algorithm will not work over integers.

Second, you are using a variable named "i" which you never declared. In fact, I think you are being confused by the choices you have made in naming your variables.

Here is a rule that will help greatly:

•

Never use one letter variable names unless:

1. The name of the variable is "not important" AND you are using it in about ten lines of code or less.

By "not important" I mean that the meaning of variable itself is ancillary. For example:
Help with Code Tags

(Toggle Plain Text)
procedure print_string( str: pchar; length: integer );
  var i: integer;
  begin
  for i := 0 to length-1 do write( str[ i ] )
  end;
procedure print_string( str: pchar; length: integer ); var i: integer; begin for i := 0 to length-1 do write( str[ i ] ) end;
While this is an obviously contrived example, it illustrates the point: the name of the variable "i" is unimportant because it is immediately apparent what the code that uses it does.

However, once you get into multiple loops, or multidimensional arrays, or anything of any complexity at all, find a name for the variable that expresses its purpose without having to figure out the code around it. So, if I were doing something complex to each element of str, I would name my index variable something like "index_into_str".

2. The name is a canonical name for something well-known and it is used obviously.

For example, you are generally safe using "x" and "y" index a pixel in an image, because it is such a common idiom that the meaning is apparent.

This rule applies equally to all no-thought variable names. So "zz" is not excluded from the rule.

The key concepts here are transparency (making the meaning of things obvious), and thought (you cannot avoid thinking about how things work and relate).

So then, you could get away with naming a row-counter "r", but you are better off just naming it "row", because when you scroll down 15 lines you don't have to remember what "r" is or does: it says "row" --then it must be the row number, so everywhere you see "row" you know you are doing something that has to do with the current row.

You have not given information about your matrix or vector types other than that the matrix is square. I presume this to mean that the matrix sits on the left-hand side of your equations and the vector on the right-hand side. There is a pretty good psuedo-code available on the Wikipedia to do exactly what you want. However, it assumes an m x n matrix, where n = m+1. I don't know why you have the rhs in a separate vector, but you can still do that if you like; you'll just have to check where you are trying to index in the algorithm and choose the matrix or vector appropriately.

Finally, while you do have some comments, your code is not clear enough to justify their sparsity. Make sure you explain clearly what it is you are trying to accomplish in the next few lines of code, and when you (or someone else) look it over you will be able to immediately see whether or not the code does what it says it is supposed to do.

I know this isn't a whole lot of help (at least, not like you expected). Is this a school assignment? or are you doing this for work? or just on your own? I can help you refine your code, but you need to think its pieces through a bit more...