%====================================================================
% Script per determinare l'inversa di una matrice triangolare sup.
% Necessita della function trisup_inv 
%====================================================================

clear all

% Ingresso dati Matrice U
U = [2 4 6 1 3; 0 7 0 2 2; 0 0 -3 -4 -12; 0 0 0 1 2; 0 0 0 0  -1 ];

% MATRICI PER TESTARE CONTROLLI

% U Rettangolare
% U = [0 4 6 1 3; 0 7 0 2 2; 0 0 -3 -4 -12; 0 0 0 1 2];

% U Singolare
% U = [0 4 6 1 3; 0 7 0 2 2; 0 0 -3 -4 -12; 0 0 0 1 2; 0 0 0 0  -1 ];

% U non triangolare superiore
% U = [2 4 6 1 3; 3 7 0 2 2; 0 0 -3 -4 -12; 0 0 0 1 2; 0 0 0 0  -1 ];

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

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

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

% calcola l'inversa di U con la funzione trisup_inv
Z = trisup_inv(U);

% calcola il prodotto Z*U
Id = Z*U;

% calcola l'inversa di U con la funzione Matlab inv
ZM = inv(U);

% Uscita risultati
fprintf('\n Matrice U \n');
for k = 1:n
    fprintf('\t %10g \t %8g \t %8g \t %8g \t %8g \n',U(k,:));
end

fprintf('\n Matrice U^-1 con function \n');
for k = 1:n
    fprintf('\t %10g \t %8g \t %8g \t %8g \t %8g \n',Z(k,:));
end

fprintf('\n Matrice U^-1 con Matlab inv \n');
for k = 1:n
    fprintf('\t %10g \t %8g \t %8g \t %8g \t %8g \n',ZM(k,:));
end

fprintf('\n Verifica che U^-1*U = I (U^-1 calcolata con function) \n');
for k = 1:n
    fprintf('\t %10g \t %8g \t %8g \t %8g \t %8.0f \n',Id(k,:));
end

% Scrive anche su file i risultati

file1= fopen('tabella.txt','w');
fprintf(file1, '\n Matrice U \n');
for k = 1:n
    fprintf(file1, '\t %10g \t %8g \t %8g \t %8g \t %8g \n',U(k,:));
end

fprintf(file1, '\n Matrice U^-1 con function \n');
for k = 1:n
    fprintf(file1, '\t %10g \t %8g \t %8g \t %8g \t %8g \n',Z(k,:));
end

fprintf(file1, '\n Matrice U^-1 con Matlab \n');
for k = 1:n
    fprintf(file1, '\t %10g \t %8g \t %8g \t %8g \t %8g \n',ZM(k,:));
end

fprintf(file1, '\n Verifica che U^-1*U = I (U^-1 calcolata con function) \n');
for k = 1:n
    fprintf(file1, '\t %10g \t %8g \t %8g \t %8g \t %8.0f \n',Id(k,:));
end

% calcola la norma della differenza
fprintf('\n \n norma della differenza delle matrici inverse  \n');
fprintf('\t %20.5e\n',norm(Z - ZM));

fprintf(file1, '\n \n norma della differenza delle matrici inverse  \n');
fprintf(file1, '\t %20.5e\n',norm(Z - ZM));
fclose(file1);