%Es werden die Verfahren von Crout und Doolittle auf eine Tridiagonalmatrix
%angewendent. Numerik SS 98, Aufgabe 31, Blatt 16.
%Alle Nichtdiagonalelemente von LD und von R
%sind dem Betrage nach <1 falls A diagonaldominant. 
%Die Berechnungsreihenfolge kann man in den Matrizen Rf (Crout)
%und RfD (Doolittle) ablesen, zuerst 1, dann -1 dann 2 dann -2 etc.
%Funktioniert gut, 26. 4. 1998
M=5;
L=zeros(M,M);
Rf=L; RfD=L;%in Rf, RfD wird die Berechnungsreihenfolge festgehalten
R=L;
LD=L;RD=L; %L,R für Crout, LD, RD für Doolittle
%A=unifrnd(-1,1,M,M);
%A=A+M*eye(M); %diagonaldominant
%A=10*A; A=round(A);
%B=diag(diag(A));
%A=A-B;A=A+diag(sum(abs(A'))); %zeilendiagonaldominant
%A=[9 -36 30;-36 192 -180;30 -180 180]; %Test aus dem Buch, S.136
A=[10 -8 0  0  0
   3   7 4  0  0
   0   1 3 -1  0
   0   0 1  4 -1 
   0   0 0  8 -9];
for K=1:M
R(K,K)=1;LD(K,K)=1;
L(K:M,K)=A(K:M,K)-L(K:M,1:K-1)*R(1:K-1,K);
Rf(K:M,K)=K*ones(size(L(K:M,K)));
RD(K,K:M)=A(K,K:M)-LD(K,1:K-1)*RD(1:K-1,K:M);
RfD(K,K:M)=K*ones(size(RD(K,K:M)));
R(K,K+1:M)=(A(K,K+1:M)-L(K,1:K-1)*R(1:K-1,K+1:M))/L(K,K);
Rf(K,K+1:M)=-K*ones(size(R(K,K+1:M)));
LD(K+1:M,K)=(A(K+1:M,K)-LD(K+1:M,1:K-1)*RD(1:K-1,K))/RD(K,K);
RfD(K+1:M,K)=-K*ones(size(LD(K+1:M,K)));
end; %for K
%Rf = (Reihenfolge bei Crout: 1 -1 2 -2 etc)
%     1    -1    -1    -1    -1
%     1     2    -2    -2    -2
%     1     2     3    -3    -3
%     1     2     3     4    -4
%     1     2     3     4     5
%RfD =  (Reihenfolge bei Doolittle: 1 -1 2 -2 etc)
%     1     1     1     1     1
%    -1     2     2     2     2
%    -1    -2     3     3     3
%    -1    -2    -3     4     4
%    -1    -2    -3    -4     5
