My Gauss-Seidel function takes forever to run

Below is my code

function [daan]=GaussSeidel(A,b,n,iter,th)

% A = Left Matrix
% b = Right Matrix
% n = Number of Grid Points
% iter = Number of Iterations
% th = Threshold Value for Stopping

x=zeros(n,iter+1); % Initialize x (equivalent to guess xj to be zero)
daan=NaN; 
for k=1:iter
    dbpn=0; % Different Between Previous and New
    for i=1:n
        s1=0;
        s2=0;
        for j=1:i-1
            s1=s1+A(i,j)*x(j,k+1); % First Sum in Equation
        end
        for j=i+1:n
            s2=s2+A(i,j)*x(j,k); % Second Sum in Equation
        end
        x(i,k+1)=(b(i)-s1-s2)/A(i,i);
        dbpn=dbpn+abs(x(i,k+1)-x(i,k)); 
    end
    if abs(dbpn)<th
        daan=x(:,k); % Determine When is Good to Stop Iterating
        break
    end
end

if isnan(daan)
    error('Pls Increase Threshold or Times of Iteration')
else
    disp(strcat('Number of Iteration =',{' '},num2str(k-1)))
    disp(strcat('Error =',{' '},num2str(dbpn)))
end

end

My TA told me that for this specific problem, it should converge after about 10 iterations, but for my code, it takes me like forever to get the results. A is typically a 10000*10000 matrix, with most elements being 0. b is 10000*1. When I reduce the size of A, it gives me similar results to A\b.

Please help me, thanks!

[John D'Errico]

But without showing what problem you tried it on, it is difficult to know if your TA knows what they are talking about. Attach the matrix as a .mat file to a comment. I'd have been surprised if Gauss-Seidel did converge in only 10 iterations.

The only comment I would make is that Gauss-Seidel can be written a bit more simply, without the use of explicit loops. That would make each iteration run far more rapidly.

Even better is if your matrix (with most elements zero) is stored in sparse format.

[Mingchen Mao]

Thanks, I will try to get rid of those j loops.

[John D'Errico]

If you look in almost any G-S reference, you will find a matrix form for G-S iteration.

https://en.wikipedia.org/wiki/Gauss–Seidel_method

It should be of the form

x_kplus1 = inv(L)*(b - U*x_k)

Here U is the purely upper triangle, above the diagonal. triu will give it to you. L is the lower triangle (including the diagonal.)

L = tril(A);
U = triu(A,1);

So you can write G-S iteration using the above equation, but I doubt you really want to use inv there. Yes, you can write it using backslash too, as

x_kplus1 = L\(b - U*x_k)

Since L is purely lower triangular, you can do that operation as a forward substitution, so a simple loop, which is itself vectorizable. That avoids the use of any implicit matrix solution.

L and U are fixed. So any iteration will look vaguely like this:

temp = b - U*x_k;
x_kplus1 = zeros(n,1);
x_kplus1(1) = temp(1)/L(1,1);
for ii = 2:n
  x_kplus1(ii) = (temp(ii) - L(ii,1:ii-1)*x_kplus1(1:ii-1))/L(ii,ii);
end

Note quite the one-liner you see above, but still reasonably well vectorized.