echo on
%This example shows the QR decomposition of a 4 by 4 magic matrix
%This creates a 4 by 4 magic matrix
A= magic(4)

I = eye(4)

%The following find the vector
%v = fist column vector - norm(fist column vector)*[1 0 0 0]'
%note that I(:,1)=[1 0 0 0]'

v = A(:,1) - norm(A(:,1))*I(:,1)

H_1 = householder(v)
A_1 = H_1*A

% I(2:4,2)=[1 0 0]'
% v = [ A_1(1,2) A_1(2,2) A_1(2,3) A_1(2,4)] 
%    -[ A_1(1,2) norm([A_(2,2) A_1(2,3) A_1(2,4)]) 0 0]
% v = [ 0  norm([A_(2,2) A_1(2,3)* [ 1 0 0 ]]
% zeros(1,1) =0 

v=[zeros(1,1); A_1(2:4,2)-norm(A_1(2:4,2))*I(2:4,2)]
H_2 = householder(v)

A_2=H_2*A_1

% v = [ 0  0 norm([A_2(3,2) A_1(4,3)])* [ 1 0 ]]
% I(3:4,3)=[1 0]'
%zeros(2,1)=[0  0]'

v= [ zeros(2,1) ; A_2(3:4,3)-norm(A_2(3:4,3)) *I(3:4,3)]

H_3 = householder(v)
A_3 = H_3*A_2

% This is the Q part
Q= H_1*H_2*H_3
 
% This is the R part
R=A_3

% This check the difference
D_1=Q*R-A

% This check if Q is unitary
D_2 = eye(4)-Q'*Q

echo off