function [w,d] = lpsubprob(S, prtlevel)
%
%  Gradient Sampling LP subproblem:  
%    max t                      min u
%    s.t.  S'S w >= t e         s.t  S'S y <= u e
%           e'w  = 1                  e'y = 1
%           w >= 0                    y >= 0
%  Return w and d = -Sw.
%  The LP is solved by a call to linprog which will invoke either the 
%  Matlab Optimization Toolbox or MOSEK software, depending on which is
%  installed and, if both are, their relative order in the Matlab Path.  
%  Generally, MOSEK is far preferable.  To select it, use "addpath" to
%  add it to the front of the path.  To see which is in use, use "which".
%  
%  Written by M. Overton (overton@cs.nyu.edu), last revised November 2004 
%
[m,N] = size(S);  % N gradients of length m
H = S'*S;
vlb = [zeros(N,1); -inf];   % lower bounds on w, t
vub = [];
c = [zeros(N,1);
     -1 ];
A = [-H ones(N,1)];   % N inequality constraints
blower = -inf*ones(N,1);
bupper = zeros(N,1);
Aeq = [ones(1,N) 0];  % 1 equality constraint
beq = 1;
lp_options = optimset('Display','off','TolX', 1e-12, 'TolFun', 1e-12, 'Simplex', 'on');
% solve LP by either Optimization Toolbox or, preferably, MOSEK
w = linprog(c, A, bupper, Aeq, beq, vlb, vub,[],lp_options);  
if isempty(w)
    error('w is empty: MOSEK license problem?')
end
checktol = 1e-7*max(1,norm(H, inf));
if prtlevel > 0 & max(A*w)> checktol
    fprintf('LP inequalities not satisfied within desired accuracy\n')
end
w = w(1:N); % discard extra variable t
d = -S*w;
if prtlevel > 0 & (min(w) < -checktol | abs(sum(w) - 1) > checktol)
    fprintf('computed w does not satisfy requirements\n')
end