Метод Gauss-Seidel не работает для больших разреженных массивов?

Question

Еще раз у меня проблема с методом Гаусса-Зейделя в Matlab. Вот он:

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 

Это действительно работает отлично для матрицы 100x100 или меньше. Однако мой наставник хочет, чтобы он работал на матрицах 100000x100000. Сначала это было трудно даже создать саму матрицу, но мне удалось сделать это с небольшой помощью здесь: Matlab Help Center

Теперь я вызываю функцию ex1_3 с А в качестве параметра, но он идет очень медленно. На самом деле это никогда не заканчивается. Как я могу заставить его работать?

Вот мой код для создания конкретной матрицы мой учитель хотел: Важной частью является только, что она отвечает этим условиям: А (я, я) = 3, А (я - 1; я) = A (I ; + 1) = -1 п = 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); 

Answer 1

для цикла работает очень медленно в Matlab, возможно, вы можете попробовать matrix form итерации:

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 

для вашего A и b, он принимает только 16 ste ps для сходимости.

tril извлекает нижнюю треугольную часть A, вы также можете получить это Q при создании матрицы. Поскольку Q уже является треугольной матрицей, вы можете solve уравнение Q*dx=r очень легко, если вам не разрешено использовать функцию \.