% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%                                                                       % 
% Loesung der Aufgabe 12, Matlab Kurs WS 06/07                          % 
%                                                                       % 
% Dieses Programm liest eine Matrix A und einen Vektor b ein und loest  %
% das Gleichungssystem Ax=b mit Hilfe der Choleskyzerlegung             %
%                                                                       % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

disp('Implementation des Cholesky-Algorithmus fuer ')
disp('positiv definite Matrizen A fuer das LGS. A x = b. ')
disp(' ')
disp(' ')
disp('Bitte quadratische, symmetrische, ')
A = input('positiv definite Matrix A eingeben. A = ');
disp(' ')
disp(' ')
b = input('Bitte Loesungsvektor b eingeben. b = ');
disp(' ')
disp(' ')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Dimension von A:
[n,m] = size(A);
if n~=m
     error('Matrix A ist nicht quadratisch !');
end

% Ist die Matrix A symmetrisch ?
if A == A'
    disp('Matrix ist symmetrisch!')
else
    error('Matrix A ist nicht symmetrisch !');
end

% Ist die Matrix A positiv definit ?
% Anm.: Eigenwerte(A) > 0 und A symmetrisch => A pos definit
if min(eig(A)) < 0
    error('Matrix A ist nicht positiv definit !');
end

% Ist b ein Spaltenvektor ?
[q,p] = size(b);
if q~=n
     error('b hat die falsche Dimension!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Cholesky-Zerlegung A = L L'
%  L ist untere Dreiecksmatrix

%  Cholesky-Verfahren:
%  (erzeugt die L-Matrix, man geht spaltenweise vor)
L=zeros(n,n);

for j=1:n
    summe = 0;
    for k=1:j-1
        summe = summe + (L(j,k))^2;
    end
    L(j,j) = sqrt(A(j,j) - summe);
    for i=j+1:n
        summe = 0;
        for k=1:j-1
            summe = summe + L(i,k)*conj(L(j,k));
        end
        L(i,j) = (A(i,j) - summe) / L(j,j);
    end
end

% Ausgabe der gewonnenen Matrix L
disp ('Die untere Dreiecksmatrix L:');
disp(L);
R = L';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Loesen des Gleichungssystems

% Schritt 2: Vorwaertssubstitution Ly=b 
y = zeros(n,1); 
for j=1:n
	summe = 0;
	for k=1:j-1
		summe = summe + L(j,k) * y(k);
	end;
	y(j) = (b(j) - summe) / L(j,j);
end


% Schritt 2: Rueckwaertssubstitution Rx=z mit R=L' 
x = zeros(n,1);
for j=n:-1:1
	summe = 0;
	for k=j+1:n
		summe = summe + R(j,k) * x(k);
	end;
	x(j) = (y(j)- summe) / R(j,j);
end



disp('Loesungsvektor x=');
disp(x);