%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% forward.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = forward(L,b)

%   L : matrice triangolare inferiore NxN
%   b : termine noto N
%   Ly = b

N = length(b);
y = zeros(N,1);

y(1) = b(1)/L(1,1);
for i=2:N
  y(i) = ( b(i) - L(i,1:i-1)*y(1:i-1) ) / L(i,i);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% backward.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = backward(U,y)

%   U : matrice triangolare superiore NxN
%   y : termine noto N
%   Ux = y

N = length(y);
x = zeros(N,1);

x(N) = y(N)/U(N,N);
for i=N-1:-1:1
  x(i) = ( y(i) - U(i,i+1:N)*x(i+1:N) ) / U(i,i);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GaussTridiag.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = GaussTridiag(d,c,b)
%
% Risoluzione del sistema lineare Ax=b.
% A   matrice tridiagonale simmetrica definita positiva
%     (o a diagonale dominante) --> non necessita di pivoting
%     di A vengono assegnati solo la diagonale e le codiagonali.
% d   vettore n-dimensionale con la diagonale principale di A.
% c   vettore (n-1)-dimensionale con la sopra/sotto-diagonale di A.
% b   vettore colonna dei termini noti.
%

n=length(d);

for k = 1:n-1
   m      = - c(k)/d(k);
   d(k+1) = d(k+1) + m*c(k);
   b(k+1) = b(k+1) + m*b(k);
end

x=zeros(n,1);

x(n)=b(n)/d(n);

for k=n-1:-1:1
   x(k)=(b(k)-x(k+1)*c(k))/d(k);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solveLU.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = solveLU(A,b);

% Esegue la fattorizzazione LU senza pivoting
% Gli elementi di A in uscita contengono gli elementi di L e U
% A  matrice da fattorizzare
% b  termine noto del sistema lineare, puo` essere una matrice contenente
%    piu` vettori colonna
% x  soluzione del sistema lineare, un vettore colonna per ogni colonna di b

[n,n] = size(A);     % determina le dimensioni di A
[bm,bn] = size(b);

x = zeros(n,bn);

if bm~=n
   disp('Errore nella dimensione di b')
   x = NaN*ones(n,bn);
   return
end

A(:,n+1:n+bn)=b; % "Allega" alla matrice A i termini noti b
                 % in modo da modificare i termini noti con
                 % le stesse trasformazioni della matrice

for k=1:n-1
   if A(k,k) ~= 0
      A(k+1:n,k) = A(k+1:n,k)/A(k,k);
      A(k+1:n,k+1:n+bn) = A(k+1:n,k+1:n+bn) - A(k+1:n,k)*A(k,k+1:n+bn);
   else
      disp('Impossibile procedere con la fattorizzazione')
      disp('elemento pivot nullo')
      x = NaN*ones(n,bn);
      return
    end
end

% Risolve i sistemi lineari

for j=1:bn
    x(n,j) = A(n,n+j)/A(n,n);
    for i=n-1:-1:1
        x(i,j) = ( A(i,n+j) - A(i,i+1:n)*x(i+1:n,j) ) / A(i,i);
    end
end

A = A(:,1:n);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LUpiv.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,U,P] = LUpiv(A);

% A    matrice da fattorizzare
% L    matrice triangolare inferiore con 1 sulla diagonale
% U    matrice triangolare superiore
% P    matrice delle permutazioni
% piv  vettore delle permutazioni


[n,n] = size(A);     % determina le dimensioni di A
piv = 1:n;           % vettore pivot iniziale
for k=1:n-1

% determina l'elemento pivot nella k-esima colonna cercandolo
% tra gli elementi di indice riga compreso tra k e n.
% r parte da 1 in corrispondenza della k riga

   [maxv,r] = max(abs(A(k:n,k)));  
   q = r+k-1;               % determina la riga corrispondente a r
   piv([k q]) = piv([q k]); % scambia gli elementi k  p del vettore piv
   A([k q],:) = A([q k],:); % scambia le rige k e p della matrice A

   if A(k,k) ~= 0
      A(k+1:n,k) = A(k+1:n,k)/A(k,k);
      A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
    end

end
L = eye(n,n) + tril(A,-1);
U = triu(A);
P = eye(n);
P = P(piv,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LUpiv_mod.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,U,P] = LUpiv_mod(A);

% A    matrice da fattorizzare
% L    matrice triangolare inferiore con 1 sulla diagonale
% U    matrice triangolare superiore
% P    matrice delle permutazioni
% piv  vettore delle permutazioni


[n,n] = size(A);     % determina le dimensioni di A
piv = 1:n;           % vettore pivot iniziale
for k=1:n-1

