function rozklady(A)

% Tato funkce umožňuje provést maticový rozklad dle výběrů:
%   1. LU rozklad 
%         - vrací matice L a U
%   2. Choleského rozklad 
%         - vrací matice T a T'
%   3. Croutův rozklad
%         - vrací matice L a U
%
% Funkce ve všech třech případech přebírá jako argument jednu matici.
%
% Tento program vznikl jako zápočtový projekt do předmětu Numerické 
% výpočty I v podzimním semestru roku 2011.
%
% Vypracoval: Martin Nováček, 393995
% Vyučující: Mgr. Jiří Zelinka, Dr.
%
% Pro více informací, prosím, čtěte dokumentaci na:
%   www.martinnovacek.cz/edu/muni/rozklady
%


[m,n]=size(A);

if m*n==0
    error('Matice je prázdná! Zadejte neprázdnou matici.');     
end

if(m~=n)                                                      
    error('Matice není čtvercová! Zadejte čtvercovou matici.');             
end

disp('Vyberte typ rozkladu (1/2/3): ')                          % main title of the selection
disp('  1. LU rozklad')                                         % 1st choice for LU decomposition
disp('  2. Choleského rozklad')                                 % 2nd choice for Cholesky decomposition
disp('  3. Croutův rozklad')                                    % 3rd choice for Crout decomposition

choice = input('>> ','s');                                      % print ">>" and wait for user entered input

switch choice                                                   % switch which enable compute chosen decompositon
    case '1'                                    % LU decomposition
        disp('LU rozklad: ');                   % visual checkpoint for input to case 1
        
        L=zeros(n);                                              % preallocating L matrix
        U=zeros(n);                                              % preallocating U matrix
                                  
        for k=1:n                                               % for loop from 1 to n
            L(k,k)=1;                                           % gives ones to diagonal of L matrix
            for l=k+1:n                                         % for loop from k+1 to n
                if abs(A(k,k))<eps
                    error('Nedefinované dělení nulou. LU rozklad zadané matice nelze provést bez výměny řádku.');                    
                end
                L(l,k)=A(l,k)/A(k,k);                           % formula for under-diagonal elements of L matrix
                for j=k+1:n                                     % for loop from k+1 to n
                    A(l,j)=A(l,j)-(L(l,k)*A(k,j));              % recompute A matrix according to L matrix
                end
            end

            for j=k:n                                           % for loop from k to n
                U(k,j)=A(k,j);                                  % formula for upper-diagonal elements of U matrix
            end
        end
        
        disp('L =');                                            % print matrix L
        disp(L);
        disp('U =');                                            % print matrix U
        disp(U);
        disp('L*U =');                                          % print matrix L*U
        disp(L*U);                  
                             
    case '2'                                    % Cholesky decomposition
        disp('Choleského rozklad: ');           % visual checkpoint for input to case 2
        
        x=sum(sum(A==A'));                                      % summary of all elements
        y=m*n;                                                  % multiply m*n

        if x==y                                                 % if symmetric matrix
                             
           T = zeros(n);                                        % preallocating T(n,n) matrix
           
           for k=1:n                                            % for loop from 1 to n
               T(k,k)=sqrt(A(k,k)-(T(k,:)*T(k,:).'));           % formula for diagonal elements
               
               for l=(k+1):n                                    % for loop from i+1 to n (elements under diagonal)
                   if (abs(real(T(k,k)))<eps) && (abs(imag(T(k,k)))<eps)
                       error('Nedefinované dělení nulou. Choleského rozklad zadané matice nelze provést.');
                   end
                   T(l,k)=(A(l,k)-T(k,:)*T(l,:).')*(1/T(k,k));  % formula for under-diagonal elements
               end              
           end               
           
           disp('T =');                                         % print matrix T
           disp(T);
           disp('T'' =');                                       % print matrix T'
           disp(T.');
           disp('T*T'' =');                                     % print matrix T*T'
           disp(T*T.');
                       
        else                                                    % error message for non-symetric matrix
            error('Matice není symetrická! Zadejte symetrickou matici.')   
        end
                                   
    case '3'                                    % Crout decomposition
        disp('Croutův rozklad: ');              % visual checkpoint for input to case 3                               
        
        if ((sum(sum(triu(A,+2)))==0) && (sum(sum(tril(A,-2)))==0))
                                                                % check if A matrix is tridiagnal                                                                    
            L=zeros(n);                                         % preallocating L matrix
            U=eye(n);                                           % preallocating U eye-matrix
                        
            for i=2:n                                           % for loop from 2 to n - access to lower-diagonal elements of L matrix    
                L(i,i-1)=A(i,i-1);                              % computes lower-diagonal elements in L matrix
            end
            
            L(1,1)=A(1,1);                                      % A(1,1) value is copied to element L(1,1)
            
            for j=1:n-1                                         % for loop from 1 to n-1 - access to:
                if abs(L(j,j))<eps                              % 
                    error('Nedefinované dělení nulou. Croutův rozklad zadané matice nelze provést.');
                end                                             %
                U(j,j+1)=A(j,j+1)/L(j,j);                       %   >> upper-diagonal elements in U matrix
                L(j+1,j+1)=A(j+1,j+1)-(L(j+1,j)*U(j,j+1));      %   >> diagonal elements in L matrix
            end                        

            disp('L =');                                        % print matrix L
            disp(L);                                            
            disp('U =');                                        % print matrix U
            disp(U);
            disp('L*U =');                                      % print matrix L*U
            disp(L*U);
            
                   
        else error('Zadaná matice není třídiagonální!');
           
        end
        
    otherwise
        error('Chybný výběr rozkladu!');      
               
end     

%-------------------------------------------------------------------------%
% END - FUNCTION rozklady(A)                                              %
%-------------------------------------------------------------------------%