clc
clear
format long; % increase precision
n=64;%Set matrix size 8
A=-1+2.0*rand(n,n);
R=hess(A); %Perform the built-in hessenberg reduction (reference)
R1=hessen(A); %Perform the hessenberg reduction (column householder)
R2=hessen1(A); %Perform the hessenberg reduction (givens rotation)
figure; matrixplot(R1);
figure; matrixplot(R2);
%Test definition of Hessenberg
norm1=norm((R-hess(A))/R) %Just for sanity, built-in is assumed to work
norm2=norm((R-R1)/R)
norm3=norm((R-R2)/R)
%Test eigen values
orig_eig=eig(A);
eval=eig(R); %E-val value of built-in hessenberg reduced matrix
eval1=eig(R1); %E-val of hessenberg reduced matrix(column householder)
eval2=eig(R2); %E-val of hessenberg reduced matrix(givens rotation)
eig_norm=norm((eval-eig(R))/eig(R))
eig_norm1=norm((eval1-eval)/eval)
eig_norm2=norm((eval2-eval)/eval)
tElapsed=zeros(1,50); %Create a time elapsed vector
tElapsed1=zeros(1,50); %Create a time elapsed vector
tElapsed2=zeros(1,50); %Create a time elapsed vector
tElapsed3=zeros(1,50); %Create a time elapsed vector
x=logspace(1, 2, 50); %create logspace for the x-axis
for i=1:50 %loop through the order of matrix
    n=round(x(i)); %round the logspace to get the closest integer order
    %Initialize coefficient matrix G with a random value
    G=-1+2.0*rand(n,n); % Random matrix of order n between -1 and 1
    tStart = tic; %Start the timer
    result=hess(G); %Use the built-in hess function
    tElapsed(i) = toc(tStart); %stop the timer and capture elapsed time
    tStart = tic; %Start the timer
    result=hessen(G); %Use the column householder function
    tElapsed1(i) = toc(tStart); %stop the timer and capture elapsed time
    tStart = tic; %Start the timer
    result=hessen1(G);
    tElapsed2(i) = toc(tStart); %stop the timer and capture elapsed time
    
end
figure; loglog(x, tElapsed, 'r-', x, tElapsed1, 'g-', x, tElapsed2, 'b-'); %Plot the time vs time elapsed
grid on;
xlabel('Matrix size');
ylabel('Time of computation');
legend('Built-in', 'HR Column', 'Givens HR');