%  Reseni soustavy rovnic metodou SOR
%
%  Jan Nejedly FAV2
%
%  10.4.97
%
%
%  Vstup:  A - matice soustavy
%          b - prava strana
%

function x = sor(A, b);

clear x_k;
clear x_k1;
clear U;
clear D;
clear L;
clear cas;
clear H;
clear g;

iterace = 0;
presnost = 1*10^(-6);
dim = size(A);
norma = 1;

tic;

D = diag(diag(A));
U = triu(A) - D;
L = tril(A) - D;

% V nasledujicich radcich je pocitana omega

H=-D\(U+L);
sp=(max(abs(eig(H))))^2;   %spektralni polomer matice

if sp > 1                  %je-li spektralni polomer vetsi nez 1
      w = 1;               %je SOR prevedena na Gauss-Seidel
   else
      w=2/(1+sqrt(1-sp));
end;

% Konec vypoctu omegy

x_k = b;                   % pocatecni aproximace je prava strana

H = (D + w * L) \ ((1 - w) * D - w * U);
g = (D + w * L) \ (w * b);

% Vlastni vypocet
while ((norma >= presnost) & (iterace<5000))
   iterace = iterace +1;
   x_k1 = H * x_k + g;
   norma = 0;
   for i = 1 : dim
     norma = norma + (x_k1(i) - x_k(i))^2;
   end;
   norma = sqrt(norma);

   x_k = x_k1;
end;

cas = toc;

x = x_k1;
disp('Reseni vypoctene programem SOR: ');
disp(x);
disp('Doba vypoctu:                   ');
disp(cas);
disp('Pocet iteraci:                  ');
disp(iterace);
disp('Relaxacni parametr:');
disp(w);

%KONEC
