Member Avatar for pukepuke

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:

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...

  • 3 Contributors
  • 3 Replies
  • 1K Views
  • 7 Years Discussion Span
  • Latest Post Latest Post by Victory_1
Member Avatar for Duoas

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.

Member Avatar for Duoas

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:

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...

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.