%====================================================================
% Script per determinare la soluzione di un sistema idraulico
% con il metodo di fattorizzazione di Cholesky
% Necessita delle function triinf e trisup
%====================================================================

disp('APPLICAZIONE METODO DI CHOLESKI')

% Dati sistema letti da file idra.dat

Ab = load('idra.dat');
A = Ab(:,1:end-1);
b = Ab(:,end);


%%%% ISTRUZIONI PER FARE TEST  (matrice NON definita positiva)
%% A=-A;
%% b=-b;


%%%% ISTRUZIONI PER FARE TEST  (matrice NON simmetrica)
%% A (2,1) = 5;

% Modifica sistema originale per avere la possibilita' che la sua matrice
% sia definita positiva
A1=-A;
b1=-b;

% Controlli generali
[mA,nA] = size(A1);
if (mA~=nA)
    error('La matrice non e'' quadrata')
end

if ~isequal(A1',A1)  
    error('La matrice non e'' simmetrica')
end

if ~all(diag(A1)>0)  
% OPPURE
% if ~all(diag(A1)==abs(diag(A1)))
% OPPURE ....
    error('La matrice non puo'' essere definita positiva')
end

if length(b1) ~= nA
    error('Le dimensioni della matrice e del termine noto non sono compatibili')
end

% Risoluzione

% fattorizzazione di Cholesky
[L,flag]=chol(A1,'lower');
if flag ~= 0  
    error('La matrice non e'' definita positiva') 
end

% calcolo di L * L^T

LLT = L*L';

% Definisce il formato di visualizzazione
format long e

% Salva nel file matrici.txt le due array A1 e LLT
save('matrici.txt','-ascii','-double','LLT','A1')

% Visualizzazione matrici A1 ed LLT
disp('matrice del sistema')
disp(A1)
disp('matrice risultante dal prodotto L*L''')
disp(LLT)

% Uso la fattorizzazione per risolvere il sistema
% Sistema LL^T x = b

% Risolvo  L y = b
y = triinf(L, b1);

% Risolvo L^T x = y
x = trisup(L', y);

% Risolvo A x = b con comando Matlab
x1 = A1\b1;

% Uscita risultati

disp('termine noto')
disp(b1)
disp('dimensione')
disp(mA)
disp('matrice L')
disp(L)
disp('soluzione con fattorizzazione')
disp(x)
disp('soluzione con comando Matlab')
disp(x1)
format
disp('norma della differenza vettori soluzione')
disp(norm(x-x1))

% ritorna al format di default
format