function [Aeq, Beq] = normab(m,n)
% normalization constraints, to make sure that the m by n matrix
% is column-normalized
rind = 0;
Aeq = zeros(0,m*n);
for c=1:n
  rind = rind+1;
  for r=1:m
    ind = coord2ind([r c],[m n]);
    Aeq(rind,ind) = 1;
  end
end
Beq = ones(rind,1);
return