function [x,T1,B1,A1,flg] = QPsolve(C,p,A,b)
%
% This is a combined Phase I/Phase II calculation for solving QP
% based on Wolfe's method. 
% We assume non-degeneracy in the entire computation.
% Written by Ming Gu for Math 170
% April 2007
%
% 
% Prepare for initial Phase I computation.
%
flg       = 0;
[m,n]     = size(A);
mask      = find(b<0);
A(mask,:) = - A(mask,:);
b(mask)   = - b(mask);
B         = (n+1:n+m)';
x         = [zeros(n,1);b];
c1        = [zeros(n,1);ones(m,1)];
y         = ones(m,1);
T         = [b eye(m); c1'*x,y'];
A1        = [A eye(m)];
%
% Starting Simplex Method
%
f = ['Starting Initial Phase I Simplex Iteration... '];
disp(f);
[ITER1,B,T,flg1] = ILPRevisedSimplex(A1,B,c1,T);
f = ['Phase I took ',num2str(ITER1),' Iterations. '];
disp(f);
if (flg1 ~= 0)
    f = ['Initial Phase I Failed, flg =  ',num2str(flg1)];
    disp(f);
    flg = flg1;
    return;
end
if (max(B)>n)
    f = ['Initial Phase I Detected Primal Infeasibility, max(B) = ',num2str(max(B))];
    disp(f);
    flg = 3;
    return;
end
f = ['Starting the Wolfe Algorithm ... '];
disp(f);
x   = zeros(n,1);
x(B)= A(:,B)\b;
z  = -(C*x+p);
D  = sign(z);
A1 =[[A, zeros(m,2*(n+m))]; [C, A',-A',-eye(n),diag(D)]];
m1 = m   + n;
n1 = 3*n + 2* m;
T1 = zeros(m1+1,m1+1);
T1(1:m1,1) = [x(B);abs(z)];
T1(1:m1,2:m1+1) = [[T(1:m,2:m+1), zeros(m,n)]; [-C(:,B)* T(1:m,2:m+1), diag(D)]];
for k=m+1:m+n
   T1(k,2:m+1) = D(k-m)* T1(k,2:m+1);
end
c1 = [zeros(2*(n+m),1);ones(n,1)];
B1 = [B(:);2*(n+m)+(1:n)'];
y1 = T1(1:m1,2:m1+1)'*c1(B1); 
x1 = zeros(n1,1);
x1(B1) =[x(B);abs(z)];
T1(end,:) = [c1'*x1,y1'];
XV = [n+2*m+(1:n)';zeros(2*m,1);(1:n)';zeros(n,1)];
[ITER,B1,T1,flg] = QPRevisedSimplex(A1,B1,c1,T1,XV);
f = ['Wolfe Algorithm took ',num2str(ITER),' Iterations. '];
disp(f);
if (flg ~= 0)
    f = ['Wolfe Algorithm Failed, flg =  ',num2str(flg)];
    disp(f);
    return;
end
x1     = zeros(n1,1);
x1(B1) = A1(:,B1)\[b;-p];
x      = x1(1:n);
u      = x1(n+1:n+m)-x1(n+m+1:n+2*m);
v      = x1(n+2*m+1:2*n+2*m);
z      = x1(2*n+2*m+1:3*n+2*m);
f = ['Wolfe Algorithm Successful'];
disp([min(x) norm(A*x-b) min(v) max(abs(x.*v)) min(z)  sum(z)]);
disp(f);
return;