% determina l'elemento pivot nella k-esima colonna cercandolo
% tra gli elementi di indice riga compreso tra k e n.
% r parte da 1 in corrispondenza della k riga

   [maxv,r] = max(abs(A(k:n,k)));  
   q = r+k-1;               % determina la riga corrispondente a r
   piv([k q]) = piv([q k]); % scambia gli elementi k  p del vettore piv
   A([k q],:) = A([q k],:); % scambia le rige k e p della matrice A

   if A(k,k) ~= 0
      A(k+1:n,k) = A(k+1:n,k)/A(k,k);
      A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
    end

    disp('Strike any key to continue')
    A
    pause

end
L = eye(n,n) + tril(A,-1);
U = triu(A);
P = eye(n);
P = P(piv,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fact_exe.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualizza la fattorizzazione PA=LU in forma razionale
% passo a passo della matrice A

format rat

A = [-2 4 -10 -1;
	4 -9 0 5;
       -4 5 -5 5;
      -8 8 -23 20]

[L,U,P] = LUpiv_mod(A);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% linsolv.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrice Tridiagonale differenze finite

echo on

N = 5;

pause % Matrice Tridiagonale differenze finite, simmetrica definita positiva
A = -diag(ones(1,N-1),-1)+2*diag(ones(1,N))-diag(ones(1,N-1),1)

pause % Termine noto random
b = rand(N,1);

pause % Soluzione del sistema lineare "secondo matlab"
x1 = A\b;

pause % E' bene controllare che il problema non sia mal condizionato
cond_a = condest(A)

pause % Se A e' singolare il suo determinante e' 0.0
det_A = det(A)


pause % Fattorizzazione LU con pivoting (non necessario P=I): L*U = P*A
[L,U,P] = lu(A)

pause % Soluzione di L*U*x=P*b
pause % Sostituzione in avanti L*y=P*b
y = forward(L,P*b);
pause % Sostituzione all'indietro U*x=y
x2 = backward(U,y);
pause % La soluzione e' la stessa!!!
norm(x1-x2,inf)

pause % La matrice data e' simmetrica definita positiva...
pause % Fattorizzazione di Cholesky: L*LT=A
LT = chol(A)
pause % Sostituzione in avanti L*y=b
y = forward(LT',b);
pause % Sostituzione all'indietro LT*x=y
x3 = backward(LT,y);
pause % La soluzione e' la stessa!!!
norm(x1-x3,inf)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% linsolvflops.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrice Tridiagonale differenze finite, confronto flops

echo off;

UP = 100;

for N=2:UP;

A = -diag(ones(1,N-1),-1)+2*diag(ones(1,N))-diag(ones(1,N-1),1);

b = rand(N,1);

flop1 = flops;
x1 = A\b;
flop2 = flops;
flop_backslash(N) = flop2 - flop1;

flop1 = flops;
%[L,U,P] = lu(A);
[L,U,P] = LUpiv(A);
flop2 = flops;
flop_factorization(N) = flop2 - flop1;

flop1 = flops;
y = forward(L,P*b);
flop2 = flops;
flop_forward(N) = flop2 - flop1;

flop1 = flops;
x2 = backward(U,y);
flop2 = flops;
flop_backward(N) = flop2 - flop1;

end

x = 1:UP;

plot(x,flop_backslash,'r-o',x,x.^3/3,'g:x')
title('Confronto flops backslash e n^3/3')

disp('Strike any key to continue')
pause

plot(x,flop_factorization,'r-o',x,x.^3/3,'g:x')
title('Confronto flops factorization PA=LU e n^3/3')

disp('Strike any key to continue')
pause

plot(x,flop_forward,'r-o',x,x.^2/2,'g:x')
title('Confronto flops forward substitution e n^2/2')

disp('Strike any key to continue')
pause

plot(x,flop_backward,'r-o',x,x.^2/2,'g:x')
title('Confronto flops backward substitution e n^2/2')

% -----------------------------------------------------

disp(' Fattorizzazione di Cholesky: L*LT=A')
pause

for N=2:UP;

A = -diag(ones(1,N-1),-1)+2*diag(ones(1,N))-diag(ones(1,N-1),1);

b = rand(N,1);

flop1 = flops;
LT = chol(A);
flop2 = flops;
flop_factorization(N) = flop2 - flop1;

flop1 = flops;
y = forward(LT',b);
flop2 = flops;
flop_forward(N) = flop2 - flop1;

flop1 = flops;
x3 = backward(LT,y);
flop2 = flops;
flop_backward(N) = flop2 - flop1;

end

plot(x,flop_factorization,'r-o',x,x.^3/6,'g:x')
title('Confronto flops factorization A=L*LT e n^3/6')

disp('Strike any key to continue')
pause

plot(x,flop_forward,'r-o',x,x.^2/2,'g:x')
title('Confronto flops forward substitution e n^2/2')

disp('Strike any key to continue')
pause

plot(x,flop_backward,'r-o',x,x.^2/2,'g:x')
title('Confronto flops backward substitution e n^2/2')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pivoting.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quando e' importante il pivoting?

echo off;
clc
disp('Senza pivoting')
disp('-----------------------------------------------------')
disp(' delta            x(1)                   x(2) ')
disp('-----------------------------------------------------')
for delta = logspace(-2,-18,9)
 A = [delta 1; 1 1];
 b = [1+delta;2];
 L = [1 0; 1/delta 1];
 U = [delta 1; 0 1-1/delta];
 y(1) = b(1);
 y(2) = b(2)-L(2,1)*y(1);
 x(2) = y(2)/U(2,2);
 x(1) = (y(1)-U(1,2)*x(2))/U(1,1);
 disp(sprintf(' %5.0e   %20.15f   %20.15f',delta,x(1),x(2)))
end


pause


disp(' ')
disp(' ')
disp('Con pivoting')
disp('-----------------------------------------------------')
disp(' delta            x(1)                   x(2) ')
disp('-----------------------------------------------------')
for delta = logspace(-2,-18,9)
 PA = [1 1; delta 1];
 b = [2;1+delta];
 L = [1 0; delta 1];
 U = [1 1; 0 1-delta];
 y(1) = b(1);
 y(2) = b(2)-L(2,1)*y(1);
 x(2) = y(2)/U(2,2);
 x(1) = (y(1)-U(1,2)*x(2))/U(1,1);
 disp(sprintf(' %5.0e   %20.15f   %20.15f',delta,x(1),x(2)))
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hilbert.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcola il numero di condizionameto della matrice di Hilbert,
% ne calcola l'inversa e verifica H*H^-1=I

echo on

N=10

pause % Costruzione della matrice di Hilbert di ordine N
H=hilb(N);

pause % Calcolo del numero di condizionamento in norma spettrale
lambda = eig(H);
cond_2_eig = max(abs(lambda))/min(abs(lambda)) % condizionameto_2 definizione

pause % Calcolo del numero di condizionamento in norma spettrale
cond_2_mat = cond(H) % condizionamento_2 calcolato da matlab

pause % Valutazione approssimata del numero di condizionamento in norma 1
cond_1_est = condest(H) % valutazione approssimata di cond_1
inv_cond_1_inv = 1.0/rcond(H) % valutazione approssimata di cond_1

I = diag(ones(N,1));

pause % Calcolo esplicito dell'inversa
H_1 = inv(H);
Ihilb = H_1*H;

pause % Norma 1 della differenza tra I ed (H^-1)*H
norm_I_Ihilb = norm(I-Ihilb,1)



pause % Calcolo dell'inversa esatta (formula esplicita)
inv_H = invhilb(N);
IHilb = inv_H*H;

pause % Norma 1 della differenza tra I ed (H^-1)*H
norm_I_InvHilb = norm(I-IHilb,1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hilbertsolv.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Effetti del mal condizionamento di un sistema lineare

echo off

UP = 8;

for i=1:UP

N = 2^i;

H = hilb(N);

x_ex = ones(N);

b = H*x_ex;

x1 = H\b;

norm_inf(i) = norm(x1-x_ex,inf);
fprintf(1,'%d\n',N)

end

disp('----------------------')
disp(' N      ||x-x_ex||_inf')
disp('----------------------')
for i=1:UP
fprintf(1,'%3d\t %.6E\n',2^i,norm_inf(i))
end

fprintf(1,'\n\n')
disp('Strike any key to continue')
pause % Termine noto lievemente perturbato

B = 1E-8;

for i=1:UP

N = 2^i;

H = hilb(N);

x_ex = ones(N);

b = H*x_ex;

norm_b(i) = norm(b,inf);

for j=1:N
b(j) = b(j)+B*sin(2*pi*j/N);
end

x1 = H\b;

norm_inf(i) = norm(x1-x_ex,inf);
fprintf(1,'%d\n',N)

condit(i) = cond(H,inf);

end

disp('-------------------------------------------------------------')
disp(' N   (||x-x_ex||/||x||)_inf  (||db||/||b||)_inf    cond(H)')
disp('-------------------------------------------------------------')
for i=1:UP
fprintf(1,'%3d \t %.6E   \t %.6E  \t %.6E\n',...
	2^i,norm_inf(i),B/norm_b(i),condit(i))
end