if any(eq(ub, -inf)) 
   error('-Inf detected in upper bound: upper bounds must be > -Inf.');
elseif any(eq(lb,inf))
   error('+Inf detected in lower bound: lower bounds must be < Inf.');
end

x = xin;
function [X,lambda,exitflag,output,how]=qpsub(H,f,A,B,lb,ub,X,neqcstr,verbosity,caller,ncstr,numberOfVariables,options,defaultopt)


% Define constant strings
NewtonStep = 'Newton';
SteepDescent = 'steepest descent';
Conls = 'lsqlin';
Lp = 'linprog';
Qp = 'quadprog';
Qpsub = 'qpsub';
Nlconst = 'nlconst';
how = 'ok'; 

% Override quadprog/linprog MaxIter default limit which is
% really just for largescale, not active set methods.
defaultopt.MaxIter = Inf;  

exitflag = 1;
output = [];
iterations = 0;
if nargin < 13
   options = [];
end

lb=lb(:); ub = ub(:);

msg = nargchk(12,12,nargin);
if isempty(verbosity), verbosity = 1; end
if isempty(neqcstr), neqcstr = 0; end

LLS = 0;
if strcmp(caller, Conls)
   LLS = 1;
   [rowH,colH]=size(H);
   numberOfVariables = colH;
end
if strcmp(caller, Qpsub)
   normalize = -1;
else
   normalize = 1;
end

simplex_iter = 0;
if  norm(H,'inf')==0 | isempty(H), is_qp=0; else, is_qp=1; end



if LLS==1
   is_qp=0;
end

normf = 1;
if normalize > 0
   % Check for lp
   if ~is_qp & ~LLS
      normf = norm(f);
      if normf > 0
         f = f./normf;
      end
   end
end

% Handle bounds as linear constraints
arglb = ~eq(lb,-inf);
lenlb=length(lb); % maybe less than numberOfVariables due to old code
if nnz(arglb) > 0     
   lbmatrix = -eye(lenlb,numberOfVariables);
   
   A=[A; lbmatrix(arglb,1:numberOfVariables)]; % select non-Inf bounds
   B=[B;-lb(arglb)];
end


argub = ~eq(ub,inf);
lenub=length(ub);
if nnz(argub) > 0
   ubmatrix = eye(lenub,numberOfVariables);
   A=[A; ubmatrix(argub,1:numberOfVariables)];
   B=[B; ub(argub)];
end 
ncstr=ncstr + nnz(arglb) + nnz(argub);

% Figure out max iteration count
% For linprog/quadprog/lsqlin/qpsub problems, use 'MaxIter' for this.
% For nlconst (fmincon, etc) problems, use 'MaxSQPIter' for this.
if isequal(caller,Nlconst)
	maxiter = optimget(options,'MaxSQPIter',defaultopt,'fast');
else
	maxiter = optimget(options,'MaxIter',defaultopt,'fast'); 
end
% Used for determining threshold for whether a direction will violate
% a constraint.
normA = ones(ncstr,1);
if normalize > 0 
   for i=1:ncstr
      n = norm(A(i,:));
      if (n ~= 0)
         A(i,:) = A(i,:)/n;
         B(i) = B(i)/n;
         normA(i,1) = n;
      end
   end
else 
   normA = ones(ncstr,1);
end
errnorm = 0.01*sqrt(eps); 

tolDep = 100*numberOfVariables*eps;      
lambda=zeros(ncstr,1);
aix=lambda;
ACTCNT=0;
ACTSET=[];
ACTIND=0;
CIND=1;
eqix = 1:neqcstr; 

%------------EQUALITY CONSTRAINTS---------------------------
Q = zeros(numberOfVariables,numberOfVariables);
R = []; 
indepInd = 1:ncstr; 

