Gauss elimination code doesn't work for sparse matrices

I wrote a matlab code for Gaussian elim to solve for linear system. I tested it on different matrices, but it seems to not be working well at all with matrices which have a lot of zeros in them. Is there any way to fix it? Here is what I have:

 function x = Gauss(A, b)

 [n, n] = size(A);     

 [n, k] = size(b);  

 for i = 1:n-1

    m = -A(i+1:n,i)/A(i,i); 

    A(i+1:n,:) = A(i+1:n,:) + m*A(i,:);

    b(i+1:n,:) = b(i+1:n,:) + m*b(i,:);

 end;

 x(n,:) = b(n,:)/A(n,n);

 for i = n-1:-1:1

    x(i,:) = (b(i,:) - A(i,i+1:n)*x(i+1:n,:))/A(i,i);

 end

Is there a better way to write a function that would solve for matrices with many zeros, for example, say,

 A=[1 2 0 0 0 1; 1 2 2 0 0 0; 0 1 3 2 0 0; 0 0 1 4 2 0;0 0 0 1 5 2; 2 0 0 0 1 6]

 b=[4 5 6 7 8 9]'

Thanks for any input!

[the cyclist]

Can you be more specific about "doesn't work well"? Do you mean it doesn't converge at all? Or seems to be far away from the correct answer?

Answer 1

[John D'Errico]

This is an example of the good, bad, and the ugly. (Ok, I'll concede to serving as the ugly.)

The good? That you wrote this to learn about the method, about when it works, when it will fail, and why it does fail.

The bad? If you intend to use it! The point is, if you think there is no code out there, written by people who truly understand the numerical analysis, the methods, etc., then you are wrong. There are codes written that you can use. I'd look first to tools like LAPACK, or its predecessor, LINPACK. I recall seeing that you can get LAPACK in C. There are repositories of code that will surely contain such a tool.

Don't write basic linear algebra routines yourself. The odds of a bug are too large. And why waste time rewriting tools that are already there to be found, and are free? This is a waste of your time, to get something often of lesser quality.

As for why the code that you wrote failed, it was because of a zero pivot, nearly in the beginning.

A=[1 2 0 0 0 1; 1 2 2 0 0 0; 0 1 3 2 0 0; 0 0 1 4 2 0;0 0 0 1 5 2; 2 0 0 0 1 6]
A =
     1     2     0     0     0     1
     1     2     2     0     0     0
     0     1     3     2     0     0
     0     0     1     4     2     0
     0     0     0     1     5     2
     2     0     0     0     1     6

Look at the first two rows. When you subtract the first row from the second row, look at the result.

A(2,:) = A(2,:) - A(1,:)
A =
     1     2     0     0     0     1
     0     0     2     0     0    -1
     0     1     3     2     0     0
     0     0     1     4     2     0
     0     0     0     1     5     2
     2     0     0     0     1     6

It killed off the element in the (2,1) position, but it also created a zero element at the (2,2) position.

Then you will end up dividing by zero. Oops. A bad thing.

In fact, Gaussian elimination DOES work for sparse matrices, IF you write it properly. In order to fix this, you need to learn about partial pivoting. Better however to just learn about it. DON'T USE THE CODE!

Anyway, for sparse matrices note that there are problems with fill-in of the zero elements in your solution. There are ways to reduce that, by re-ordering the matrix.

And finally, this matrix is NOT a good example of a sparse matrix, as it is not even stored in sparse form, and it is not at all sparse. Nearly half the elements of A are non-zero. That is a waste of time to call it sparse.

Answer 2

[the cyclist]

Is there a reason you are not using a built-in MATLAB function to solve this problem? For example,

x = A\b

See this documentation for details.

[ivordes greenleaf]

I know of this function, I just want to know why the Gaussian elim doesn't work. And how does the mldivide command work anyway? (I guess it might be because I used to use C, and we need to write everything there by hand, so I still have the same habit). Thank you for your answer though.

[the cyclist]

Of course, there is nothing wrong with writing your own code and trying to figure out why it isn't working as expected. I just wanted to make sure you weren't wasting your time, possibly ignorant of an existing solution!

Regarding how mldivide works ... did you look at the documentation link? In particular, there is a flow chart of showing exactly which algorithms are used for different types of inputs.

Unlike many MATLAB functions, you can't necessarily get in with "edit filename.m" and see the actual code, though. But, these are deeply core MATLAB, highly optimized, and presumably pretty bullet-proof at this point (if that is your concern).

Answer 3

[the cyclist]

The presence of the zeros suggests to me that you could have a problem with pivoting.

(Disclaimer: I don't really know anything about this. I just googled enough to be dangerous.)

Answer 4

[Lucas Aloisio]

function resultados=ejercicio1(matriz,ladoderecho)

%Esta funcion recibe como variable de entrada a la matriz de los

%coeficientes y al lado derecho

n=size(matriz,1); %Calculamos al numero de filas que posee el sistema lineal a analizar

ceros=find(matriz==0);

j=1;

if length(ceros)~=0 %#ok<*ISMT>

for i=1:size(matriz,1)%Recorremos a la matriz fila por fila

aux=find(matriz(i,:)==0); %Calculamos al vector auxliar, el cual busca encontrar a la ubicacion geografica de aquellos elementos perntenecientes a la fila i de la matriz y atods sus columnas, que cumplan con el hecho de ser numericamente iguales a 0

if length(aux)~=0 %Si la longitud del vector auxliar es diferente de cero, implica que existen elementos en la fila i de la matriz analizada que cumplen con el hecho de ser numericamente iguales a 0

fila(j)=i;

j=j+1;

c=matriz(i,:);

matriz(i,:)=matriz(n,:);

matriz(n,:)=c;

end

end

end

if length(fila)~=0

for t=1:length(fila)-1

a=length(find(matriz(fila(t),:)==0));

for k=1:length(fila)

b=length(find(matriz(fila(t+1),:)==0));

if a>b

c=matriz(t,:);

matriz(t,:)=matriz(n,:);

matriz(n,:)=c;

end

end

end

end

for k=1:n-1 %Recorremos a la matriz desde la primera fila hasta la ultima

for i=k+1:n

factor = matriz(i,k)/matriz(k,k); %Calculamos al coeficiente con el que multiplicaremos a la fila

for j=k+1:n %Recorremos columna por columna a la matriz de coeficientes

matriz(i,j) = matriz(i,j) - factor*matriz(k,j); %Calculamos a los nuevos coeficientes para todos y cada uno de los elementos que conforman a la fila analizada y almacenamos a dicho resultado en la nuevamente en la matriz de componentes

end

ladoderecho(i)=ladoderecho(i)-factor*ladoderecho(k); %Los elementos del lado derecho tambien se modifican como consecuencia de las operaciones que realizamos sobre la fila analizada

end

end

%Comienzo del proceso de sustitucion inversa

resultados(n)=ladoderecho(length(ladoderecho))/matriz(size(matriz,1),size(matriz,2)); %La componente del vector de resultados es numericamente igual la division entre el ultimo elemento del lado derecho y el ultimo elemento de la matriz de coeficientes

i=n-1;

k=1;

while k==1

if i==0

break

end

suma=sum(ladoderecho(i));

for j=i+1:n

suma=suma-matriz(i,j)*resultados(j);

end

resultados(i)=suma/(matriz(i,i));

i=i-1;

end