%
% Risolve il sistema lineare A*x = b
%
function x = linearsistemsolve( A, b )

  [nr,nc] = size(A) ;

  % controllo dimensioni
  if nr ~= nc
    error( sprintf( 'la matrice A (%d x %d) non e` quadrata ', nr,nc) ) ;
  end

  neq = length(b) ;

  if nr ~= neq
    error( sprintf( 'la matrice A (%d x %d) non e` compatibile\ncon la lunghezza del vettore b, length(b) = %d ', nr, nc, neq ) ) ;
  end

  M = [ A b ] ; % costruisco la matrice aumentata
  for i=1:neq-1
    % cerco l'elemento di massimo modulo colonna i-esima
    [V,k] = max(abs(M(i:end,i))) ;

    % rimetto l'indice k alla giusta posizione
    k = k + i-1 ;

    % scambio la riga i con la riga k
    BF     = M(i,:) ;
    M(i,:) = M(k,:) ;
    M(k,:) = BF ;

    % controllo singolaritÃ 
    if  V == 0 
      error( 'Matrice singolare' ) ;  
    end

    % costruisco la matrice di Frobenius
    L             = eye(neq) ;
    L(i+1:end,i)  = -M(i+1:end,i)./M(i,i) ;

    % azzero gli elementi della colonna i
    M             = L * M ;

  end
    
  % alloco il vettore soluzione
  x = zeros(neq,1) ;

  % a questo punto la matrice Ã¨ in forma triangolare
  % applichiamo le operazioni di ritorno
  x(neq) = M(neq,neq) / M(neq,neq+1) ;
  for i=neq-1:-1:1
    x(i) = ( M(i,neq+1) - M(i,i+1:end-1) * x(i+1:end) ) / M(i,i) ;
  end