gaussian.sci

function [singular, x] = gaussian(A, b)
  rozmiar = size(A, 1);
  zmiany = zeros(rozmiar, 1);
  e = 1E-12;
  eglowny = 0;
  bool = 1;
  
  for i=1:rozmiar
      zmiany(i,1) = i;
     //wybor elementu glownego
      eglowny = abs(A(i, i));
      xx = i;
      y = i;
      for j = i : rozmiar
        for k = i : rozmiar
            if (eglowny < abs(A(j, k))) then
                eglowny = abs(A(j, k));
                y = j;
                xx = k;
            end;
        end;
      end;
      
      
      //warunek e
      if (eglowny <= e) then
          bool=0;
          break;
      end;
      
      //Zamiana wierszy
      [Aa, bb] = rotacjaW(A, b, y, i);
	  A = Aa;
	  b = bb;
      
      //Zamiana kol
      [Aaa, zm] = rotacjaK(A, zmiany, xx, i);
    
	  A = Aaa;
	  zmiany = zm;
	  
	  //eliminacja
	  for j=(i+1) : rozmiar
          b(j, :) = b(j, :) - b(i,:)*A(j, i)/A(i, i);
          A(j, :) = A(j, :) - A(i,:)*A(j, i)/A(i, i);
      end    
	  
end;
	
//det
  dett=1;
  for i=1 : rozmiar
    dett=dett*A(i,i);
  end;
 
 //sprawdzenie
 if (bool==0 | dett==0) then 
   singular=1;
   x=0;
else	   
	   
 
 //obliczanie x
  for i=1 : rozmiar
      m = rozmiar + 1 - i;
      b(m, 1) = b(m, 1)*(A(m, m)^(-1));
      A(m, :) = A(m, :)*(A(m, m)^(-1));
      for j=i : (rozmiar-1)
          y = rozmiar-j;
          b(y, 1) = b(y, 1) - A(y, m)*b(m, 1);
          A(y, m) = 0;
      end;
  end;
 
 //uwzglednienie rotacji kol w wyniku 
    for i=1 : rozmiar
        if ((rozmiar-i+1) <> zmiany((rozmiar-i+1), 1)) then 
          wiersz = b((rozmiar-i+1), 1);
          b((rozmiar-i+1), 1) = b(zmiany((rozmiar-i+1), 1), 1);
          b(zmiany((rozmiar-i+1), 1), 1) = wiersz;
        end;
    end;
	
  singular = 0;
  x = b;
  end;
endfunction

function [Aa,bb] = rotacjaW(A, b, y, i)
  if (y <> i) then
          wiersz = A(i, :);
          A(i, :) = A(y, :);
          A(y, :) = wiersz;
          wiersz = b(i, :);
          b(i, :) = b(y, :);
          b(y, :) = wiersz;
      end;
  
  Aa = A;
  bb = b;
endfunction

function [Aaa, zm] = rotacjaK(A, zmiany, x, i)
  if (x <> i) then
          kol = A(:, i);
          A(:, i) = A(:, x);
          A(:, x) = kol;
          zmiany(i,1) = x;
      end;
	  
	  Aaa = A;
	  zm = zmiany;
endfunction