Gauss elimination code doesn't work for sparse matrices
7 views (last 30 days)
Show older comments
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!
1 Comment
the cyclist
on 24 Apr 2015
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?
Answers (4)
John D'Errico
on 24 Apr 2015
Edited: John D'Errico
on 24 Apr 2015
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.
0 Comments
the cyclist
on 24 Apr 2015
Is there a reason you are not using a built-in MATLAB function to solve this problem? For example,
x = A\b
2 Comments
the cyclist
on 24 Apr 2015
Edited: the cyclist
on 24 Apr 2015
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).
the cyclist
on 24 Apr 2015
(Disclaimer: I don't really know anything about this. I just googled enough to be dangerous.)
0 Comments
Lucas Aloisio
on 15 Sep 2021
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
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!