function [x2, duals_rlx, idx_workc, idx_bindc] = LPrelax(a, f, b, nequs, vlb, vub, idx_workc, mpopt)

%   MATPOWER
%   $Id: LPrelax.m,v 1.12 2009/12/11 16:07:13 ray Exp $
%   by Deqiang (David) Gan, PSERC Cornell & Zhejiang University
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

%% options
alg     = mpopt(11);

if alg == 320		%% dense LP
    idx_workc = find(b < 0.001);
end


converged = 0;
while converged == 0

    atemp = a(idx_workc, :);
    btemp = b(idx_workc);

    [x2, duals] = mp_lp(f, atemp, btemp, vlb, vub, [], nequs, -1);

    diffs = b - a * x2;                 % diffs should be normalized by what means? under development
    idx_bindc = find(diffs < 1.0e-8);
    
    if isempty(find(diffs < -1.0e-8))
        converged = 1;
    else
        flag = zeros(length(b), 1);         % set up flag from scratch
        flag(idx_workc) = ones(length(idx_workc), 1);   % enforce historical working constraints


        idx_add = find(diffs < 0.001);
        flag(idx_add) = ones(length(idx_add), 1);   % enforce violating constraints



        flag(1:nequs) = ones(nequs, 1);         % enforce original equality constraints
        idx_workc_new = find(flag);

        if length(idx_workc) == length(idx_workc_new)   % safeguard step
            if isempty(find(idx_workc - idx_workc_new))
                converged = 1;
            end
        end
        idx_workc = idx_workc_new;

    end

end

duals_rlx = zeros(length(b), 1);
duals_rlx(idx_workc) = duals(1:length(btemp));

