function [Q R] = BlockHess (A, r)
% Block Hessenberg Reduction
% p - number of blocks
% r - number of extra diagonals
% pr x r block structure
[numRows, numCols]  = size(A)
numParts = max(2, floor(numCols/r));

a = floor(numCols/numParts);
partition = ones(1, numParts)*a;

% Split matrix rows into partition, storing result in a cell array
R = mat2cell(A, numRows, partition)
Q=eye(numRows);
V=eye(numRows);
for k = 1:numParts
    for j=1:r
        s=j+(k-1)*r;
        v=house(R{k}(s:numRows, j))
        V(s:numRows,j)=v;
        R{k}(s:numRows, j:r) = R{k}(s:numRows, j:r)-2*v*v'*R{k}(s:numRows, j:r);
    end
    Y=V(:,1);
    W=-2*V(:,1);
    for j=2:r
        z=-2*(eye(length(W))+W*Y')*V(:,j);
        W=[W z];
        Y=[Y V(:,j)]
    end
    if k<numParts
        R_tmp=horzcat(R{k+1:numParts});
        R_tmp=(eye(length(Y))+Y*W')*R_tmp;
        R_tmp=horzcat(R{1:k}, R_tmp);
        R = mat2cell(R_tmp, numRows, partition);
        Q=Q*(eye(length(Y))+W*Y');
    end
end

R=cell2mat(R);
end