%% Cholesky-Verfahren
% Von der  symmetischen n\times n-MatrixA wird mit dem Cholesky-Verfahren eine Cholesky-Zerlegung A=LDL^T bestimmt.
%% Beispielmatrizen
% Beispiel 3.42 (Buch)
% Auf Grundlage der Matrix aus Beispiel 3.24 erzeugte Matrix A^TA
% Eine symmetrische, aber nicht positiv definite Matrix

bsp = 1;
A=example_A(bsp)

% 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=1;

% Erzeugung der Zerlegung

[L,D] = cholesky_zerlegung(A,stepwise)

%
%% Aktivität:
% Wenden Sie das Cholesky-Verfahren auf die drei Beispiele an.  Interpretieren Sie die Ergebnisse.
%
%

function A = example_A(k)  
    switch k
        case 1
            A = [2 6 -2; 6 21 0;-2 0 16];  
        case 2
            B =[ 2 -1 -3 3; 4 0 -3 1; 6 1 -1 6; -2 -5 4 1];
            A = B'*B;
        case 3
            A = [1 2 1; 2 3 3;1 3 2];
    end
end

function [L,D] = cholesky_zerlegung(A,stepwise)
    n=length(A);
    D=zeros(n);
    L=diag(ones(n,1));
    for j=1:n
        if stepwise==1 
            disp('-------------------');
            fprintf('\n \n%d. Spalte',j)
        end
        s = 0;
        for k=1:j-1
            s=s+L(j,k)*L(j,k)*D(k,k);
        end
        diagelem = A(j,j)-s;
        if diagelem < 10^(-5)*A(j,j) 
            fprintf('\nDiagonaleintrag hat den Wert %f.\n',diagelem)
            fprintf('Berechneter Diagonaleintrag ist nicht positiv (im Sinne des Abbruchkriteriums).\n') 
            fprintf('Die Berechnung der Cholesky-Zerlegung wird abgebrochen.\n')
            fprintf('Die Matrix ist nicht positiv definit.\n \n')
            L=NaN;    
            D=NaN;
            return
        end
        A(j,j)=diagelem;
        D(j,j)=diagelem;
        if stepwise==1                    
            D
        end        
        for i=j+1:n
            s=0;
            for k=1:j-1
                s=s+L(i,k)*D(k,k)*L(j,k);
            end
            L(i,j) = (A(i,j) - s)/D(j,j);           
            A(i,j) = (A(i,j) - s)/D(j,j);

            if stepwise==1                    
                fprintf('\n \n%d. Zeile',i)
                L
            end
        end   

    end
    disp('-------------------');
    if stepwise==1
        fprintf('\n\nEs folgen die berechneten Matrizen der Zerlegung:  \n \n')
    end
end