function Z = trisup_inv(U)
%TRISUP_INV Inverte una matrice triangolare superiore
%
% Z = TRISUP_INV(U)
%
% Dati di ingresso:
%   U:    matrice quadrata triangolare superiore
%
% Dati di uscita:
%   Z:    matrice triangolare superiore inversa della U

% Controlli
[m,n] = size(U);
if (m~=n)
    error('La matrice U non e'' quadrata')
end

if ~isequal(U,triu(U))
    error('La matrice U non e'' triangolare superiore')
end

if prod(diag(U)) == 0
    error('La matrice U e'' singolare e quindi non invertibile')
end

% Inizializzo la diagonale in modo vettorizzato
% per non inserirla nel ciclo e per preallocarla
Z = diag(1./diag(U));

for i = n:-1:1
    for j = i-1:-1:1
        somma = 0;
        for k =j+1:i
            somma = somma +U(j,k)*Z(k,i);
        end
        Z(j,i) = -1/U(j,j)*somma;
    end
end