CODEWITHSUNDEEP
matlabmatrixlinear-algebranumerical-methods
Sofia
Sofia

Reputation: 456

Gauss-Seidel method doesn't work for large sparse arrays?

Once again I have a problem with the Gauss-Seidel Method in Matlab. Here it is: 

function [x] = ex1_3(A,b)

format long 

sizeA=size(A,1);

x=zeros(sizeA,1);  

%Just a check for the conditions of the Gauss-Seidel Method (if it has dominant diagonal)

for i=1:sizeA
    sum=0;
    for j=1:sizeA
        if i~=j 
            sum=sum+abs(A(i,j));
        end
    end
    if abs(A(i,i))<sum
        fprintf('\nGauss-Seidel''s conditions not met!\n');
        return     
    end
end

%Actual Gauss-Seidel Method

max_temp=10^(-6); %Pass first iteration
while max_temp>(0.5*10^(-6))
    xprevious=x;
    for i=1:sizeA
        x(i,1)=b(i,1);
        for j=1:sizeA
            if i~=j
              x(i,1)=x(i,1)-A(i,j)*x(j,1); 
            end
        end
        x(i,1)=x(i,1)/A(i,i);
    end
    x
    %Calculating infinite norm of vector x-xprevious 

    temp=x-xprevious;
    max_temp=temp(1,1);
    for i=2:sizeA
       if abs(temp(i,1))>max_temp
           max_temp=abs(temp(i,1));
       end
    end
end

It actually works fine for a 100x100 matrix or smaller. However, my tutor wants it to work for 100000x100000 matrices. At first it was difficult to even create the matrix itself, but I managed to do it with a little help from here: Matlab Help Center

Now, I call the ex1_3 function with A as a parameter, but it goes really slow. Actually it never ends. How can I make it work?

Here's my code for creating the specific matrix my tutor wanted: The important part is just that it meets these conditions: A(i; i) = 3, A(i - 1; i) = A(i; i + 1) = -1 n=100000

b=ones(100000,1);
b(1,1)=2;
b(100000,1)=2;

i=zeros(299998,1); %Matrix with the lines that we want to put nonzero elements 
j=zeros(299998,1);  %Matrix with the columns that we want to put nonzero elements 
s=zeros(299998,1); %Matrix with the nonzero elements. 
number=1; 
previousNumberJ=0;
numberJ=0;
for k=1:299998 %Our index in i and j matrices
    if mod((k-1),3)==0
        s(k,1)=3;
    else
        s(k,1)=-1;
    end
    if k==1 || k==2
        i(k,1)=1;
        j(k,1)=k;
    elseif k==299997 || k==299998   
        i(k,1)=100000;
        j(k,1)=(k-200000)+2;
    else
        if mod(k,3)==0
            number=number+1;
            numberJ=previousNumberJ+1;
            previousNumberJ=numberJ;
        end
        i(k,1)=number;
        j(k,1)=numberJ;
        numberJ=numberJ+1;
    end
end

A=sparse(i,j,s); %Creating the sparse array

x=ex1_3(A,b);

Upvotes: 0

Views: 1832

Answers (1)

lennon310
lennon310

Reputation: 12689

the for loop works very slowly in Matlab, perhaps you may want to try the matrix form of the iteration:

function x=gseidel(A,b)
    max_temp=10^(-6); %Pass first iteration
    x=b;
    Q=tril(A);
    r=b-A*x;

    for i=1:100
        dx=Q\r; 
        x=x+1*dx; 
        r=b-A*x; 

        % convergence check
       if all(abs(r)<max_temp) && all(abs(dx)<max_temp), return; end
    end

For your A and b, it only takes 16 steps to converge.

tril extracts the lower triangular part of A, you can also obtain this Q when you build up the matrix. Since Q is already the triangular matrix, you can solve the equation Q*dx=r very easily if you are not allowed to use \ function.

Upvotes: 1

Related Questions