%% Gauß-Elimination
% Das Gleichungssystem Ax=b wird gelöst mittels Gauß-Elimination mit/ohne Zeilenskalierung und mit/ohne Spaltenpivotisierung. Die Zwischenresultate in jedem Schritt des Eliminationsverfahren können angezeigt werden.
%%  Beispielmatrizen
% Beispiel 3.24 (Buch)
% Ohne Pivotisierung wird Pivotelement gleich Null erzeugt
% Ohne Pivotisierung nicht stabil
% Matrix mit Rangdefizit (= kein voller Rang)
% Das Beispiel bsp (1-4) wird gewählt:

bsp=1;
format shorte;
[A,b] = example(bsp)
format short;

% Zeilenskalierung
% Der Parameterwert scaling soll auf 0 (keine Skalierung) oder 1 (mit Skalierung) gesetzt werden.

scaling=0;

% Spaltenpivotisierung
% Der Parameterwert pivoting soll auf 0 (keine Pivotisierung) oder 1 (mit Pivotisierung) gesetzt werden.

pivoting=1;

% Schrittweiser Durchlauf
% Wenn das Verfahren nicht schrittweise durchlaufen werden soll, ist der folgende Parameter stepwise auf 0 zu setzen. Ansonsten ist der Wert 1 zuzuweisen.

stepwise=0;

%
%% Das Gauß-Eliminationsverfahren \bigl([A\mid b]\bigr)
% Die Gauß-Elimination wird durchgeführt. Die Matrix A und  rechte Seite b werden in eine n\times (n+1)-Matrix Ab geschrieben  und synchron modifiziert.

Ab = [A b]
n=length(A);
Konditionszahl= cond(Ab(:,1:n),'inf')

% Zeilenskalierung
% Die Zeilenskalierung verbessert die Kondition des Problems bzgl. der Maximumnorm.

Ab=zeilenskalierung(Ab,scaling);
Konditionszahl= cond(Ab(:,1:n),'inf')

% Erzeugung der Dreiecksgestalt
% Hier wird nun die Dreiecksgestalt von A erzeugt, eventuell unter Zurhilfenahme von Pivotisierung.

for i=1:length(Ab)-1 
    Ab=pivotisierung(Ab,i,pivoting,stepwise);
    if test_diagonalelement(Ab,i)
       return
    end
    Ab=annullierung(Ab,i,stepwise);    
end

% Dies ist nun die Matrix in Dreiecksgestalt plus die transformierte rechte Seite.

Ab

% Rückwärtseinsetzen
% Es wird das Gleichungssystem durch Rückwärtseinsetzen gelöst.

 disp('-------------------');
x=rueckwaerts_einsetzen(Ab,stepwise)

%
%% Aktivität:
% Wählen Sie die unterschiedlichen Beispiele und variieren Sie die Parameterwerte  scaling und pivoting.
% Interpretieren Sie die Ergebnisse.
%

function [A,b] = example(k)  
    switch k
        case 1
            A =[ 2 -1 -3 3; 4 0 -3 1; 6 1 -1 6; -2 -5 4 1];
            b= [ 1; -8; -16; -12];
        case 2
            A = [1 2 3 4; -2 1 0 -1;1 -3 -3 5; 0 10 21 24];
            b= [ 24; 4; -12; 151];
        case 3
            A = [10^(-17) 3 5; 6 2 2; 1 3 4];
            b = [11; 9; 10.5];                  
        case 4
            A = [1 3 8 4;2 6 22 15;-1 -3 1 -2;3 9 29 22];
            b = ones(4,1); 
    end 
end

function [A] = zeilenskalierung(Ab,scaling)
    
    n=length(Ab)-1;
    
    if scaling==1
    
        for i=1:n
            d=sum(abs(Ab(i,1:n)));
            Ab(i,:)=(1/d)*Ab(i,:);
        end
        fprintf('Skalierte Matrix')
        Ab
    else
        fprintf('Zeilenskalierung ausgeschaltet')
    end
    A=Ab;
end

function [A] = pivotisierung(Ab,i,pivoting,stepwise)
     disp('-------------------');
% Die Pivotisierung, wie sie hier durchgeführt wird, sucht 
% in der Spalte vom behandelten Diagonalelement bis zur letzten 
% Zeile den betragsmäßig größten Eintrag    
    n=length(Ab)-1;
    if pivoting==1                
        if i==1        
            fprintf('Pivotisierung ist eingeschaltet.\n');
        end
% Matlab liefert hier das Maximum und die Zeile des Maximums.        
        [amax,pos] = max(abs(Ab(i:n,i)));
        if amax==0
            fprintf('Pivotelement ist 0!\n');
            pos=1;
        end
% Der lokale Index pos liegt im Bereich 1...n-i+1 und muss
% in den globalen Index p im Bereich i...n umgerechnet werden.
        p = pos+i-1;
        if p>i
% Im Falle, dass das Diagonalelement nicht das größte Element ist,
% wird die Zeile mit dem behandelten Diagonalelement mit der Zeile 
% in der das betragsmäßige größte Element steht, vertauscht.            
           Ab([p i],:) = Ab([i p],:);
           if stepwise==1            
               fprintf('Zeile %d vertauscht mit Zeile %d',i,p)
               Ab
           end
        end
    else
        if i==1        
            fprintf('Pivotisierung ist ausgeschaltet.\n');
        end
    end
    
    A=Ab;
end

function A = annullierung(Ab,i,stepwise)
    n=length(Ab)-1;
    m=n+1;
    for k=i+1:n
        faktor = Ab(k,i)/Ab(i,i); % = L(k,i) aus LR-Zerlegung
        Ab(k,i)=0;
        Ab(k,i+1:m)=Ab(k,i+1:m)-faktor*Ab(i,i+1:m);
        
        if stepwise==1
            fprintf('\nAnnulierung des Elements (%d,%d)\n',k,i)
            Ab
        end
    end
    A=Ab;
end

function testresult = test_diagonalelement(Ab,i)
% Wenn das Diagonalelement gleich 0 ist, kann der Algorithmus nicht weiter 
% durchgeführt werden, da es sonst zu einer Division durch käme. Sollte 
% die Pivotisierung nicht verwendet worden sein, kann sie in den Fällen, 
% in denen die Matrix regulär ist, Abhilfe schaffen.    
    n=length(Ab)-1;
    if Ab(i,i)~=0
        testresult=  false;
    else
        testresult=  true;
        fprintf('Das Diagonalelement Ab(%d,%d) ist 0\n',i,i);
        amax =max(abs(Ab(i:n,i)));
        if amax==0 
            fprintf('Das Gleichungssystem besitzt keine eindeutige Lösung\n')
        else
            fprintf('Pivotisierung kann helfen.\n')
        end
        fprintf('Der Algorithmus bricht ab.')
        Ab
    end
end

function [x] = rueckwaerts_einsetzen(Ab,stepwise)
    m=length(Ab);
    n=m-1;    
    x=zeros(n,1);
    for j=n:-1:1
        if j<n
            s=Ab(j,j+1:n)*x(j+1:n,1);
        else
            s=0;
        end       
        x(j,1) = (Ab(j,m)-s)/Ab(j,j);
        if stepwise==1
            fprintf('Bestimmung von x(%d)',j)
            x 
        end
    end
end