disp(' ');
disp(' spectral_decomp.m ');
disp(' version 1.1  March 30, 2010 ');
disp(' ');
disp(' By Tom Irvine ');
disp(' ');
disp(' Spectral Decomposition:  A=LDLT ');
disp(' ');
disp(' The input matrix must be symmetric. ');
disp(' ');
%
clear AA;
clear A;
clear G;
clear L;
clear U;
clear D;
%
disp(' ');
disp(' Select file input method ');
disp('   1=file preloaded into Matlab ');
disp('   2=Excel file ');
file_choice = input('');
%
disp(' ');
disp(' Mass Matrix ');
%
if(file_choice==1)
        AA = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%
        AA = xlsread(xfile);
%
end
%
A=AA;
nn=size(A);
n=nn(1,1);
%
L=zeros(n,n);
U=zeros(n,n);
D=zeros(n,n);
%
%  Triangular Decomposition
%
L=eye(n,n);
%
i=1;
%
U(1,1)=A(1,1);
for i=2:n
    L(i,1)=A(i,1)/U(1,1);
    U(1,i)=A(1,i);
end
%
for i=2:n
%
    sum=0.;
    for k=1:i-1
        sum=sum+U(k,i)*L(i,k);
    end
    U(i,i)=A(i,i)-sum;
%
    for j=i+1:n
%        
        sum=0.;
        for k=1:j-1
            sum=sum+U(k,i)*L(j,k);
        end
        L(j,i)=(A(j,i)-sum)/U(i,i);
%
        sum=0.;
        for k=1:j-1
            sum=sum+U(k,j)*L(i,k);
        end
        U(i,j)=(A(i,j)-sum);
    end
end
%
% Solve for Diagonal Matrix D
%
for i=1:n
    D(i,i)=U(i,i);
end
L
D
LT=L';
LT