clear;
//
//    Esercitazione 9
//
// Studieremo il sistema di equazioni lineari
//
//    A x = b
//
// di ordine N definito da una matrice A notoriamente insidiosa
// per il procedimento di risoluzione basato su EGP e dalla colonna
// b = A ones(N,1).
//
// La soluzione del sistema è: x = ones(N,1).
//
// Cercheremo approssimazioni della soluzione del sistema sia con il
// procedimento basato su EGP che con quello basato sull'uso della
// fattorizzazione QR.
//
// ***************************
//
//     Procedure SI ed SA
//
// ***************************
//
function [x,info] = SI(T,c)
//
// Funzione SI: per T matrice N x N triangolare superiore invertibile e c colonna
// di N componenti, calcola la colonna soluzione del sistema di equazioni lineari
//
//       T x = c
//
// x: se T invertibile, la colonna x (di N numeri reali) approssima la soluzione
//    del sistema, altrimenti x = [].
// info = 0 se T non invertibile,
//      = 1 se T invertibile.
//
N = size(T,'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   info = 0; x = [];
   error("SI: elementi nulli sulla diagonale.");
   else
       info = 1;
       //
       // *** Calcolo soluzione
       //
       x = zeros(N,1);
       //
       x(N) = c(N) / T(N,N);
       for k = N-1:-1:1,
           x(k) = (c(k) - T(k,k+1:N)*x(k+1:N)) / T(k,k);
       end;
end;
endfunction
//
function [x,info] = SA(T,c)
//
// Funzione SA: per T matrice N x N triangolare inferiore invertibile e c colonna
// di N componenti, calcola la colonna soluzione del sistema di equazioni lineari
//
//       T x = c
//
// x: se T invertibile, la colonna x (di N numeri reali) approssima la soluzione
//    del sistema, altrimenti x = [].
// info = 0 se T non invertibile,
//      = 1 se T invertibile.
//
N = size(T, 'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   info = 0; x = [];
   error("SA: elementi nulli sulla diagonale.");
   else
       info = 1;
       //
       // *** Calcolo soluzione
       //
       x = zeros(N,1);
       //
       x(1) = c(1) / T(1,1);
       for k = 2:N,
           x(k) = (c(k) - T(k,1:k-1)*x(1:k-1)) / T(k,k);
       end;
end
endfunction
//
// ************************************************************************
//
N = input('Scegliere la dimensione del sistema:   N = ');
printf('\nN = %.0f\n', N);
A = -ones(N,N); A = tril(A,-1) + eye(A); A(:,N) = ones(N,1);
x = ones(N,1);
b = A*x;
//
// vediamo le proprietà di condizionamento del sistema:
COND_A = cond(A);
// calcolo un'approssimazione della soluzione con il procedimento che usa
// la fattorizzazione LR:
[S,D,P] = lu(A);
c = SA(S,P*b);
x_egpp = SI(D,c);
//
// errore commesso approssimando x con x_egpp:
EPS_xegpp = norm(x_egpp - x) / norm(x);
// interpreto x_egpp come soluzione del sistema perturbato A x = b + res_egpp
res_egpp = A*x_egpp - b;
// uso il Teorema di condizionamento, I:
//
//    EPS_d  <=  cond(A) * EPS_b
//
EPS_b_egpp = norm(res_egpp) / norm(b);
//
printf('\nEGPP:   EPS_d = %.3e  ,  cond(A) * EPS_b = %.3e\n',...
        EPS_xegpp,COND_A * EPS_b_egpp);
//
// calcolo adesso un'approssimazione della soluzione con il procedimento che usa
// la fattorizzazione QR:
[U,T] = qr(A);
c = U'*b;
x_qr = SI(T,c);
//
// errore commesso approssimando x con x_qr:
EPS_xqr = norm(x_qr - x) / norm(x);
// interpreto x_qr come soluzione del sistema perturbato A x = b + res_qr
res_qr = A*x_qr - b;
// uso il Teorema di condizionamento, I:
//
//    EPS_d  <=  cond(A) * EPS_b
//
EPS_b_qr = norm(res_qr) / norm(b);
//
printf('\nQR:     EPS_d = %.3e  ,  cond(A) * EPS_b = %.3e\n',...
        EPS_xqr,COND_A * EPS_b_qr);
abort;
//
// **************************************************************************
//
// Modello più dettagliato dei numeri di macchina:
//
//       F(base,precisione,esp_minimo,esp_massimo)
//
base = number_properties('radix');
precisione = number_properties('digits');
esp_minimo = number_properties('minexp');
esp_massimo = number_properties('maxexp');
//
printf('\nModello più dettagliato dei numeri di macchina:\n\n');
printf('      F(base,precisione,esp_minimo,esp_massimo)\n\n');
printf('      F(   %.0f,        %.0f,     %.0f,       %.0f)',...
           base,precisione,esp_minimo,esp_massimo);
//
DENORM = number_properties('denorm');
massimo = number_properties('huge');
minimo_normalizzato = number_properties('tiny');
minimo = number_properties('tiniest');
abort;
nearfloat('succ',0) == minimo
//
// Curiosità: come farsi mostrare esattamente da Scilab il valore della variabile
// massimo:
fd = mopen(TMPDIR+'/text.txt','wt');mfprintf(fd,'%.0f',massimo);mclose(fd);
edit(TMPDIR+'/text.txt')
//
