% Funkcija demonstrira postupak nalazenja ortogonalno-trougaone 
% dekompozicije realne matrice A, koriscenjem Givensovih matrica 
% rotacije
% Algoritam se sastoji u uzastopnom mnozenju matrice A sa leve 
% strane matricama rotacije, sve dok ne svedemo matricu A na 
% gornje trougaonu matricu R tj.
%       Un*Un-1*..*U1*A=R
% Odavde je A=U1'*U2'*..Un'*R
% tj. Q=U1'*...*Un' (Ovo je ortogonalna matrica kao proizvod 
% ortogonalnih matrica  

function [Q,R]=GivensQR(A)
  % Odredjujemo dimenziju matrice A
  [m,n]=size(A);
    
   % Matricu Q postavljamo na jedinicnu  
  Q=eye(m,n);

  % Anuliramo jedan po jedan element ispod dijagonale 
  % matrice A 
  for k=1:n-1
    for l=k+1:n
        if A(l,k)~=0 
            % Izracunavamo elemente matrice rotacije
            % koja anulira A(l,k).
            % Obratiti paznju na indekse u ovim formulama koji se razlikuju
            % od indeksa u GivensHessenberg1.m. Razlog je sto se tamo
            % anuliraju svi elementi ispod poddijagonale (svodjenje u
            % slicnu matricu u gornjeHesenbergovoj formi), dok se ovde
            % anuliraju svi elementi ispod dijagonale (svodjenje na matricu
            % koja je gorenjetrougaona i nije slicbna polaznoj).
            alfa=A(k,k)/(sqrt(A(k,k)^2+A(l,k)^2));
            beta=-A(l,k)/(sqrt(A(k,k)^2+A(l,k)^2));
            

            % Gradimo matricu rotacije umetanjem u jedinicnu
            U=eye(size(A));
            U(k,k)=alfa;
            U(l,l)=alfa;
            U(k,l)=-beta;
            U(l,k)=beta;
                       
            % Anuliramo odgovarajuci elemenat matrice A 
            % (matrica nece biti slicna polaznoj, te se nemnozi sa obe
            % strane)
            A=U*A;
            % Akumuliramo proizvod matrica rotacije u Q
            Q=Q*U';
            
        end
    end
end

% Matrica R je dobijena gornjetrougaona matrica
R=A;

