%
%
%
% F.Marcuzzi  13/05/2005
%
clc
disp('Esercitazione sui sistemi con matrice sparsa :');
disp(' ');
disp('(ad ogni interruzione del programma, premi un tasto per continuare)');
disp(' '); disp(' ');
disp('Prenderemo in considerazione la fattorizzazione LU di una matrice.');
disp(' ');
disp('Sappiamo già che la fattorizzazione LU di una matrice nxn richiede in generale circa n^3 flops.');
disp(' ');
disp('Quando però ci sono elementi nulli nella matrice, si creano notevoli possibilità di ridurre');
disp('il numero di operazioni necessarie. ');
pause

clc
disp('Ad es. consideriamo la matrice A (vedi Figura 1): ');
n=100; A=rand(n); A=(A+A')./2+5*diag(diag(ones(n))); I=find(A<0.85); A(I)=0; origA=A; 
figure(1); clf, spy(A); title(['matrice A     (isp = ' num2str( nnz(A)/(n^2) ) ')']);
figure(2); clf, spy(A); title(['matrice A     (isp = ' num2str( nnz(A)/(n^2) ) ')']); pause
disp(' ');
disp('e vediamo il numero di flops necessari alla fattorizzazione LU');
disp(' ');
disp('senza tenere conto degli zeri:');
A=origA;
nf = 0;
for i=1:n-1
  A(i+1:n,i) = A(i+1:n,i)./A(i,i);
  A(i+1:n,i+1:n) = A(i+1:n,i+1:n) - A(i+1:n,i)*A(i,i+1:n);
  nf = nf + (n-i) + 2*(n-i)^2;
  figure(2); clf, hold on, spy(A); 
             plot(i,i,'ro')
             title(['matrice L+U     (isp = ' num2str( nnz(A)/(n^2) ) ')']);
  tic, while toc<0.5  figure(2); drawnow, end;
end;
disp(['flops necessari = ' num2str(nf)]);
disp(' ');
disp('ed, invece, tenendo conto degli zeri:');
Az=origA;
nf = 0;
for i=1:n-1
  I=find([zeros(i,1); Az(i+1:n,i)]);
  if ~isempty(I)
    Az(I,i)=Az(I,i)./Az(i,i);
    Az(I,I)=Az(I,I) - Az(I,i)*Az(i,I);
    nf = nf + length(I) + 2*length(I)^2;
  end;
end;
disp(['flops necessari = ' num2str(nf)]);
pause
disp(' '); disp(' '); disp(' '); disp(' ');
disp('E'' da notare,pero'', che la fattorizzazione ha prodotto due fattori (L e U) che complessivamente');
disp('possiedono più elementi diversi da zero di quanti ne avesse la matrice di partenza:');
disp(['elementi non-zero di A   = ' num2str(nnz(origA))]);
disp(['elementi non-zero di L+U = ' num2str(nnz(A))]);
disp(' ');
disp('questo fenomeno è chiamato "fill-in"');
disp(' ');
disp('Il fill-in è un fenomeno da tenere in grande considerazione, in quanto la sua presenza produce');
disp('un aumento considerevole della quantità di flops necessari durante la fattorizzazione.');
disp(' ');
disp('Per rendersene conto, basta guardare le ultime righe di codice scritte qui sopra e notare ad es. che');
disp('il prodotto esterno "Az(I,i)*Az(i,I)" richiede una frazione considerevole del numero di flops complessivi');
disp(' ');
pause
pause
pause

clc
disp('Il fill-in dipende fortemente dall''ordine con cui vengono eliminate le variabili');
disp('durante il processo di fattorizzazione.');
disp('Vediamo a titolo di esempio la seguente matrice (vedi Figura 1) :');
n=12; A=diag(diag(ones(n))); A(n,1:n-1)=repmat(0.1,1,n-1); A(1:n-1,n)=repmat(0.1,n-1,1); Afreccia=A; 
figure(1); clf, spy(A); title('matrice A'); pause,
disp(' ');
disp('e vediamo il numero di flops necessari alla fattorizzazione LU :');
Az=Afreccia;
nf = 0;
for i=1:n-1
  I=find([zeros(i,1); Az(i+1:n,i)]);
  if ~isempty(I)
    Az(I,i)=Az(I,i)./Az(i,i);
    Az(I,I)=Az(I,I)-Az(I,i)*Az(i,I);
    nf = nf + length(I) + 2*length(I)^2;
  end;
end;
disp(['flops necessari = ' num2str(nf)]);
disp(' ');
figure(2); clf, spy(Az); title('matrice L+U');
pause
disp('mentre, utilizzando l''ordinamento esattamente opposto (vedi Figura 3) : ');
Az2=Afreccia(n:-1:1,n:-1:1);  
figure(3); clf, spy(Az2); title('matrice A riordinata nel senso opposto'); pause
nf = 0;
for i=1:n-1
  I=find([zeros(i,1); Az2(i+1:n,i)]);
  if ~isempty(I)
    Az2(I,i)=Az2(I,i)./Az2(i,i);
    Az2(I,I)=Az2(I,I)-Az2(I,i)*Az2(i,I);
    nf = nf + length(I) + 2*length(I)^2;
  end;
end;
disp(['flops necessari = ' num2str(nf)]);
figure(4); clf, spy(Az2); title('matrice L+U');
pause
disp(' '); disp(' '); disp(' '); disp(' ');
disp('Vediamo quale è il fill-in risultante :');
disp(['elementi non-zero di A                                      = ' num2str(nnz(Afreccia))]);
disp(['elementi non-zero di L+U secondo l''ordinamento diretto di A = ' num2str(nnz(Az))]);
disp(['elementi non-zero di L+U secondo l''ordinamento opposto      = ' num2str(nnz(Az2))]);
disp(' ');
disp('Notare che utilizzando l''ordinamento opposto, i fattori L e U risultano essere addirittura densi!');
disp(['(vedi Figura 4 ed osserva che ' num2str(nnz(Az2)) ' = ' num2str(n) 'x' num2str(n) ' )']);
pause
close all;

clc
pause
disp('In generale è molto difficile sapere a-priori quale è l''ordinamento ottimale, ovvero');
disp('quello che crea meno fill-in in assoluto. ');
disp('E'' un problema di complessità "NP-complete" (cresce esponenzialmente con la dimensione della matrice).');
disp(' ');
disp('Per questo motivo si utilizzano di solito degli algoritmi di pre-ordinamento di tipo euristico,');
disp('basati sulla Teoria dei Grafi (perchè ad una matrice sparsa viene generalmente associato un grafo,');
disp('i cui vertici sono le variabili ed i cui lati sono gli elementi non-nulli).');
disp(' ');
disp('Di questi ne citiamo alcuni (che ritroviamo implementati nel seguito):');
disp(' ');
disp('1) l''algoritmo di Cuthill-McKee (symrcm), il cui obbiettivo è quello di rendere minima la banda');
disp('   della matrice riordinata.');
disp(' ');
disp('2) l''algoritmo del grado minimo (colmmd e symmmd), che agisce scegliendo per l''eliminazione la colonna');
disp('   (o riga) di grado minimo (intendendo per grado il numero di lati uscenti dal vertice corrispondente alla');
disp('   colonna nel grafo associato).');
disp(' ');
disp('3) l''algoritmo del grado minimo approssimato (COLAMD), che aggiorna il grado dei vertici rimanenti in modo');
disp('   approssimato e rappresenta il miglior risultato ottenuto ad oggi in questa classe di algoritmi.');
disp('   Esso e'' disponibile gratuitamente al sito http://www.cise.ufl.edu/research/sparse/colamd/ ');
pause

clc
scelta=1; 
while scelta~=0 
  clc
  n=input('inserisci il numero di righe (colonne) "n" della matrice : ');
  s=input('inserisci un numero da 0 a 1 che indichi l''indice di densità della matrice\n1=matrice densa, 0=matrice nulla : ');
  A=rand(n); A=(A+A')./2+5*diag(diag(ones(n))); I=find(A<(1-s)); A(I)=0; origA=A; 
  stessamat=1;
  while stessamat==1
    sA=sparse(origA);
    disp('Algoritmi di ordinamento:');
    disp('1) colmmd');
    disp('2) symmmd');
    disp('3) symrcm');
    disp('4) colperm');
    disp('5) dmperm');
    disp('6) colamd (COLAMD)');
    disp('7) symamd (COLAMD)');
    disp('8) nessun pre-ordinamento');    
    m=input('scegli l''algoritmo (0 per terminare) : ');
    switch m
      case 0, break;
      case 1, order=colmmd(sA);
      case 2, order=symmmd(sA);
      case 3, order=symrcm(sA);
      case 4, order=colperm(sA);
      case 5, order=dmperm(sA);
      case 6, order=colamd(sA);
      case 7, order=symamd(sA);
      case 8, order=1:n;
      otherwise error('numero di algoritmo inesistente !');
    end;
    figure(1); clf, spy(A); title('matrice A'); pause,
    reordA=sA(order,order); figure(2); spy(reordA); title('matrice A secondo l''ordinamento scelto');
    nf = 0;
    for i=1:n-1
      I=find([zeros(i,1); reordA(i+1:n,i)]);
      if ~isempty(I)
        reordA(I,i)=reordA(I,i)./reordA(i,i);
        reordA(I,I)=reordA(I,I)-reordA(I,i)*reordA(i,I);
        nf = nf + length(I) + 2*length(I)^2;
      end;
      figure(3); clf, spy(reordA,'m'); title('matrice L+U'), hold on, spy(sA(order,order));
                 plot(i,i,'ro')
      tic, while toc<0.5 figure(3); drawnow, end;
    end;
    disp(['flops necessari = ' num2str(nf)]);
    disp(['elementi non-zero di L+U = ' num2str(nnz(reordA))]);
    stessamat=input('come vuoi continuare : 1-stessa matrice, 2-cambia matrice, 0-uscire   --> ');
    if stessamat==0  scelta=0; end;
  end;
end;