
function [H] = BlockHess1 (A, r, p) 
% Block Hessenberg Reduction
% p - number of blocks
% r - number of extra diagonals
% pr x r block structure
n=length(A);
Q=eye(n);
k = n/r; % blocks in the first block column. 
for j = 1:r:n-r
k = max(k - 1, 2); % forces at least 2 blocks
l = ceil((k-1)/(p-1));% steps required for the reduction.
i = n;
for step = 1:l
nb = min(p*r, i-j-r+1);
% Annihilation of p r × r blocks in block column j of A. 
[q, A(i-nb+1:i,j:j+r-1)] = qr(A(i-nb+1:i,j:j+r-1)); 
A(i-nb+1:i,j+r:n) = q'*A(i-nb+1:i,j+r:n)*q; 
Q(:,i-nb+1:i) = Q(:,i-nb+1:i)*q; 
i = i - nb + r; %next block.
end
end
H=A;
end