if neqcstr>0
   % call equality constraint solver
   [Q,R,A,B,CIND,X,Z,actlambda,how,...
         ACTSET,ACTIND,ACTCNT,aix,eqix,neqcstr,ncstr,remove,exitflag]= ...
      eqnsolv(A,B,eqix,neqcstr,ncstr,numberOfVariables,LLS,H,X,f,normf,normA,verbosity, ...
      aix,how,exitflag);   
   
   if ~isempty(remove)
      indepInd(remove)=[];
      normA = normA(indepInd);
   end
   
   if ACTCNT >= numberOfVariables - 1  
      simplex_iter = 1; 
   end
   [m,n]=size(ACTSET);
   
   if strcmp(how,'infeasible')
      % Equalities are inconsistent, so X and lambda have no valid values
      % Return original X and zeros for lambda.
      output.iterations = iterations;
      return
   end
   
   err = 0;
   if neqcstr > numberOfVariables
      err = max(abs(A(eqix,:)*X-B(eqix)));
      if (err > 1e-8)  % Equalities not met
         how='infeasible';
         % was exitflag = 7; 
         exitflag = -1;
         if verbosity > 0 
            disp('Exiting: The equality constraints are overly stringent;')
            disp('         there is no feasible solution.')
         end
         % Equalities are inconsistent, X and lambda have no valid values
         % Return original X and zeros for lambda.
         output.iterations = iterations;
         return
      else % Check inequalities
         if (max(A*X-B) > 1e-8)
            how = 'infeasible';
            % was exitflag = 8; 
            exitflag = -1;
            if verbosity > 0
               disp('Exiting: The constraints or bounds are overly stringent;')
               disp('         there is no feasible solution.')
               disp('         Equality constraints have been met.')
            end
         end
      end
      if is_qp
         actlambda = -R\(Q'*(H*X+f));
      elseif LLS
         actlambda = -R\(Q'*(H'*(H*X-f)));
      else
         actlambda = -R\(Q'*f);
      end
      lambda(indepInd(eqix)) = normf * (actlambda ./normA(eqix));
      output.iterations = iterations;
      return
   end
   if isempty(Z)
      if is_qp
         actlambda = -R\(Q'*(H*X+f));
      elseif LLS
         actlambda = -R\(Q'*(H'*(H*X-f)));
      else
         actlambda = -R\(Q'*f);
      end
      lambda(indepInd(eqix)) = normf * (actlambda./normA(eqix));
      if (max(A*X-B) > 1e-8)
         how = 'infeasible';
         % was exitflag = 8; 
         exitflag = -1;
         if verbosity > 0
            disp('Exiting: The constraints or bounds are overly stringent;')
            disp('         there is no feasible solution.')
            disp('         Equality constraints have been met.')
         end
      end
      output.iterations = iterations;
      return
   end
   
   
   % Check whether in Phase 1 of feasibility point finding. 
   if (verbosity == -2)
      cstr = A*X-B; 
      mc=max(cstr(neqcstr+1:ncstr));
      if (mc > 0)
         X(numberOfVariables) = mc + 1;
      end
   end
else
   Z=1;
end

% Find Initial Feasible Solution
cstr = A*X-B;
mc=max(cstr(neqcstr+1:ncstr));
if mc>eps
   A2=[[A;zeros(1,numberOfVariables)],[zeros(neqcstr,1);-ones(ncstr+1-neqcstr,1)]];
   quiet = -2;
   options = struct('MaxIter',Inf);
   defaultopt = options;
   [XS,lambdaS,exitflagS,outputS] = qpsub([],[zeros(numberOfVariables,1);1],A2,[B;1e-5], ...
      [],[],[X;mc+1],neqcstr,quiet,Qpsub,size(A2,1),numberOfVariables+1,options,defaultopt);
   X=XS(1:numberOfVariables);
   cstr=A*X-B;
   if XS(numberOfVariables+1)>eps 
      if XS(numberOfVariables+1)>1e-8 
         how='infeasible';
         % was exitflag = 4; 
         exitflag = -1;
         if verbosity > 0
            disp('Exiting: The constraints are overly stringent;')
            disp('         no feasible starting point found.')
         end
      else
         how = 'overly constrained';
         % was exitflag = 3; 
         exitflag = -1;
         if verbosity > 0
            disp('Exiting: The constraints are overly stringent;')
            disp(' initial feasible point found violates constraints ')
            disp(' by more than eps.');
         end
      end
      lambda(indepInd) = normf * (lambdaS((1:ncstr)')./normA);
      output.iterations = iterations;
      return
   end
end

if (is_qp)
   gf=H*X+f;
   %  SD=-Z*((Z'*H*Z)\(Z'*gf));
   [SD, dirType] = compdir(Z,H,gf,numberOfVariables,f);
   
   % Check for -ve definite problems:
   %  if SD'*gf>0, is_qp = 0; SD=-SD; end
elseif (LLS)
   HXf=H*X-f;
   gf=H'*(HXf);
   HZ= H*Z;
   [mm,nn]=size(HZ);
   if mm >= nn
      %   SD =-Z*((HZ'*HZ)\(Z'*gf));
      [QHZ, RHZ] =  qr(HZ,0);
      Pd = QHZ'*HXf;
      % Now need to check which is dependent
      if min(size(RHZ))==1 % Make sure RHZ isn't a vector
         depInd = find( abs(RHZ(1,1)) < tolDep);
      else
         depInd = find( abs(diag(RHZ)) < tolDep );
      end  
   end
   if mm >= nn & isempty(depInd) % Newton step
      SD = - Z*(RHZ(1:nn, 1:nn) \ Pd(1:nn,:));
      dirType = NewtonStep;
   else % steepest descent direction
      SD = -Z*(Z'*gf);
      dirType = SteepDescent;
   end
else % lp
   gf = f;
   SD=-Z*Z'*gf;
   dirType = SteepDescent; 
   if norm(SD) < 1e-10 & neqcstr
      % This happens when equality constraint is perpendicular
      % to objective function f.x.
      actlambda = -R\(Q'*(gf));
      lambda(indepInd(eqix)) = normf * (actlambda ./ normA(eqix));
      output.iterations = iterations;
      return;
   end
end

oldind = 0; 

% The maximum number of iterations for a simplex type method is when ncstr >=n:
% maxiters = prod(1:ncstr)/(prod(1:numberOfVariables)*prod(1:max(1,ncstr-numberOfVariables)));

%--------------Main Routine-------------------
while iterations < maxiter
   iterations = iterations + 1;
   if isinf(verbosity)
      curr_out = sprintf('Iter: %5.0f, Active: %5.0f, step: %s, proc: %s',iterations,ACTCNT,dirType,how);
      disp(curr_out); 
   end
   
   % Find distance we can move in search direction SD before a 
   % constraint is violated.
   % Gradient with respect to search direction.
   GSD=A*SD;
   
   % Note: we consider only constraints whose gradients are greater
   % than some threshold. If we considered all gradients greater than 
   % zero then it might be possible to add a constraint which would lead to
   % a singular (rank deficient) working set. The gradient (GSD) of such
   % a constraint in the direction of search would be very close to zero.
   indf = find((GSD > errnorm * norm(SD))  &  ~aix);
   
   if isempty(indf) % No constraints to hit
      STEPMIN=1e16;
      dist=[]; ind2=[]; ind=[];
   else % Find distance to the nearest constraint
      dist = abs(cstr(indf)./GSD(indf));
      [STEPMIN,ind2] =  min(dist);
      ind2 = find(dist == STEPMIN);
      % Bland's rule for anti-cycling: if there is more than one 
      % blocking constraint then add the one with the smallest index.
      ind=indf(min(ind2));
      % Non-cycling rule:
      % ind = indf(ind2(1));
   end
   %-----Update X-------------
   
   % Assume we do not delete a constraint
   delete_constr = 0;   
   
   if ~isempty(indf)& isfinite(STEPMIN) % Hit a constraint
      if strcmp(dirType, NewtonStep)
         % Newton step and hit a constraint: LLS or is_qp
         if STEPMIN > 1  % Overstepped minimum; reset STEPMIN
            STEPMIN = 1;
            delete_constr = 1;
         end
         X = X+STEPMIN*SD;
      else
         % Not a Newton step and hit a constraint: is_qp or LLS or maybe lp
         X = X+STEPMIN*SD;          
      end              
   else %  isempty(indf) | ~isfinite(STEPMIN)
      % did not hit a constraint
      if strcmp(dirType, NewtonStep)
         % Newton step and no constraint hit: LLS or maybe is_qp
         STEPMIN = 1;   % Exact distance to the solution. Now delete constr.
         X = X + SD;
         delete_constr = 1;
      else % Not a Newton step: is_qp or lp or LLS
         if is_qp
            % Is it semi-def, neg-def or indef?
            eigoptions.disp = 0;
            ZHZ = Z'*H*Z;
            if numberOfVariables < 400 % only use EIGS on large problems
               [VV,DD] = eig(ZHZ);
               [smallRealEig, eigind] = min(diag(DD));
               ev = VV(:,eigind(1));
            else
               [ev,smallRealEig,flag] = eigs(ZHZ,1,'sr',eigoptions);
               if flag  % Call to eigs failed
                  [VV,DD] = eig(ZHZ);
                  [smallRealEig, eigind] = min(diag(DD));
                  ev = VV(:,eigind(1));
               end
            end
            
         else % define smallRealEig for LLS
            smallRealEig=0;
         end
         
         if (~is_qp & ~LLS) | (smallRealEig < -100*eps) % LP or neg def: not LLS
            % neg def -- unbounded
            if norm(SD) > errnorm
               if normalize < 0
                  STEPMIN=abs((X(numberOfVariables)+1e-5)/(SD(numberOfVariables)+eps));
               else 
                  STEPMIN = 1e16;
               end
               X=X+STEPMIN*SD;
               how='unbounded'; 
               % was exitflag = 5; 
               exitflag = -1;