function Aeq = HBlkDiag_congress(ny,nu,bsy,bsu)
cbsu = cumsum(bsu);
cbsy = cumsum(bsy);
Hmin = zeros(nu,ny);
nsu = size(bsu,2);
nsy = size(bsy,2);
for k = 1:nu
    Hmin(k,1:ny) = [ny*(k-1)+1:ny*k];
end
zy = reshape(Hmin,1,nu*ny);

if (size(bsu,2)~= size(bsy,2))
    display('The number of blocks on U and Y should be same');
end

if (sum(bsu)~= nu)
    display('The number of inputs should match the sum of number of inputs in blocks')
end

if (sum(bsy)~= ny)
    display('The number of measurements should match the sum of number of measurements in blocks')
end

if (nsu~= nsy)
    display('The number of measurement blocks should match the number of input blocks')
end

Hf = [];
for zu = 1:size(bsu,2)
    if (zu == 1)
        Hm = [zeros(cbsu(1),cbsy(1)) Hmin(zu:cbsu(zu),cbsy(zu)+1:cbsy(zu+1))];
    end
    if (zu == size(bsu,2))
        Hm = [Hmin(cbsu(zu-1)+1:cbsu(zu),1:cbsy(zu))];
    end
    Hf = [Hf;Hm];
end

HBlk = reshape(Hf',nu*ny,1);
ind = find(HBlk==0);

for zzy = 1:length(ind)
    i = ind(zzy);
    P(zzy,:) = [zeros(1,i-1) ones(1,1) zeros(1,nu*ny-i)];
end
Aeq = P;

% syms x x1
% x = [x1];
% for l = 2:nu*ny
%     x = [x;sprintf('x%d',l)];
% end
% H = conj(reshape(x,ny,nu)');
% Hchk = Aeq * x;
% for i = 1:length(Hchk)
%     [row,col] = find(H==(Hchk(i)));
%     H(row,col) = 0;
% end
% H


